#version 3.8; global_settings {assumed_gamma 1.0 } #include "colors.inc" #include "functions.inc" #include "math.inc" #declare SS = 50; #declare DataScale = 5; #declare NP = 20; #declare Seed = seed (123); camera { location //location <0, 50, 0> right x*image_width/image_height up y look_at //look_at <0, 0, 0.01> } light_source {<0, 10, -50> rgb 1 shadowless} // tiling pattern in the xz plane plane {z, 0.25 pigment {rgb 1}} /* plane {y, 0 pigment{ tiling 15 // 1~24 Pattern_Number color_map{ [ 0.0 color rgb<1,1,1>*1] [ 0.5 color rgb<1,0.5,0>*1] [ 1.0 color rgb<1,1,1>*0] } // end color_map scale 0.5 // rotate<-90,0,0> // align to xy plane } // end pigment } */ // prod (i, b, n, a) for i from b to n, multiply a-sub-i's // sum (i, b, n, a) for i from b to n, add a-sub-i's //------------------------------------------------------------------------------ // Fast NChooseM #declare FastNCM = function (n, k) {prod(i, k+1, n, i) / prod(i, 1, n-k, i)} //------------------------------------------------------------------------------ // Cycles through array elements backwards or forwards // -1 of an array with dimension size N gives element N-1 // N of an array with dimension size N gives element 0 //#declare Cycle = function (_N, _Size) {select (_N, (_Size-1)+mod(_N+1, _Size), mod(_N, _Size) )} //------------------------------------------------------------------------------ #macro Choose (_N, _Array) // returns an element from an array // uses mod () to cycle through array no matter the size of _N or _Array // *** Add optional dimension to access other than 1 #local _Size = dimension_size (_Array,1); #local __N = select (_N, (_Size + _N), 0, _N); #local _Item = mod (__N, _Size); #local _Result = _Array [_Item]; _Result #end //------------------------------------------------------------------------------ #declare Round = function (Value) {ceil(Value - 0.5)} //------------------------------------------------------------------------------ #declare PowerFn = function(s, p) { select(p, pow(s, p), 1, pow(s, p)) } //------------------------------------------------------------------------------ #declare Bernstein = function (Deg, K, S) {FastNCM(Deg, K) * PowerFn(S, K) * PowerFn (1 - S, Deg - K)} //------------------------------------------------------------------------------ #macro Point (VG, UG, Uc, Vc) // called during interpolation loop due to using Current [V][U] array //#debug concat(Scalar(VG, 0, 0), ", ", Scalar(UG, 0, 1), "\n") #local Degree = 3; #local P = #for (j, 0, 3) #for (i, 0, 3) //#debug concat(Scalar(VG, 0, 0), ", ", Scalar(UG, 0, 0), ", ", Scalar(i, 0, 0), ", ", Scalar(j, 0, 1), "\n") Bernstein(Degree, i, Vc) * Bernstein(Degree, j, Uc) * StitchedArray [VG][UG][i][j] + #end #end 0; P #end // generate a function that allows the equivalent of returning the value of a 1-dimensional array // ( Lagrange Polynomial Interpolation) /* functions cannot be used to access array elements, however, for any array index, a function can be written to return a scalar value of an array. This can also be used to emulate a switch-case block. Array elements cannot be accessed from within a function, so a function cannot be written to directly sum the products of the terms derived from the array elements, however a MACRO _can_ be used to assemble a Lagrange Polynomial Interpolation function from the array's elements. Evaluating that function then serves to return the abscissa that is equivalent to the array element for any ordinate index value. */ // But can we use a SPLINE function to build the equation _directly? Probably yes. // Advantages: no control points to worry about. // With a spline function, you don't need to set the array size. #declare Points = array; #for (NNN, 0, NP-1) #local Points[NNN] = ; #end #declare N = dimension_size (Points, 1)-1; #macro LPI (Array) #local _N = dimension_size (Array, 1)-1; #local Equation = function (X) { //#debug "\n Equation = function (X) {" #for (i, 0, _N) #local Xi = Array [i].x; #local Yi = Array [i].y; Yi //#debug str (Yi, 0, 0) #for (j, 0, _N) #local Xj = Array [j].x; //#debug concat(Scalar(VG, 0, 0), ", ", Scalar(UG, 0, 0), ", ", Scalar(i, 0, 0), ", ", Scalar(j, 0, 1), "\n") #if (j = i) // don't add this term to equation #else * (X - Xj)/(Xi - Xj) //#debug concat ("*(X-",str (Xj, 0, 0), ")/(", str (Xi, 0, 0), "-", str (Xj, 0, 0), ")") #end // end if j = i #end // end for j + //#debug "+ " #end // end for i 0} //#debug "0} \n\n" Equation #end #declare PointSpline = function { spline { linear_spline #for (NNN, 0, NP-1) NNN, Points[NNN], #end } } #declare SSS = function (_Q) {select ((_Q!=3)-1, -1, 0)} #for (QQ, 0, 5) #local Ans = SSS(QQ); #debug concat (str (QQ, 0, 2), " returns ", str (Ans, 0, 2), "\n") #end #debug "\n" //define interpolation function directly #declare LPInterpolate = function (NN, X) { sum (i, 0, NN, PointSpline(i).y * prod (j, 0, NN, select ((j!=i)-1, 1, (X - PointSpline(j).x) / (PointSpline(i).x - PointSpline(j).x)) ) ) } #for (XX, 0, N) #debug concat (str (XX, 0, 2), " = ",str (PointSpline(XX).x, 0, 2), ", ", str (PointSpline (XX).y, 0, 2), ", ", str (PointSpline (XX).z, 0, 2) "\n") #end #debug "\n" // invoke macro to define interpolation equation by that method #declare Fitted = LPI (Points); #declare YScale = 1; union { #for (XX, 0, N) sphere {Points [XX] 0.5 pigment {rgb 0}} #end scale <1, YScale, 1> } #debug "Macro generated function data: \n" union { #for (XX, 0, N, 0.05) #local Current = ; sphere {Current 0.125 pigment {rgb <1, 0, 0>}} #if (XX >0) cylinder {Old, Current 0.125 pigment {rgb <1, 0, 0>}} #end #declare Old = Current; //#debug concat (str (XX, 0, 2), ", ", str (Fitted(XX), 0, 2), "\n") #end scale <1, YScale, 1> } #debug "Pure spline function data: \n" union { #for (XX, 0, N, 0.05) #local Current = ; sphere { 0.125 pigment {rgb <0, 0, 1>} translate y*0.25} #if (XX >0) cylinder {Old, Current 0.125 pigment {rgb <0, 0, 1>} translate y*0.25} #end #declare Old = Current; //#debug concat (str (XX, 0, 2), ", ", str (Fitted(XX), 0, 2), "\n") #end scale <1, YScale, 1> }