#version 3.8; global_settings {assumed_gamma 1.0} camera { location <0.5, 0.5, -1> right x*image_width/image_height up y look_at <0.5, 0.3, 0.5> rotate -x*80 } light_source {<0.5, 0.5, -5> rgb 0.75 shadowless} sky_sphere {pigment {rgb 1}} //#include "BezierInclude.inc" #include "math.inc" #declare E = 0.0000001; #declare Lump = -z*1; #declare del = 0; #declare Patch = 2; // ******* CONTROL POINTS FOR BEZIER PATCH ******************* #switch (Patch) #case (0) #declare P1 = <0.00, 0.00, 0.00>; // uv vector <0, 0> #declare P5 = < 1/3, 0.00, 0.00>; #declare P9 = < 2/3, 0.00, 0.00>; #declare P13 = <1.00, 0.00, 0.00>; // uv vector <1, 0> #declare P2 = <0.00, 1/3+del, 0.00>; #declare P6 = < 1/3+del, 1/3+del, 0.00> - Lump; #declare P10 = < 2/3-del, 1/3+del, 0.00> + Lump; #declare P14 = <1.00, 1/3+del, 0.00>; #declare P3 = <0.00, 2/3-del, 0.00>; #declare P7 = < 1/3+del, 2/3-del, 0.00> + Lump; #declare P11 = < 2/3-del, 2/3-del, 0.00> - Lump; #declare P15 = <1.00, 2/3-del, 0.00>; #declare P4 = <0.00, 1.00, 0.00>; // uv vector <0, 1> #declare P8 = < 1/3, 1.00, 0.00>; #declare P12 = < 2/3, 1.00, 0.00>; #declare P16 = <1.00, 1.00, 0.00>; // uv vector <1, 1> #break #case (1) #declare P1 = <0.00, 0.00, 0.00>; // uv vector <0, 0> #declare P5 = < 1/3, 0.00, 0.00>; #declare P9 = < 2/3, 0.00, 0.00>; #declare P13 = <1.00, 0.00, 0.00>; // uv vector <1, 0> #declare P2 = <0.00, 1/3, 0.00>; #declare P6 = < 1/3, 1/3, 0.00>; #declare P10 = < 2/3, 1/3, 0.00>; #declare P14 = <1.00, 1/3, 0.00>; #declare P3 = <0.00, 2/3, 0.00>; #declare P7 = < 1/3, 2/3, 0.00>; #declare P11 = < 2/3, 2/3, 0.00>; #declare P15 = <1.00, 2/3, 0.00>; #declare P4 = <0.00, 1.00, 0.00> + Lump; // uv vector <0, 1> #declare P8 = < 1/3, 1.00, 0.00>; #declare P12 = < 2/3, 1.00, 0.00>; #declare P16 = <1.00, 1.00, 0.00>; // uv vector <1, 1> #break #case (2) #declare P1 = <0.00, 0.00, 0.00>; // P[0][0] uv vector <0, 0> #declare P5 = <0.00, 1/3+del, 0.00>; // P[0][1] #declare P9 = <0.00, 2/3-del, 0.00>; // P[0][2] #declare P13 = <0.00, 1.00, 0.00>; // P[0][3] uv vector <0, 1> #declare P2 = < 1/3, 0.00, 0.00>; // P[1][0] #declare P6 = < 1/3+del,1/3+del, 0.00> - Lump*0; // P[1][1] #declare P10 = < 1/3+del,2/3-del, 0.00> + Lump; // P[1][2] #declare P14 = < 1/3, 1.00, 0.00>; // P[1][3] #declare P3 = < 2/3, 0.00, 0.00>; // P[2][0] #declare P7 = < 2/3-del,1/3+del, 0.00> - Lump; // P[2][1] #declare P11 = < 2/3-del,2/3-del, 0.00> - Lump*0; // P[2][2] #declare P15 = < 2/3, 1.00, 0.00>; // P[2][3] #declare P4 = <1.00, 0.00, 0.00>; // P[3][0] uv vector <1, 0> #declare P8 = <1.00, 1/3+del, 0.00>; // P[3][1] #declare P12 = <1.00, 2/3-del, 0.00>; // P[3][2] #declare P16 = <1.00, 1.00, 0.00>; // P[3][3] // uv vector <1, 1> #break #else #end #declare ControlPoints = array [4][4] { { P1, P5, P9, P13}, { P2, P6, P10, P14}, { P3, P7 , P11, P15}, { P4, P8 , P12, P16} } #declare PP = array [4][4] { //[0][0-3] { "P[0,0]", "P[0,1]", "P[0,2]", "P[0,3]" } //[1][0-3] { "P[1,0]", "P[1,1]", "P[1,2]", "P[1,3]" } //[2][0-3] { "P[2,0]", "P[2,1]", "P[2,2]", "P[2,3]" } //[3][0-3] { "P[3,0]", "P[3,1]", "P[3,2]", "P[3,3]" } } #declare FastNCM = function (n, k) {prod(i, k+1, n, i) / prod(i, 1, n-k, i)} #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)} // I would like general macros that support derivatives 0-2 and NxM patches #macro MakeVFn_BezierPatch (_Array, _Element) // Assembles formula for one element of Bezier patch // (need to run 3 times for x, y, and z) #local Degree = dimension_size (_Array, 1)-1; function (_U, _V) { #for (j, 0, Degree) #for (i, 0, Degree) #local _P = select (_Element, _Array [i][j].x, _Array [i][j].y, _Array [i][j].z); //#debug concat(Scalar(i, 0, 0), ", ", Scalar(j, 0, 1), "\n") Bernstein(Degree, i, _V) * Bernstein(Degree, j, _U) * _P + #end #end 0} // terminates equation for sum of all Bernstein polynomials in the patch #end // ****** Bernstein polynomials, 1st derivatives, 2nd partial derivatives, and mixed partial derivatives // Bernstein polynomials (16) #declare F = array [4][4] { // U[0]V[0-3] // P[U][V] { function (_U, _V) { 1 * pow (1-_U, 3) * 1 * pow (1-_V, 3)}, // P[0][0] function (_U, _V) { 1 * pow (1-_U, 3) * 3 * _V * pow (1-_V, 2)}, // P[0][1] function (_U, _V) { 1 * pow (1-_U, 3) * 3 * pow (_V, 2) * (1-_V) }, // P[0][2] function (_U, _V) { 1 * pow (1-_U, 3) * pow (_V, 3) } // P[0][3] } // U[1]V[0-3] { function (_U, _V) {3 * _U * pow (1-_U, 2) * 1 * pow (1-_V, 3)}, // P[1][0] function (_U, _V) {3 * _U * pow (1-_U, 2) * 3 * _V * pow (1-_V, 2)}, // P[1][1] function (_U, _V) {3 * _U * pow (1-_U, 2) * 3 * pow (_V, 2) * (1-_V) }, // P[1][2] function (_U, _V) {3 * _U * pow (1-_U, 2) * pow (_V, 3) } // P[1][3] } // U[2]V[0-3] { function (_U, _V) {3 * pow (_U, 2) * (1-_U) * 1 * pow (1-_V, 3)}, // P[2][0] function (_U, _V) {3 * pow (_U, 2) * (1-_U) * 3 * _V * pow (1-_V, 2)}, // P[2][1] function (_U, _V) {3 * pow (_U, 2) * (1-_U) * 3 * pow (_V, 2) * (1-_V) }, // P[2][2] function (_U, _V) {3 * pow (_U, 2) * (1-_U) * pow (_V, 3) } // P[2][3] } // U[3]V[0-3] { function (_U, _V) { pow (_U, 3) * 1 * pow (1-_V, 3)}, // P[3][0] function (_U, _V) { pow (_U, 3) * 3 * _V * pow (1-_V, 2)}, // P[3][1] function (_U, _V) { pow (_U, 3) * 3 * pow (_V, 2) * (1-_V) }, // P[3][2] function (_U, _V) { pow (_U, 3) * pow (_V, 3) } // P[3][3] } } // First partial derivatives with respect to U #declare Fu = array [4][4] { // U[0]V[0-3] // P[U][V] { function (_U, _V) { -3 * pow (1-_U, 2) * 1 * pow (1-_V, 3)}, // P[0][0] function (_U, _V) { -3 * pow (1-_U, 2) * 3 * _V * pow (1-_V, 2)}, // P[0][1] function (_U, _V) { -3 * pow (1-_U, 2) * 3 * pow (_V, 2) * (1-_V) }, // P[0][2] function (_U, _V) { -3 * pow (1-_U, 2) * pow (_V, 3) } // P[0][3] } // U[1]V[0-3] { function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * 1 * pow (1-_V, 3)}, // P[1][0] function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * 3 * _V * pow (1-_V, 2)}, // P[1][1] function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * 3 * pow (_V, 2) * (1-_V) }, // P[1][2] function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * pow (_V, 3) } // P[1][3] } // U[2]V[0-3] { function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * 1 * pow (1-_V, 3)}, // P[2][0] function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * 3 * _V * pow (1-_V, 2)}, // P[2][1] function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * 3 * pow (_V, 2) * (1-_V) }, // P[2][2] function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * pow (_V, 3) } // P[2][3] } // U[3]V[0-3] { function (_U, _V) { 3 * pow (_U, 2) * 1 * pow (1-_V, 3)}, // P[3][0] function (_U, _V) { 3 * pow (_U, 2) * 3 * _V * pow (1-_V, 2)}, // P[3][1] function (_U, _V) { 3 * pow (_U, 2) * 3 * pow (_V, 2) * (1-_V) }, // P[3][2] function (_U, _V) { 3 * pow (_U, 2) * pow (_V, 3) } // P[3][3] } } // First partial derivatives with respect to V #declare Fv = array [4][4] { // U[0]V[0-3] // P[U][V] { function (_U, _V) { 1 * pow (1-_U, 3) * (-3 * pow (1-_V, 2))}, // P[0][0] function (_U, _V) { 1 * pow (1-_U, 3) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[0][1] function (_U, _V) { 1 * pow (1-_U, 3) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[0][2] function (_U, _V) { 1 * pow (1-_U, 3) * 3 * pow (_V, 2) } // P[0][3] } // U[1]V[0-3] { function (_U, _V) { 3 * _U * pow (1-_U, 2) * (-3 * pow (1-_V, 2))}, // P[1][0] function (_U, _V) { 3 * _U * pow (1-_U, 2) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[1][1] function (_U, _V) { 3 * _U * pow (1-_U, 2) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[1][2] function (_U, _V) { 3 * _U * pow (1-_U, 2) * 3 * pow (_V, 2) } // P[1][3] } // U[2]V[0-3] { function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * (-3 * pow (1-_V, 2))}, // P[2][0] function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[2][1] function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[2][2] function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * 3 * pow (_V, 2) } // P[2][3] } // U[2]V[0-3] { function (_U, _V) { pow (_U, 3) * (-3 * pow (1-_V, 2))}, // P[3][0] function (_U, _V) { pow (_U, 3) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[3][1] function (_U, _V) { pow (_U, 3) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[3][2] function (_U, _V) { pow (_U, 3) * 3 * pow (_V, 2) } // P[3][3] } } // Second partial derivatives with respect to U #declare Fuu = array [4][4] { // U[0]V[0-3] // P[U][V] { function (_U, _V) {( -6 * _U + 6) * 1 * pow (1-_V, 3)}, // P[0][0] function (_U, _V) {( -6 * _U + 6) * 3 * _V * pow (1-_V, 2)}, // P[0][1] function (_U, _V) {( -6 * _U + 6) * 3 * pow (_V, 2) * (1-_V) }, // P[0][2] function (_U, _V) {( -6 * _U + 6) * pow (_V, 3) } // P[0][3] } // U[1]V[0-3] { function (_U, _V) {( 18 * _U - 12) * 1 * pow (1-_V, 3)}, // P[1][0] function (_U, _V) {( 18 * _U - 12) * 3 * _V * pow (1-_V, 2)}, // P[1][1] function (_U, _V) {( 18 * _U - 12) * 3 * pow (_V, 2) * (1-_V) }, // P[1][2] function (_U, _V) {( 18 * _U - 12) * pow (_V, 3) } // P[1][3] } // U[2]V[0-3] { function (_U, _V) {( -18 * _U + 6) * 1 * pow (1-_V, 3)}, // P[2][0] function (_U, _V) {( -18 * _U + 6) * 3 * _V * pow (1-_V, 2)}, // P[2][1] function (_U, _V) {( -18 * _U + 6) * 3 * pow (_V, 2) * (1-_V) }, // P[2][2] function (_U, _V) {( -18 * _U + 6) * pow (_V, 3) } // P[2][3] } // U[3]V[0-3] { function (_U, _V) { 6 * _U * 1 * pow (1-_V, 3)}, // P[3][0] function (_U, _V) { 6 * _U * 3 * _V * pow (1-_V, 2)}, // P[3][1] function (_U, _V) { 6 * _U * 3 * pow (_V, 2) * (1-_V) }, // P[3][2] function (_U, _V) { 6 * _U * pow (_V, 3) } // P[3][3] } } // Second partial derivatives with respect to V #declare Fvv = array [4][4] { // U[0]V[0-3] // P[U][V] { function (_U, _V) { 1 * pow (1-_U, 3) * ( -6 * _V + 6) }, // P[0][0] function (_U, _V) { 1 * pow (1-_U, 3) * ( 18 * _V - 12)}, // P[0][1] function (_U, _V) { 1 * pow (1-_U, 3) * ( -18 * _V + 6) }, // P[0][2] function (_U, _V) { 1 * pow (1-_U, 3) * 6 * _V } // P[0][3] } // U[1]V[0-3] { function (_U, _V) { 3 * _U * pow (1-_U, 2) * ( -6 * _V + 6) }, // P[1][0] function (_U, _V) { 3 * _U * pow (1-_U, 2) * ( 18 * _V - 12)}, // P[1][1] function (_U, _V) { 3 * _U * pow (1-_U, 2) * ( -18 * _V + 6) }, // P[1][2] function (_U, _V) { 3 * _U * pow (1-_U, 2) * 6 * _V } // P[1][3] } // U[2]V[0-3] { function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * ( -6 * _V + 6) }, // P[2][0] function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * ( 18 * _V - 12)}, // P[2][1] function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * ( -18 * _V + 6) }, // P[2][2] function (_U, _V) { 3 * pow (_U, 2) * (1-_U) * 6 * _V } // P[2][3] } // U[3]V[0-3] { function (_U, _V) { pow (_U, 3) * ( -6 * _V + 6) }, // P[3][0] function (_U, _V) { pow (_U, 3) * ( 18 * _V - 12)}, // P[3][1] function (_U, _V) { pow (_U, 3) * ( -18 * _V + 6) }, // P[3][2] function (_U, _V) { pow (_U, 3) * 6 * _V } // P[3][3] } } // Mixed partial derivatives with respect to U & V #declare Fuv = array [4][4] { // U[0]V[0-3] // P[U][V] { function (_U, _V) { -3 * pow (1-_U, 2) * (-3 * pow (1-_V, 2))}, // P[0][0] function (_U, _V) { -3 * pow (1-_U, 2) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[0][1] function (_U, _V) { -3 * pow (1-_U, 2) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[0][2] function (_U, _V) { -3 * pow (1-_U, 2) * 3 * pow (_V, 2) } // P[0][3] } // U[1]V[0-3] { function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * (-3 * pow (1-_V, 2))}, // P[1][0] function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[1][1] function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[1][2] function (_U, _V) {(-12 * _U + 9 * pow (_U, 2) + 3) * 3 * pow (_V, 2) } // P[1][3] } // U[2]V[0-3] { function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * (-3 * pow (1-_V, 2))}, // P[2][0] function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[2][1] function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[2][2] function (_U, _V) {( 6 * _U - 9 * pow (_U, 2)) * 3 * pow (_V, 2) } // P[2][3] } // U[3]V[0-3] { function (_U, _V) { 3 * pow (_U, 2) * (-3 * pow (1-_V, 2))}, // P[3][0] function (_U, _V) { 3 * pow (_U, 2) * (-12 * _V + 9 * pow (_V, 2) + 3) }, // P[3][1] function (_U, _V) { 3 * pow (_U, 2) * ( 6 * _V - 9 * pow (_V, 2)) }, // P[3][2] function (_U, _V) { 3 * pow (_U, 2) * 3 * pow (_V, 2) } // P[3][3] } } #macro MakeFn_BezierPatch_du (_Array, _Element) // Assembles formula for one element of 1st derivative of Bezier spline // (need to run 3 times for x, y, and z) [-1 = x, 0 = y, 1 = z] #local Degree = dimension_size (_Array, 1)-1; function (_U, _V) { #for (j, 0, Degree) #for (i, 0, Degree) #local _P = select (_Element, _Array [i][j].x, _Array [i][j].y, _Array [i][j].z); // select vector component Fu[i][j](_V,_U) * _P + #end #end 0} // terminates equation #end #macro MakeFn_BezierPatch_dv (_Array, _Element) // Assembles formula for one element of 1st derivative of Bezier spline // (need to run 3 times for x, y, and z) [-1 = x, 0 = y, 1 = z] #local Degree = dimension_size (_Array, 1)-1; function (_U, _V) { #for (j, 0, Degree) #for (i, 0, Degree) #local _P = select (_Element, _Array [i][j].x, _Array [i][j].y, _Array [i][j].z); // select vector component Fv[i][j](_V,_U) * _P + #end #end 0} // terminates equation #end #macro MakeFn_BezierPatch_d2u (_Array, _Element) // Assembles formula for one element of 1st derivative of Bezier spline // (need to run 3 times for x, y, and z) [-1 = x, 0 = y, 1 = z] #local Degree = dimension_size (_Array, 1)-1; function (_U, _V) { #for (j, 0, Degree) #for (i, 0, Degree) #local _P = select (_Element, _Array [i][j].x, _Array [i][j].y, _Array [i][j].z); // select vector component Fuu[i][j](_U,_V) * _P + #end #end 0} // terminates equation #end #macro MakeFn_BezierPatch_d2v (_Array, _Element) // Assembles formula for one element of 1st derivative of Bezier spline // (need to run 3 times for x, y, and z) [-1 = x, 0 = y, 1 = z] #local Degree = dimension_size (_Array, 1)-1; function (_U, _V) { #for (j, 0, Degree) #for (i, 0, Degree) #local _P = select (_Element, _Array [i][j].x, _Array [i][j].y, _Array [i][j].z); // select vector component Fvv[i][j](_U,_V) * _P + #end #end 0} // terminates equation #end #macro MakeFn_BezierPatch_d2uv (_Array, _Element) // Assembles formula for one element of 1st derivative of Bezier spline // (need to run 3 times for x, y, and z) [-1 = x, 0 = y, 1 = z] #local Degree = dimension_size (_Array, 1)-1; function (_U, _V) { #for (j, 0, Degree) #for (i, 0, Degree) #local _P = select (_Element, _Array [i][j].x, _Array [i][j].y, _Array [i][j].z); // select vector component Fuv[i][j](_U,_V) * _P + #end #end 0} // terminates equation #end // Bezier surface functions #declare VFn_BezierX = MakeVFn_BezierPatch (ControlPoints, -1) #declare VFn_BezierY = MakeVFn_BezierPatch (ControlPoints, 0) #declare VFn_BezierZ = MakeVFn_BezierPatch (ControlPoints, 1) // First partial derivatives #declare Fn_BezierPatch_duX = MakeFn_BezierPatch_du (ControlPoints, -1) #declare Fn_BezierPatch_duY = MakeFn_BezierPatch_du (ControlPoints, 0) #declare Fn_BezierPatch_duZ = MakeFn_BezierPatch_du (ControlPoints, 1) #declare Fn_BezierPatch_dvX = MakeFn_BezierPatch_dv (ControlPoints, -1) #declare Fn_BezierPatch_dvY = MakeFn_BezierPatch_dv (ControlPoints, 0) #declare Fn_BezierPatch_dvZ = MakeFn_BezierPatch_dv (ControlPoints, 1) // Second partial derivatives #declare Fn_BezierPatch_d2uX = MakeFn_BezierPatch_d2u (ControlPoints, -1) #declare Fn_BezierPatch_d2uY = MakeFn_BezierPatch_d2u (ControlPoints, 0) #declare Fn_BezierPatch_d2uZ = MakeFn_BezierPatch_d2u (ControlPoints, 1) #declare Fn_BezierPatch_d2vX = MakeFn_BezierPatch_d2v (ControlPoints, -1) #declare Fn_BezierPatch_d2vY = MakeFn_BezierPatch_d2v (ControlPoints, 0) #declare Fn_BezierPatch_d2vZ = MakeFn_BezierPatch_d2v (ControlPoints, 1) #declare Fn_BezierPatch_d2uvX = MakeFn_BezierPatch_d2uv (ControlPoints, -1) #declare Fn_BezierPatch_d2uvY = MakeFn_BezierPatch_d2uv (ControlPoints, 0) #declare Fn_BezierPatch_d2uvZ = MakeFn_BezierPatch_d2uv (ControlPoints, 1) // Outer product (OP) // inner product (IP) {dot product} // trace {sum of elements of diagonal of OP} /* Wedge product = OP (u, v) - OP (v, u) [0 u1v2-v1u2 u1v3-v1u3 u2v1-v2u1 0 u2v3-v2u3 u3v1-v3u1 u3v2-v3u2 0 ] = [yz zy xy ax ay az bz by bz] */ // SFn = Scalar function // VFn = Vector function // DOT product #declare SFn_vdot = function (ax, ay, az, bx, by, bz) {ax*bx + ay*by + az*bz} // CROSS product #declare VFn_vcrossX = function (ax, ay, az, bx, by, bz) {ay*bz - az*by} #declare VFn_vcrossY = function (ax, ay, az, bx, by, bz) {az*bx - ax*bz} #declare VFn_vcrossZ = function (ax, ay, az, bx, by, bz) {ax*by - ay*bx} // Length or Norm #declare SFn_vlength = function (ax, ay, az) {sqrt(ax*ax + ay*ay + az*az)} // Surface normal = du X dv #declare VFn_SurfaceNormalX = function (_U, _V) { VFn_vcrossX ( // cross product of Fn_BezierPatch_duX (_U, _V), Fn_BezierPatch_duY (_U, _V), Fn_BezierPatch_duZ (_U, _V), // du Fn_BezierPatch_dvX (_U, _V), Fn_BezierPatch_dvY (_U, _V), Fn_BezierPatch_dvZ (_U, _V) // dv ) } #declare VFn_SurfaceNormalY = function (_U, _V) { VFn_vcrossY ( // cross product of Fn_BezierPatch_duX (_U, _V), Fn_BezierPatch_duY (_U, _V), Fn_BezierPatch_duZ (_U, _V), // du Fn_BezierPatch_dvX (_U, _V), Fn_BezierPatch_dvY (_U, _V), Fn_BezierPatch_dvZ (_U, _V) // dv ) } #declare VFn_SurfaceNormalZ = function (_U, _V) { VFn_vcrossZ ( // cross product of Fn_BezierPatch_duX (_U, _V), Fn_BezierPatch_duY (_U, _V), Fn_BezierPatch_duZ (_U, _V), // du Fn_BezierPatch_dvX (_U, _V), Fn_BezierPatch_dvY (_U, _V), Fn_BezierPatch_dvZ (_U, _V) // dv ) } /* // experimenting with a gradient-style method //----------------------------------- #declare d = 0.0001; #declare VFn_SurfaceNormalX = function (_U, _V) { VFn_vcrossX ( (VFn_BezierX (_U+d, _V) - VFn_BezierX (_U-d, _V))/(2*d), (VFn_BezierY (_U+d, _V) - VFn_BezierY (_U-d, _V))/(2*d), (VFn_BezierZ (_U+d, _V) - VFn_BezierZ (_U-d, _V))/(2*d), (VFn_BezierX (_U, _V+d) - VFn_BezierX (_U, _V-d))/(2*d), (VFn_BezierY (_U, _V+d) - VFn_BezierY (_U, _V-d))/(2*d), (VFn_BezierZ (_U, _V+d) - VFn_BezierZ (_U, _V-d))/(2*d) ) } #declare VFn_SurfaceNormalY = function (_U, _V) { VFn_vcrossY ( (VFn_BezierX (_U+d, _V) - VFn_BezierX (_U-d, _V))/(2*d), (VFn_BezierY (_U+d, _V) - VFn_BezierY (_U-d, _V))/(2*d), (VFn_BezierZ (_U+d, _V) - VFn_BezierZ (_U-d, _V))/(2*d), (VFn_BezierX (_U, _V+d) - VFn_BezierX (_U, _V-d))/(2*d), (VFn_BezierY (_U, _V+d) - VFn_BezierY (_U, _V-d))/(2*d), (VFn_BezierZ (_U, _V+d) - VFn_BezierZ (_U, _V-d))/(2*d) ) } #declare VFn_SurfaceNormalZ = function (_U, _V) { VFn_vcrossZ ( (VFn_BezierX (_U+d, _V) - VFn_BezierX (_U-d, _V))/(2*d), (VFn_BezierY (_U+d, _V) - VFn_BezierY (_U-d, _V))/(2*d), (VFn_BezierZ (_U+d, _V) - VFn_BezierZ (_U-d, _V))/(2*d), (VFn_BezierX (_U, _V+d) - VFn_BezierX (_U, _V-d))/(2*d), (VFn_BezierY (_U, _V+d) - VFn_BezierY (_U, _V-d))/(2*d), (VFn_BezierZ (_U, _V+d) - VFn_BezierZ (_U, _V-d))/(2*d) ) } */ // Calculate surface unit normal vector #declare VFn_UnitNormalX = function (_U, _V) {VFn_SurfaceNormalX (_U, _V) / // normal SFn_vlength (VFn_SurfaceNormalX (_U, _V), VFn_SurfaceNormalY (_U, _V), VFn_SurfaceNormalZ (_U, _V))} // length #declare VFn_UnitNormalY = function (_U, _V) {VFn_SurfaceNormalY (_U, _V) / // normal SFn_vlength (VFn_SurfaceNormalX (_U, _V), VFn_SurfaceNormalY (_U, _V), VFn_SurfaceNormalZ (_U, _V))} // length #declare VFn_UnitNormalZ = function (_U, _V) {VFn_SurfaceNormalZ (_U, _V) / // normal SFn_vlength (VFn_SurfaceNormalX (_U, _V), VFn_SurfaceNormalY (_U, _V), VFn_SurfaceNormalZ (_U, _V))} // length // coefficients of the First Fundamental Form // Fu dot Fu #declare SFn_Iuu = function (_U, _V) { SFn_vdot ( // dot product of Fn_BezierPatch_duX (_U, _V), Fn_BezierPatch_duY (_U, _V), Fn_BezierPatch_duZ (_U, _V), // du Fn_BezierPatch_duX (_U, _V), Fn_BezierPatch_duY (_U, _V), Fn_BezierPatch_duZ (_U, _V) // du ) } // Fv dot Fv #declare SFn_Ivv = function (_U, _V) { SFn_vdot ( // dot product of Fn_BezierPatch_dvX (_U, _V), Fn_BezierPatch_dvY (_U, _V), Fn_BezierPatch_dvZ (_U, _V), // dv Fn_BezierPatch_dvX (_U, _V), Fn_BezierPatch_dvY (_U, _V), Fn_BezierPatch_dvZ (_U, _V) // dv ) } // Fu dot Fv #declare SFn_Iuv = function (_U, _V) { SFn_vdot ( // dot product of Fn_BezierPatch_duX (_U, _V), Fn_BezierPatch_duY (_U, _V), Fn_BezierPatch_duZ (_U, _V), // du Fn_BezierPatch_dvX (_U, _V), Fn_BezierPatch_dvY (_U, _V), Fn_BezierPatch_dvZ (_U, _V) // dv ) } // coefficients of the Second Fundamental Form // Fuu dot n #declare SFn_IIuu = function (_U, _V) { SFn_vdot ( // dot product of Fn_BezierPatch_d2uX (_U, _V), Fn_BezierPatch_d2uY (_U, _V), Fn_BezierPatch_d2uZ (_U, _V), // d2u VFn_UnitNormalX (_U, _V), VFn_UnitNormalY (_U, _V), VFn_UnitNormalZ (_U, _V) // unit normal ) } // Fvv dot n #declare SFn_IIvv = function (_U, _V) { SFn_vdot ( // dot product of Fn_BezierPatch_d2vX (_U, _V), Fn_BezierPatch_d2vY (_U, _V), Fn_BezierPatch_d2vZ (_U, _V), // d2v VFn_UnitNormalX (_U, _V), VFn_UnitNormalY (_U, _V), VFn_UnitNormalZ (_U, _V) // unit normal ) } // Fuv dot n #declare SFn_IIuv = function (_U, _V) { SFn_vdot ( // dot product of Fn_BezierPatch_d2uvX (_U, _V), Fn_BezierPatch_d2uvY (_U, _V), Fn_BezierPatch_d2uvZ (_U, _V), // duv VFn_UnitNormalX (_U, _V), VFn_UnitNormalY (_U, _V), VFn_UnitNormalZ (_U, _V) // unit normal ) } #declare GaussianCurvature = function (_U, _V) { ( SFn_IIuu(_U, _V)*SFn_IIvv(_U, _V) - pow(SFn_IIuv(_U, _V),2) )/ ( SFn_Iuu (_U, _V)*SFn_Ivv (_U, _V) - pow(SFn_Iuv (_U, _V),2) ) } #declare MeanCurvature = function (_U, _V) { ( SFn_Iuu(_U, _V)*SFn_IIvv(_U, _V) -2*SFn_Iuv(_U, _V)*SFn_IIuv(_U, _V) + SFn_Ivv(_U, _V)*SFn_IIuu(_U, _V))/ ( SFn_Iuu(_U, _V)*SFn_Ivv (_U, _V) - pow(SFn_Iuv(_U, _V),2) ) } // Scene union { #declare Normals = yes; #declare UV = yes; #if (Normals) #declare Step = 0.05; #declare MinC = 0; #declare MaxC = 0; #for (V, 0, 1, Step) #for (U, 0, 1, Step) #local P = ; #local N = ; #local du = /50; #local dv = /50; sphere {P, 0.0025 pigment {rgb 1}} //cylinder {0, N/25 0.0025 pigment {rgb <0, 1, 0>} translate P} //#local K = GaussianCurvature (U, V); #local K = MeanCurvature (U, V); //cylinder {0, -(N/25) 0.0025 pigment {rgb } translate P} cylinder {0, du 0.00125 pigment {rgb x} translate P} cylinder {0, dv 0.00125 pigment {rgb y} translate P} cylinder {0, -(N/50) 0.00125 pigment {rgb 0.5} translate P} #declare MinC = min (K, MinC); #declare MaxC = max (K, MaxC); #end #end #debug concat ("Min = ", str (MinC, 0, 16), "\n") #debug concat ("Max = ", str (MaxC, 0, 16), "\n") #end // PLOT CONTROL POINTS #for (V, 0, 3) #for (U, 0, 3) #local P = ControlPoints[U][V]; sphere {P, 0.01 pigment {rgb z}} text { ttf "timrom.ttf", PP[U][V], 0.04, 0.0 scale 0.03 pigment {rgb 0} translate P translate -z*0.02} #end #debug "\n" #end #if (UV) #declare Range = MaxC-MinC; #declare T_Curvature = texture { //pigment {function { ((GaussianCurvature (x, y) - MinC)/Range)*0.99} pigment {function { ((MeanCurvature (x, y) - MinC)/Range)*0.99} color_map { [ 0.0 rgb z ] [ 0.55 rgb y ] //[ 0.66 rgb x+y ] [ 1.0 rgb x ] } // end color map } finish {diffuse 1} } #end bicubic_patch { type 1 flatness 0 u_steps 5 v_steps 5 uv_vectors <0,0>,<1,0>,<1,1>,<0,1> P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12, P13, P14, P15, P16 #if (UV) uv_mapping texture {T_Curvature} #else texture {pigment {rgbt <0.6, 0.6, 0.4, 0.0>} finish {specular 0.4}} #end } // end bicubic_patch 1 union { //text { ttf "timrom.ttf", "Gaussian curvature mapping of a Bezier surface", 0.04, 0.0 scale 0.03 pigment {rgb 0} translate 0.8*y} text { ttf "timrom.ttf", "Mean curvature mapping of a Bezier surface", 0.04, 0.0 scale 0.03 pigment {rgb 0} translate 0.8*y} text { ttf "timrom.ttf", "Bill Walker (\"Bald Eagle\") 2022", 0.04, 0.0 scale 0.03 pigment {rgb 0} translate -0.1*y} rotate -x*80 } //rotate x*90 } //box {<-0.25, -0.25, 0.01>, <1.25, 1.25, 0.02> texture {pigment {checker rgb 0 rgb 1 scale 0.125} finish {diffuse 0.7 specular 0.2}}}