Permalink
Browse files

Add MALeaf - this combines MALeaf with MAleafViens

and supporting functions in MAEquations.h
as well as MARotateZ to add variation
  • Loading branch information...
sambler sambler
sambler authored and sambler committed Jan 25, 2013
1 parent b5c274c commit 74f2b16668316bf7c9a44485e6e5fa791bae4c84
View
@@ -0,0 +1,211 @@
+/*
+ * MAEquations.h by Michel J. Anders (c)2013
+ * from https://github.com/sambler/osl-shaders
+ *
+ * license: cc-by-sa
+ *
+ * original script from -
+ * http://blenderthings.blogspot.nl/p/a-small-osl-libray.html
+ *
+ */
+
+// cubic roots adapted for OSL from http://van-der-waals.pc.uni-koeln.de/quartic/quartic.html
+// so now we have a translation from Fortran -> C -> OSL. It doesn't look that well, but it works :-)
+
+float CBRT(float Z) { return abs(pow(abs(Z),1.0/3.0)) * sign(Z); }
+
+/*-------------------- Global Function Description Block ----------------------
+*
+* ***CUBIC************************************************08.11.1986
+* Solution of a cubic equation
+* Equations of lesser degree are solved by the appropriate formulas.
+* The solutions are arranged in ascending order.
+* NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
+* ******************************************************************
+* A(0:3) (i) vector containing the polynomial coefficients
+* X(1:L) (o) results
+* L (o) number of valid solutions (beginning with X(1))
+* ==================================================================
+* 17-Oct-2004 / Raoul Rausch
+* Conversion from Fortran to C
+* 09-Jan-2012 / Michel Anders (varkenvarken)
+* Conversion from C to OSL
+*
+*-----------------------------------------------------------------------------
+*/
+
+int cubic(float A[4], float X[3], int L)
+{
+ float PI = 3.1415926535897932;
+ float THIRD = 1./3.;
+ float U[3],W, P, Q, DIS, PHI;
+ int i;
+
+ // ====determine the degree of the polynomial ====
+
+ if (A[3] != 0.0)
+ {
+ //cubic problem
+ W = A[2]/A[3]*THIRD;
+ P = pow((A[1]/A[3]*THIRD - pow(W,2)),3);
+ Q = -.5*(2.0*pow(W,3)-(A[1]*W-A[0])/A[3] );
+ DIS = pow(Q,2)+P;
+ if ( DIS < 0.0 )
+ {
+ //three real solutions!
+ //Confine the argument of ACOS to the interval [-1;1]!
+ PHI = acos(min(1.0,max(-1.0,Q/sqrt(-P))));
+ P=2.0*pow((-P),(5.e-1*THIRD));
+ for (i=0;i<3;i++)
+ U[i] = P*cos((PHI+2*((float)i)*PI)*THIRD)-W;
+ X[0] = min(U[0], min(U[1], U[2]));
+ X[1] = max(min(U[0], U[1]),max( min(U[0], U[2]), min(U[1], U[2])));
+ X[2] = max(U[0], max(U[1], U[2]));
+ L = 3;
+ }
+ else
+ {
+ // only one real solution!
+ DIS = sqrt(DIS);
+ X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
+ L=1;
+ }
+ }
+ else if (A[2] != 0.0)
+ {
+ // quadratic problem
+ P = 0.5*A[1]/A[2];
+ DIS = pow(P,2)-A[0]/A[2];
+ if (DIS > 0.0)
+ {
+ // 2 real solutions
+ X[0] = -P - sqrt(DIS);
+ X[1] = -P + sqrt(DIS);
+ L=2;
+ }
+ else
+ {
+ // no real solution
+ L=0;
+ }
+ }
+ else if (A[1] != 0.0)
+ {
+ //linear equation
+ X[0] =A[0]/A[1];
+ L=1;
+ }
+ else
+ {
+ //no equation
+ L=0;
+ }
+/*
+* ==== perform one step of a newton iteration in order to minimize
+* round-off errors ====
+*/
+ for (i=0;i < L;i++)
+ {
+ X[i] = X[i] - (A[0]+X[i]*(A[1]+X[i]*(A[2]+X[i]*A[3])))
+ /(A[1]+X[i]*(2.0*A[2]+X[i]*3.0*A[3]));
+ }
+
+ return 0;
+}
+
+
+point cubicspline(float t, point P0, point P1, point P2, point P3)
+{
+ return (1-t)*(1-t)*(1-t)*P0
+ + 3*(1-t)*(1-t)*t*P1
+ + 3*(1-t)*t*t*P2
+ + t*t*t*P3;
+
+ // how to group all factors of the powers of t
+ //
+ // (1-2t+t2-t+2t2+t3)P0 => (1-3t+3t2+t3)P0
+ // (3t-6t2+3t3)P1
+ // (3t2-3t3)P2
+ // t3P3
+ //
+ // 1 * ( P0 )
+ // t * ( 3P0+3P1 )
+ // t2* ( 3P0-6P1+3P2 )
+ // t3* ( P0+3P1-3P2+P3)
+}
+
+// calculate the closest distance from Pos to a quadratic bezier curve
+// the curve is defined by the 3 points P0, P1 and P2
+
+int splinedist(point p0, point p1, point p2, point Pos, float d, float tc)
+{
+
+ point P0 = p0;
+ point P1 = p1;
+ point P2 = p2;
+
+ // following definitions are for the four polynomic coefficients for the well known
+ // equation dB/dt . (Pos-B) (i.e. the inproduct of the tangent to the bezier and the
+ // difference vector from the point under considertion to the Bezier curve.
+ // If the difference vector is perpendicular to the tangent we have found a closest point
+ // on the Bezier curve.
+ // The stuff below is generated by a script and no effort is spent on collecting factors.
+ // We let the OSL compiler worry about that :-)
+
+ float t0 = -2*P0[0]*P0[0]+-2*P1[0]*Pos[0]+2*P0[0]*P1[0]+2*P0[0]*Pos[0]
+ +-2*P0[1]*P0[1]+-2*P1[1]*Pos[1]+2*P0[1]*P1[1]+2*P0[1]*Pos[1]
+ +-2*P0[2]*P0[2]+-2*P1[2]*Pos[2]+2*P0[2]*P1[2]+2*P0[2]*Pos[2];
+
+ float t1 = -4*P0[0]*P1[0]+-4*P0[0]*P1[0]+-4*P0[0]*P1[0]+-2*P0[0]*Pos[0]
+ +-2*P2[0]*Pos[0]+2*P0[0]*P0[0]+2*P0[0]*P2[0]+4*P0[0]*P0[0]
+ +4*P1[0]*P1[0]+4*P1[0]*Pos[0]+-4*P0[1]*P1[1]+-4*P0[1]*P1[1]
+ +-4*P0[1]*P1[1]+-2*P0[1]*Pos[1]+-2*P2[1]*Pos[1]+2*P0[1]*P0[1]
+ +2*P0[1]*P2[1]+4*P0[1]*P0[1]+4*P1[1]*P1[1]+4*P1[1]*Pos[1]
+ +-4*P0[2]*P1[2]+-4*P0[2]*P1[2]+-4*P0[2]*P1[2]+-2*P0[2]*Pos[2]
+ +-2*P2[2]*Pos[2]+2*P0[2]*P0[2]+2*P0[2]*P2[2]+4*P0[2]*P0[2]
+ +4*P1[2]*P1[2]+4*P1[2]*Pos[2];
+
+ float t2 = -8*P1[0]*P1[0]+-4*P0[0]*P0[0]+-4*P0[0]*P2[0]+-4*P1[0]*P1[0]
+ +-2*P0[0]*P0[0]+-2*P0[0]*P2[0]+2*P0[0]*P1[0]+2*P1[0]*P2[0]
+ +4*P0[0]*P1[0]+4*P0[0]*P1[0]+4*P1[0]*P2[0]+8*P0[0]*P1[0]
+ +-8*P1[1]*P1[1]+-4*P0[1]*P0[1]+-4*P0[1]*P2[1]+-4*P1[1]*P1[1]
+ +-2*P0[1]*P0[1]+-2*P0[1]*P2[1]+2*P0[1]*P1[1]+2*P1[1]*P2[1]
+ +4*P0[1]*P1[1]+4*P0[1]*P1[1]+4*P1[1]*P2[1]+8*P0[1]*P1[1]
+ +-8*P1[2]*P1[2]+-4*P0[2]*P0[2]+-4*P0[2]*P2[2]+-4*P1[2]*P1[2]
+ +-2*P0[2]*P0[2]+-2*P0[2]*P2[2]+2*P0[2]*P1[2]+2*P1[2]*P2[2]
+ +4*P0[2]*P1[2]+4*P0[2]*P1[2]+4*P1[2]*P2[2]+8*P0[2]*P1[2];
+
+ float t3 = -4*P0[0]*P1[0]+-4*P0[0]*P1[0]+-4*P1[0]*P2[0]+-4*P1[0]*P2[0]
+ +2*P0[0]*P0[0]+2*P0[0]*P2[0]+2*P0[0]*P2[0]+2*P2[0]*P2[0]
+ +8*P1[0]*P1[0]+-4*P0[1]*P1[1]+-4*P0[1]*P1[1]+-4*P1[1]*P2[1]
+ +-4*P1[1]*P2[1]+2*P0[1]*P0[1]+2*P0[1]*P2[1]+2*P0[1]*P2[1]
+ +2*P2[1]*P2[1]+8*P1[1]*P1[1]+-4*P0[2]*P1[2]+-4*P0[2]*P1[2]
+ +-4*P1[2]*P2[2]+-4*P1[2]*P2[2]+2*P0[2]*P0[2]+2*P0[2]*P2[2]
+ +2*P0[2]*P2[2]+2*P2[2]*P2[2]+8*P1[2]*P1[2];
+
+ float A[4] = {t0,t1,t2,t3};
+ float T[3] ;
+ int n ;
+ cubic(A,T,n);
+ d = 1e6;
+ // cubic() will return 0 , 1 or 3 values for t
+ // we are only interested in values that lie in the interval [0,1]
+ // for those we calculate the position on the curve and check whether
+ // we have found the shortest distance.
+ int found = 0;
+ while(n>0){
+ n--;
+ if(T[n]>=0 && T[n]<=1){
+ float t = T[n];
+ found = 1;
+ float dd = distance((1-t)*(1-t)*P0 + 2*(1-t)*t *P1 + t*t*P2, Pos);
+ if (dd < d) {
+ d = dd;
+ tc = t;
+ }
+ }
+ }
+ return found;
+}
+
+
View
@@ -0,0 +1,86 @@
+/*
+ * MALeaf.osl by Michel J. Anders (c)2013
+ * from https://github.com/sambler/osl-shaders
+ *
+ * license: cc-by-sa
+ *
+ * original script from -
+ * http://blenderthings.blogspot.com.au/2013/01/a-osl-leaf-shape-shader-for-cycles.html
+ *
+ */
+
+#include "stdosl.h"
+#include "MAEquations.h"
+
+shader leaf(
+ point Vector = P,
+ float BaseAngle = 107.0,
+ float BaseCurve = 0.7,
+ float TipAngle = 44.0,
+ float TipCurve = 1.2,
+ output float Leaf = 0 )
+{
+ // calculate the four control point of the cubic spline
+ float x1,y1,x2,y2;
+ sincos(radians(BaseAngle),y1,x1);
+ sincos(radians(TipAngle),y2,x2);
+ point P0 = point(0, 0, 0);
+ point P1 = point(x1, y1, 0)*BaseCurve;
+ point P2 = point(1-x2*TipCurve, y2*TipCurve, 0);
+ point P3 = point(1, 0, 0);
+
+ // to determin the y value(s) of the spline at the x position we
+ // are located, we want to solve spline(t) - x = 0
+ // we therefore gather all factors and solve the cubic equation
+ float tfactor[4] = { P0[0]-Vector[0],
+ 3*P0[0]+3*P1[0],
+ 3*P0[0]-6*P1[0]+3*P2[0],
+ P0[0]+3*P1[0]-3*P2[0]+P3[0] };
+ float t[3];
+ int nrealroots;
+ cubic(tfactor, t, nrealroots);
+
+ // at this point, the array t holds up to 3 real roots
+ // remove any real root that is not in range [0,1]
+ int i=0;
+ while(i < nrealroots){
+ if ((t[i] < 0) || (t[i] > 1)) {
+ int j=i;
+ while(j < (nrealroots-1)){
+ t[j]=t[j+1];
+ j++;
+ }
+ nrealroots--;
+ }
+ i++;
+ }
+
+ // note that a cubic funtion can have 3 real roots,
+ // but in this case we ignore such very warped curves
+ // TODO: w. 3 real roots w could set leaf = 1, if y < y0 OR y between y1,y2
+ // TODO: seration, possible by determining the closest
+ // distance (if inside leaf) to the spline and
+ // determining if w are within some periodic funtion f(t)
+
+ // we generate the shape mirrored about the x-axis
+ float y = Vector[1];
+ if(y<0) y = -y;
+
+ if(nrealroots > 0){
+ point Sy0 = cubicspline(t[0],P0,P1,P2,P3);
+ if(nrealroots > 1){
+ // if we have 2 roots we calculate and order the y values
+ // and check whether the current y values is between them
+ point Sy1 = cubicspline(t[1],P0,P1,P2,P3);
+ if ( Sy1[1] < Sy0[1] ){
+ if( (y > Sy1[1]) && (y < Sy0[1]) ) Leaf = 1;
+ }else{
+ if( (y > Sy0[1]) && (y < Sy1[1]) ) Leaf = 1;
+ }
+ }else{
+ // with a single value we check if we are below the y value
+ if( y < Sy0[1] ) Leaf = 1;
+ }
+ }
+}
+
@@ -0,0 +1,92 @@
+/*
+ * MALeafVeins.osl by Michel J. Anders (c)2013
+ * from https://github.com/sambler/osl-shaders
+ *
+ * license: cc-by-sa
+ *
+ * original script from -
+ * http://blenderthings.blogspot.com.au/2013/01/osl-leaf-veins-shader-for-cycles.html
+ *
+ */
+
+#include "stdosl.h"
+#include "MAEquations.h"
+
+shader arcuateveins(
+ point Vector = P,
+
+ float BaseAngle = 107.0,
+ float BaseCurve = 0.7,
+ float TipAngle = 44.0,
+ float TipCurve = 1.2,
+
+ int Veins = 7,
+ int Seed = 42,
+ float Variance = 0,
+ float InnerWidth = 0.05,
+ float NWidth = 0.25, // size of the reticulated area
+
+ float Curve = 0.5, // distribution of endpoints on edge
+ float Curve2 = 0.5, // distribution of controlpoints
+ float Spacing = 0.5, // distribution of starting points
+ float Up = 0.5,
+
+ output float Vein = 0,
+ output float Net = 0,
+ output float Fac = 0 )
+{
+
+ float delta = 1.0/((float)Veins+1);
+ float delta2= delta/2;
+ float delta4= delta/4;
+
+ // calculate the four control points of the cubic spline that defines the leaf edge
+ float x1,y1,x2,y2;
+ sincos(radians(BaseAngle),y1,x1);
+ sincos(radians(TipAngle),y2,x2);
+ point P0 = point(0, 0, 0);
+ point P1 = point(x1, y1, 0)*BaseCurve;
+ point P2 = point(1-x2*TipCurve, y2*TipCurve, 0);
+ point P3 = point(1, 0, 0);
+
+ point P0q = point(P0[0],P0[1]*Up,P0[2]);
+ point P1q = point(P1[0],P1[1]*Up,P1[2]);
+ point P2q = point(P2[0],P2[1]*Up,P2[2]);
+ point P3q = point(P3[0],P3[1]*Up,P3[2]);
+
+ int i;
+ for(i=0;i < Veins;i++){
+
+ // determine the starting points of the veins
+ float x = (i*delta+delta2*Variance*cellnoise(i+10+Seed))*Spacing;
+ float dx = (delta4*Variance*cellnoise(i+17+Seed))*Spacing;
+ point P0up = point(delta2+x+dx,0,0);
+ point P0down = point(delta2+x,0,0);
+ // determine the endpoints on the leaf edge
+ float t=(i*delta+delta2)*Curve+1-Curve;
+ point P2up = cubicspline(t,P0,P1,P2,P3);
+ point P2down = point(P2up[0],-P2up[1],P2up[2]);
+ // the veins are quadratic splines, so need one additional control point
+ t=(i*delta+delta2)*Curve2+1-Curve2;
+ point P1up = cubicspline(t,P0q,P1q,P2q,P3q);
+ point P1down = point(P1up[0],-P1up[1],P1up[2]);
+
+ float r;
+ int f = splinedist(P0up, P1up, P2up, Vector, r, t);
+ if ( f && (r < NWidth ) ) Net = 1 ;
+ if ( f && (r < InnerWidth * ( 1- t) * (1-Vector[0]) ) ) {
+ Vein = 1; Fac = sqrt(1-r/InnerWidth); break;
+ }
+ f = splinedist(P0down, P1down, P2down, Vector, r , t);
+ if ( f && (r < NWidth ) ) Net = 1 ;
+ if ( f && (r < InnerWidth * ( 1- t) * (1-Vector[0]) ) ) {
+ Vein = 1; Fac = sqrt(1-r/InnerWidth); break;
+ }
+ }
+
+ // the central vein
+ float d = distance(point(0,0,0),point(1,0,0),Vector);
+ if ( d < NWidth ) Net = 1 ;
+ if (d < (InnerWidth * (1-Vector[0])) ) { Vein = 1; Fac = sqrt(1-d/InnerWidth);}
+}
+
Oops, something went wrong.

0 comments on commit 74f2b16

Please sign in to comment.