/* Include file: curve.inc Spline and curve interpolation for many variants. Based upon the works of Steve H. Noskowicz ! http://web.archive.org/web/20151002232205/http://home.comcast.net/~k9dci/site/?/page/Piecewise_Polynomial_Interpolation/ Get the book PDF and the appendix. Author : Ingo Janssen Version: 0.1 Alpha, 2019-06-16 initial Version: 0.2 Alpha, 2021-05-24 complete re-write from earlier attempt. cubic Lagrange still not right. Version: 0.3 Alpha, 2021-06-06 yet again a complete re-write. Simplification. cubic Lagrange still not right, left out. version: 0.4 Alpha, 2021-06-26 Fixed TCB spline bug More refactoring Added printing macro for curve object (3d vectors only) Added CIt curve, to kind of simulate the way splines work with the t-value */ #version 3.7; #ifndef(CURVE_INC_TEMP) #declare CURVE_INC_TEMP = version; #ifdef(View_POV_Include_Stack) #debug "including curve.inc\n" #end //enum CICurve items #declare _cu_items = array[4]{"Ctrl", "Point", "Curmult", "Matrix"} #declare _cu_ctrl = 0; #declare _cu_point = 1; #declare _cu_curmult = 2; #declare _cu_matrix = 3; //enum curve ctrl items #declare _cu_ctrl_items = array[5]{"DimSeg", "Step", "Len", "Segments", "CurrSeg"} #declare _cu_dimseg = 0; #declare _cu_step = 1; #declare _cu_len = 2; #declare _cu_segments = 3; #declare _cu_curseg = 4; //enum CITCurve items #declare _cu_tcurve = 0; #declare _cu_ccurve = 1; #macro _arrFill(Len, Val) /*:Create and fill an array of length Len with value Val*/ #local Arr = array[Len]; #for(i,0,Len-1) #local Arr[i] = Val; #end Arr #end #macro CIprint(CIarr) /*:Print the control values, the matrix and the point arrays Only works for CIcurves with 3D vectors! */ #debug concat(_cu_items[0]," : \n") #for(i, 0, dimension_size(CIarr[_cu_ctrl],1)-1) #debug concat(" ", _cu_ctrl_items[i]," : ", str(CIarr[_cu_ctrl][i],0,0),"\n") #end #debug concat(_cu_items[1]," : \n array{") #for(i, 0, dimension_size(CIarr[_cu_point],1)-1) #debug concat("<", vstr(3,CIarr[_cu_point][i],", ",0,-1),">, ") #end #debug "}\n" #debug concat(_cu_items[2]," : \n array{") #for(i, 0, dimension_size(CIarr[_cu_curmult],1)-1) #debug concat("<", vstr(3,CIarr[_cu_curmult][i],", ",0,-1),">, ") #end #debug "}\n" #debug concat(_cu_items[3]," : \n array{") #for(i, 0, dimension_size(CIarr[_cu_matrix],1)-1) #debug "\n array{" #for(j, 0, dimension_size(CIarr[_cu_matrix][i],1)-1) #debug concat("<", vstr(3,CIarr[_cu_matrix][i][j],", ",0,-1),">, ") #end #debug "}" #end #debug "\n }\n" #end #macro _CIinit(CIarr, DimSeg, Step, Matrix) #local Len = dimension_size(CIarr, 1); #local Lsd = Len - DimSeg; #if(Lsd < 0) #error "Too few points in curve array" #end #if(mod(Lsd, Step) = 0) #local Segments = div(Lsd, Step) + 1; #else #error "Not the right amount of points in curve array" #end #local CurrSeg = -1; #local CICurve = array[4]; #local CICurve[_cu_ctrl] = array[5]{DimSeg, Step, Len, Segments, CurrSeg}; #local CICurve[_cu_point] = CIarr; #local CICurve[_cu_curmult] = _arrFill(dimension_size(Matrix,1), <0,0,0>); #local CICurve[_cu_matrix] = Matrix; CICurve #end #macro _CIdivMatrix(Matrix, Div) #for(i, 0, dimension_size(Matrix,1)-1) #for(j, 0, dimension_size(Matrix[i],1)-1) #local Matrix[i][j] = Matrix[i][j] / Div; #end #end Matrix #end #macro CIt(tCurve, cCurve) /*: The tCurve modulates T for the cCurve. Use CITget te get a value. tCurve is interpolated [0 - 1], its output is consumed by cCurve. The modulated T value changes the speed/acceleration of the cCurve. parameter: tCurve (curve, float): any curve based on an array with floats [0-1] cCurve (curve, vector): any curve based on an array with vectors result: An array with the two curve arrays. */ #local CITCurve = array[2]{ tCurve, cCurve }; CITCurve #end #macro CIlinear(CIarr) /*:Linear interpolation*/ #local DimSeg = 2; #local Step = 1; #local Div = 1; #local Matrix = array[2]{ array[2]{-1, 1}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIsimple3(CIarr) /*:Simple Cubic, linear sections*/ #local DimSeg = 2; #local Step = 1; #local Div = 1; #local Matrix = array[4]{ array[2]{ 2,-2}, array[2]{-3, 3}, array[2]{ 0, 0}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIsimple3ESC(CIarr, ESC) /*:Simple Cubic, Endpoint Slope Control*/ #local DimSeg = 2; #local Step = 1; #local Div = 1; #local Matrix = array[4]{ array[2]{ 2 - (2 * ESC), (2 * ESC) - 2}, array[2]{ (3 * ESC) - 3, 3 - (3 * ESC)}, array[2]{-ESC, ESC}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIsimple3MSC(CIarr, MSC) /*:Simple Cubic, Midpoint Slope Control*/ #local DimSeg = 2; #local Step = 1; #local Div = 1; #local Matrix = array[4]{ array[2]{ (4 * MSC) - 4, 4 - (4 * MSC)}, array[2]{ 6 - (6 * MSC), (6 * MSC) - 6}, array[2]{ (2 * MSC) - 3, 3 - (2 * MSC)}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIsimple3IESC(CIarr, ESC0, ESC1) /*:Simple Cubic, Independent Endpoint Slope Control*/ #local DimSeg = 2; #local Step = 1; #local Div = 1; #local Matrix = array[4]{ array[2]{ 2 - ESC0 - ESC1, ESC0 + ESC1 - 2}, array[2]{ (2 * ESC0) + ESC1 - 3, 3 - (2 * ESC0) - ESC1}, array[2]{ -ESC0, ESC0}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIbezier2(CIarr) //:Quadratic Bezier curve #local DimSeg = 3; #local Step = 2; #local Div = 1; #local Matrix = array[3]{ array[3]{ 1,-2, 1}, array[2]{-2, 2}, array[1]{ 1}, } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIbezier2EVC(CIarr, EVC) //:Quadratic Bezier curve with End Velocity Control #local DimSeg = 3; #local Step = 2; #local Div = 1; #local Matrix = array[4]{ array[3]{-EVC, 0.0, EVC}, array[3]{ 1 + (2 * EVC), -2 - EVC, 1 - EVC}, array[2]{-2 - EVC, 2 + EVC}, array[1]{ 1}, } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIbezier3(CIarr) //:Cubic Bezier curve #local DimSeg = 4; #local Step = 3; #local Div = 1; #local Matrix = array[4]{ array[4]{-1, 3,-3, 1}, array[3]{ 3,-6, 3}, array[2]{-3, 3}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIbezier3EVC(CIarr, EVC) //:Cubic Bezier curve with End Velocity Control #local DimSeg = 4; #local Step = 3; #local Div = 1; #local Matrix = array[4]{ array[4]{-1.0 - EVC, 3.0 + EVC, -3.0 - EVC, 1.0 + EVC}, array[4]{ 3.0 + (2 * EVC), -6.0 - (2 * EVC), 3.0 + EVC, -EVC}, array[2]{-3.0 - EVC, 3.0 + EVC}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIbspline2(CIarr) //:Cubic b-spline #local DimSeg = 3; #local Step = 1; #local Div = 2; #local Matrix = array[3]{ array[3]{ 1,-2, 1}, array[2]{-2, 2}, array[2]{ 1, 1} } #local Matrix = _CIdivMatrix(Matrix, Div); _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIbspline3(CIarr) //:Cubic b-spline #local DimSeg = 4; #local Step = 1; #local Div = 6; #local Matrix = array[4]{ array[4]{-1, 3,-3, 1}, array[3]{ 3,-6, 3}, array[3]{-3, 0, 3}, array[3]{ 1, 4, 1} } #local Matrix = _CIdivMatrix(Matrix, Div); _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIgolden(CIarr) //:Golden curve #local DimSeg = 3; #local Step = 2; #local Div = 1; #local Matrix = array[3]{ array[3]{ 2,-4, 2}, array[3]{-3, 4,-1}, array[1]{1}, } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CI4p3(CIarr) //:Four Point Cubic #local DimSeg = 4; #local Step = 4; #local Div = 2; #local Matrix = array[4]{ array[4]{ -9, 27,-27, 9}, array[4]{ 18,-45, 36,-9}, array[4]{-11, 18, -9, 2}, array[1]{ 2} } #local Matrix = _CIdivMatrix(Matrix, Div); _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIhermite(CIarr) //:Cubic Hermite spline #local DimSeg = 4; #local Step = 3; #local Div = 1; #local Matrix = array[4]{ array[4]{ 2, 1, 1,-2}, array[4]{-3,-2,-1, 3}, array[2]{ 0, 1}, array[1]{ 1} } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIcatmull(CIarr) //:Catmull Rom spline #local DimSeg = 4; #local Step = 1; #local Div = 2; #local Matrix = array[4]{ array[4]{-1, 3,-3, 1}, array[4]{ 2,-5, 4,-1}, array[3]{-1, 0, 1}, array[2]{0,2} } #local Matrix = _CIdivMatrix(Matrix, Div); _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIlagrange2(CIarr, T2) //:Quadratic Lagrange curve #local DimSeg = 3; #local Step = 2; #local Div = 1; #local N1 = 1/T2; #local N2 = 1/(T2 * (T2 - 1)); #local N3 = 1/(1 - T2); #local Matrix = array[3]{ array[3]{ N1, N2, N3}, array[3]{-N1 - 1, -N2, -N3 * T2}, array[1]{ 1}, } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIbeta(CIarr, T, B) //:Beta-spline #local DimSeg = 4; #local Step = 1; #local Div = (T + (2 * pow(B, 3)) + (4 * pow(B,2)) + (4 * B) + 2); #local Matrix =array[4]{ array[4]{ -2 * pow(B, 3), 2 * (T + pow(B, 3) + pow(B, 2) + B), -2 * (T + pow(B, 2) + B + 1), 2 }, array[3]{ 6 * pow(B, 3), -3 * (T + (2 * pow(B, 3) + (2 * pow(B, 2)))), 3 * (T + (2 * pow(B, 2))) }, array[3]{ -6 * pow(B, 3), 6 * (pow(B, 3) - B), 6 * B }, array[3]{ 2 * pow(B, 3), T + (4 * (pow(B, 2) + B)), 2 } } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CItcb(CIarr, Tension, Cont, Bias) //:Kochanek Bartels aka TCB spline, Tension, Continuity, Bias #local DimSeg = 4; #local Step = 1; #local Div = 2; #local A = (1 - Tension) * (1 + Cont) * (1 + Bias); #local B = (1 - Tension) * (1 - Cont) * (1 - Bias); #local C = (1 - Tension) * (1 - Cont) * (1 + Bias); #local D = (1 - Tension) * (1 + Cont) * (1 - Bias); #local Matrix = array[4]{ array[4]{ -A, 4 + A - B - C, -4 + B + C - D, D }, array[4]{ 2 * A, -6 - (2 * A) + (2 * B) + C, 6 - (2 * B) - C + D, -D }, array[3]{ -A, A - B, B }, array[2]{0,2} } #local Matrix = _CIdivMatrix(Matrix, Div); _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIpalmer(CIarr, T, B) //:Palmer #local DimSeg = 3; #local Step = 2; #local Div = 1; #local N0 = 1 / B; #local N1 = 1 / (B * (1 - B)); #local N2 = 1 / (1 - B); #local Matrix = array[4]{ array[3]{ -2 * T, 0, 2 * T, }, array[3]{ N0 + (3 * T), -N1, N2 - (3 * T), }, array[3]{ -1 - N0 - T, N1, T - (N2 / N0), }, array[1]{ 1 } } _CIinit(CIarr, DimSeg, Step, Matrix) #end #macro CIget(CICurve, T) /*: Gets the point on the curve at T parameter: CICurve (array, vector): an array of vectors and other parameters as created via one of the spline initialisation macros. T (float, [0,1]): position at which to interpolate the curve. return: (vector) returns a point on the curve at T. use example: #declare TheCurve = CIlinear(PointArray); #declare P = CIget(TheCurve, 0.1); */ /* T goes from 0 to 1 for the whole curve. From current T find the proper curve segment. Tseg also goes from 0 to 1. Once the proper segment is found, get the proper Tseg. */ #local Len = CICurve[_cu_ctrl][_cu_len]; #local Segments = CICurve[_cu_ctrl][_cu_segments]; #local CImult = CICurve[_cu_curmult]; #if(T = 1) #local Segment = Segments - 1; #local Tseg = 1; #else #local Segment = floor(T * Segments); #local Tseg = (T * Segments) - Segment; #end #if(Segment != CICurve[_cu_ctrl][_cu_curseg]) #local CICurve[_cu_ctrl][_cu_curseg] = Segment; #local CurrPos = Segment * CICurve[_cu_ctrl][_cu_step]; #local CImatrix = CICurve[_cu_matrix]; #local MaxPos = Len - CICurve[_cu_ctrl][_cu_dimseg] + CICurve[_cu_ctrl][_cu_step]; #if(CurrPos < MaxPos) #for(i, 0, dimension_size(CImatrix, 1) - 1) #local CImult[i] = <0, 0, 0>; #for(j, 0, dimension_size(CImatrix[i], 1) - 1) #local jSegment = j + CurrPos; #local CImult[i] = CImult[i] + CImatrix[i][j] * CICurve[_cu_point][jSegment]; #end #end #end #local CICurve[_cu_curmult] = CImult; #end #local F = CImult[0]; #for(i, 1, dimension_size(CImult, 1) - 1) #local F = (F * Tseg) + CImult[i]; #end F #end #macro CITget(CITCurve, T) /*: Gets the point on the curve at T, after being modulated parameter: CITCurve (CICurve, CICurve): an array of two CICurves. See CIt() macro T (float, [0,1]): position at which to interpolate the curve. return: (vector) returns a point on the curve at modulated T. */ #local nT = CIget(CITCurve[_cu_tcurve], T).x; CIget(CITCurve[_cu_ccurve], nT) #end #if(input_file_name="curve.inc") //tests //#global_settings{ assumed_gamma 1.0 } //#default{ finish{ ambient 0.1 diffuse 0.9 emission 0}} #declare Lin = CIlinear( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Simp3 = CIsimple3( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Bspl2 = CIbspline2( array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Bez2 = CIbezier2( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Gold = CIgolden( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Bez3 = CIbezier3( array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Four3 = CI4p3( array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Herm = CIhermite( array[4]{<0, 0, 0>, <0, 2, 0>, <1, 2, 0>, <1, 0, 0>} ); #declare Bspl3 = CIbspline3( array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>} ); #declare Catmull = CIcatmull( array[6]{<-1, 0, 0>,<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>,<2,0,0>} ); #declare ESC = CIsimple3ESC( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}, 3 //End point Slope control parameter ); #declare MSC = CIsimple3MSC( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}, 0.5 //Mid point Slope control parameter ); #declare IESC = CIsimple3IESC( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}, -.5, 2 //Independent End point Slope control parameters ); #declare Bez2EVC = CIbezier2EVC( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}, -3 //End point Velocity Control parameter ); #declare Lag2 = CIlagrange2( array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}, 0.3 ); #declare Beta = CIbeta( array[6]{<0,0,0>, <0,1,0>, <1,1,0>, <1,0,0>, <2,0,0>, <2,1,0>} 0, 1 //Tension, Bias ) #declare TCB = CItcb( array[6]{<0,0,0>, <0,1,0>, <1,1,0>, <1,0,0>, <2,0,0>, <2,1,0>} 0, 1, 1 //Tension, Continuity, Bias ) #declare Bez3EVC = CIbezier3EVC( array[7]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>, <1,-1,0>, <2,-1,0>, <2,0,0>}, -2 //End point Velocity Control parameter ); #declare Palmer = CIpalmer( array[5]{<0, 0, 0>, <0, 1, 0>, <1, 0.5, 0>, <2,1,0>,<2,0,0>}, 0, 0.5 //Tension, Bias ); // Set curve for test render // #declare Test = Bez3; //CIprint(Test) //get the points on the curve #for(T, 0, 1, 0.01) sphere{CIget(Test, T), 0.02 pigment{rgb 1}} #end //control points sphere{Test[1][0] 0.03 pigment{rgb <1,0,0>}} #for(i, 1, dimension_size(Test[1],1) - 1) sphere{Test[1][i] 0.03 pigment{rgb <1,0,0>}} cylinder{Test[1][i-1], Test[1][i], 0.01 pigment{rgb <1,0,0>}} #end camera { perspective angle 95 location <0.5, 0.5,-2> right x*image_width/image_height look_at <0.5, 0.5, 0.0> } light_source{<0, 0,-3000> color rgb 1} #end #version CURVE_INC_TEMP; #end