wiggles.osl

// cubic, quartic and quitic 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 pretty, but it works well enough :-)

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
*
*   21-Nov-2013 / Michel Anders
*       add refine and sort parameters to make these actions optional
*-----------------------------------------------------------------------------
*/

int cubic(float A[4], float X[3], int L, int refine, int sort)
{
	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;
            if(sort){
			  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]));
			}else{
              X[0] = U[0];
              X[1] = U[1];
              X[2] = 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 ====
*/
    if(refine){
	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;
}

#define DOT(a,b) (a[0]*b[0]+a[1]*b[1])

// 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 d00=DOT(p0,p0), d2p=DOT(p2,Pos), d1p=DOT(p1,Pos), d0p=DOT(p0,Pos);
float d01=DOT(p0,p1), d02=DOT(p0,p2), d11=DOT(p1,p1);
float d12=DOT(p1,p2), d22=DOT(p2,p2);

float t0=-2*(d00+d1p-d01-d0p);

//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 t1=-12*d01-2*d0p+-2*d2p+6*d00 +2*d02+4*d11+4*d1p;

//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 t2=-8*d11-6*d00-4*d02-4*d11-2*d02+18*d01+6*d12;
float t2=-8*d11-6*(d00-d12)-4*(d02+d11)-2*d02+18*d01;

//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 t3=-8*d01-8*d12+2*d00+4*d02+2*d22+8*d11;
float t3=-8*(d01+d12-d11)+2*(d00+d22)+4*d02;


   float A[4] = {t0,t1,t2,t3};
   float T[3] ;
   int n ;
   cubic(A,T,n, 0,0);
   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;
}

#define SUB(a,b) vector(a[0]-b[0],a[1]-b[1],0)

// determine if point M is inside a rectangle with a margin
int in_rectangle(point M, point a, point b,
  vector u, float W, vector v, float linewidth){
  point A=a+linewidth*(-u-v);
  point B=b+linewidth*(u-v);
  point D=B+(W+2*linewidth)*v;
  vector AM=SUB(M,A);
  vector AD=SUB(D,A);
  vector AB=SUB(B,A);
  float dotamad=DOT(AM, AD);
  float dotadad=DOT(AD, AD);
  float dotamab=DOT(AM, AB);
  float dotabab=DOT(AB, AB);
  return (dotamad > 0 && dotamad < dotadad) && 
         (dotamab > 0 && dotamab < dotabab);
}

#define CELL noise("cell", cp, seed++)
#define CELL2 vector(CELL, CELL, 0)

shader wiggles(
  point Pos=P,
  float Scale=1,

  int Number=1,

  float Length=0.5,
  float LengthVar=0,
  float Kink=0,
  float Curl=0.2,
  float Wave=30,    // degrees
  int Steps=2,
  float StepsVar=0,

  float Width=0.02,
  float WidthVar=0,

  int Seed=0,
  
  output float Fac=0
){
  point p = Pos * Scale;
  p[2]=0;
  point ip= point(floor(p[0]),floor(p[1]),0);

  int nn=1+(int)ceil(Steps*Length);

  for(int xx=-nn; xx <= nn; xx++){
    for(int yy=-nn; yy <= nn; yy++){
      int seed=Seed;
      point cp = ip + vector(xx, yy, 0);
      for(int wiggle=0; wiggle<Number; wiggle++){
        vector start = cp + CELL2;
        start[2]=0;
        vector dir = CELL2 - 0.5;
        dir[2]=0;
        dir = normalize(dir);
          
        vector perp = vector(dir[1],-dir[0],0);
        float k=0.5 + Kink * (CELL-0.5);
        float c=Curl*(CELL-0.5);
        point p1=start+k*dir+c*perp;
        for(int step=0; step < Steps; step++){
          vector ldir = dir;
          ldir *= Length + LengthVar*CELL;
          point end=start+ldir;
          if(in_rectangle(p, start, end, dir, c, perp, Width+WidthVar)){
            float d,t;
            if(splinedist(start, p1, end, p, d, t)){
              float localwidth = Width+WidthVar*noise("uperlin",start,t);
              if(d < localwidth){
                Fac = (localwidth - d)/localwidth;
                return;
              }
            }
          }

          if(CELL < StepsVar){
            break;
          }else{
            p1 = end + (end - p1)*(1+noise("perlin",end)*Kink); 
            start = end;
            dir = rotate(dir, radians(Wave*noise("perlin", start)), vector(0,0,0), vector(0,0,1));
          }
        }
      }
    }
  }
}