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)`

A[3] * t**3 + A[2] * t**2 + A[1] * t + A[0] == 0and 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