A small OSL libray

The code below is meant to be saved as equations.h in a location that is looked in by OSLs #include directive. Currently is contains the following two functions (but more will be added as I develop other shaders):

  • cubic(float A[4], float X[3], int L)
  • point cubicspline(float t, point P0, point P1, point P2, point P3)
  • int splinedist(point p0, point p1, point p2, point Pos, float d, float tc)
The first one will find the real roots of the cubic equation
A[3] * t**3 + A[2] * t**2 + A[1] * t + A[0] == 0
and will return the roots in X and the number of (real) roots in L. I converted the code from http://van-der-waals.pc.uni-koeln.de/quartic/quartic.html to OSL, which is supprisingly simple. There is other cool stuff on that site, so check it out.

The second one will return a point on a cubic bezier spline that starts at P0, heads in the direction of P1 and ends at P3 coming from the direction of P2. (note that the spline() only works for uniformly distributed points and the documentation is a bit unclear so I thought I'd better implement one I understood)

The third one will return the closest distance to a cubic spline with control points P0, P1 and P2, returning the result in d and the corresponding t value in tc. The function returns 1 if a closest distance is found with a t value in the range [0,1], zero otherwise.

An example of how to include this code is shown below

#include "equation.h"

shader myshader( ... )
{
 ...
 point p = cubicspline(t,P0,P1,P2,P3);
 ...
}

equations.h

This is the actual library. Save it for example as scripts\addons\cycles\shaders\equations.h in your Blender installation directory or in the same directory as your main .osl file.

// 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
*
*-----------------------------------------------------------------------------
*/

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;
}

No comments:

Post a Comment