// 3D extension of Bernstein polynomials for deforming meshes // by Bill Walker "Bald Eagle" // Based on an idea / request of Josh English // Takes 64 control points for "4 Bezier bicubic patches", // shows the points and control grid, // uses Bernstein polynomials to create an array of points to use as triangle corners, // calculates the vertex normals, // and creates 4 meshes {} of smooth_triangles #version 3.8; global_settings {assumed_gamma 1.0} camera { location <0, 2, -12> right x up y look_at <0, 0.5, 1> } //light_source {<0, 5, -1> rgb 1.0} light_source {<1, 0, -15> rgb 0.5} sky_sphere {pigment {rgb <0.6, 0.6, 1.0>}} #declare Ext = 2.5; #declare Square = union { triangle {<-Ext, 0, -Ext>, <-Ext, 0, Ext>, } triangle {, <-Ext, 0, Ext>, } texture {checker} } #declare Shade = 0.2; #declare _Red = texture {pigment {rgb <1, 0, 0>*Shade} finish {diffuse 0.7}} #declare _Green = texture {pigment {rgb <0, 1, 0>*Shade} finish {diffuse 0.7}} #declare _Blue = texture {pigment {rgb <0, 0, 1>*Shade} finish {diffuse 0.7}} #declare Axes = union { #declare Line = 0.005; #declare Base = Line*2; #declare Length = 5; #declare Ext = 0.25; cylinder {<0, 0, 0>, , Line texture {_Red}} cone {, Base, , 0 texture {_Red} } cylinder {<0, 0, 0>, <0, Length, 0>, Line texture {_Green}} cone {<0, Length, 0>, Base, <0, Length+Ext, 0>, 0 texture {_Green}} union { cylinder {<0, 0, 0>, <0, 0, Length>, Line texture {_Blue}} cone {<0, 0, Length>, Base, <0, 0, Length+Ext>, 0 texture {_Blue}} translate -z*(Length+Ext) } //no_shadow } object {Axes} //========================================== #declare Sr = 0.0075; // sphere radius #declare T_S = texture {pigment {rgb <0.5, 0.0, 0.5>} finish {specular 0.4}} // sphere texture #declare Cr = Sr/3; #declare T_C = texture {pigment {rgb <0.0, 0.75, 0.75>} finish {specular 0.4}} // cylinder texture #declare T_F = texture {pigment {rgb <1.0, 0.3, 0.0>} finish {diffuse 0.7 specular 0.3 emission 0.1} normal {agate bump_size 0.1 scale 0.04}} // face texture #declare T_FI = texture {pigment {rgb <0.0, 0.1, 0.0>} finish {diffuse 0.6 specular 0.4} normal {agate bump_size 0.1 scale 0.1}} // face interior texture #declare Control_Points = yes; #declare CPcolor = <1.0, 0.0, 0.0>; #declare ControlGrid = yes; #declare CGcolor = <1.0, 0.25, 0.25>; #declare Vertices = yes; #declare Edges = yes; #declare Faces = yes; #declare Pattern = yes; #declare UStep = 3; #declare VStep = 3; #declare WStep = 3; #declare Lump = -z*1; // deviation from planar patch #declare delta = 0; // shift of center control points from equidistant positions //========================================== //############################################################## // Control points for tricubic Bezier parallelpiped in xyz space // Patches are aligned in xz plane, and stacked in the y direction // Bottom layer // Front row #declare P1 = <0.00, 0.00, 0.00>; // uvw vector <0, 0, 0> #declare P2 = < 1/3, 0.00, 0.00>; #declare P3 = < 2/3, 0.00, 0.00>; #declare P4 = <1.00, 0.00, 0.00>; // uvw vector <1, 0, 0> // Second row #declare P5 = <0.00, 0.00, 1/3>; #declare P6 = < 1/3, 0.00, 1/3>; #declare P7 = < 2/3, 0.00, 1/3>; #declare P8 = <1.00, 0.00, 1/3>; // Third row #declare P9 = <0.00, 0.00, 2/3>; #declare P10 = < 1/3, 0.00, 2/3>; #declare P11 = < 2/3, 0.00, 2/3>; #declare P12 = <1.00, 0.00, 2/3>; // Rear row #declare P13 = <0.00, 0.00, 1.00>; // uvw vector <0, 0, 1> #declare P14 = < 1/3, 0.00, 1.00>; #declare P15 = < 2/3, 0.00, 1.00>; #declare P16 = <1.00, 0.00, 1.00>; // uvw vector <1, 0, 1> // Second Layer // Front row #declare P17 = <0.00, 1/3, 0.00>; // uvw vector <0, 1/3, 0> #declare P18 = < 1/3, 1/3, 0.00>; #declare P19 = < 2/3, 1/3, 0.00>; #declare P20 = <1.00, 1/3, 0.00>; // uvw vector <1, 1/3, 0> // Second row #declare P21 = <0.00, 1/3, 1/3>; #declare P22 = < 1/3, 1/3, 1/3>; #declare P23 = < 2/3, 1/3, 1/3>; #declare P24 = <1.00, 1/3, 1/3>; // Third row #declare P25 = <0.00, 1/3, 2/3>; #declare P26 = < 1/3, 1/3, 2/3>; #declare P27 = < 2/3, 1/3, 2/3>; #declare P28 = <1.00, 1/3, 2/3>; // Rear row #declare P29 = <0.00, 1/3, 1.00>; // uvw vector <0, 1/3, 1> #declare P30 = < 1/3, 1/3, 1.00>; #declare P31 = < 2/3, 1/3, 1.00>; #declare P32 = <1.00, 1/3, 1.00>; // uvw vector <1, 1/3, 1> // Third Layer // Front row #declare P33 = <0.00, 2/3, 0.00>; // uvw vector <0, 2/3, 0> #declare P34 = < 1/3, 2/3, 0.00>; #declare P35 = < 2/3, 2/3, 0.00>; #declare P36 = <1.00, 2/3, 0.00>; // uvw vector <1, 2/3, 0> // Second row #declare P37 = <0.00, 2/3, 1/3>; #declare P38 = < 1/3, 2/3, 1/3>; #declare P39 = < 2/3, 2/3, 1/3>; #declare P40 = <1.00, 2/3, 1/3>; // Third row #declare P41 = <0.00, 2/3, 2/3>; #declare P42 = < 1/3, 2/3, 2/3>; #declare P43 = < 2/3, 2/3, 2/3>; #declare P44 = <1.00, 2/3, 2/3>; // Rear row #declare P45 = <0.00, 2/3, 1.00>; // uvw vector <0, 2/3, 1> #declare P46 = < 1/3, 2/3, 1.00>; #declare P47 = < 2/3, 2/3, 1.00>; #declare P48 = <1.00, 2/3, 1.00>; // uvw vector <1, 2/3, 1> // Top Layer // Front row #declare P49 = <0.00, 1.00, 0.00>; // uvw vector <0, 1, 0> #declare P50 = < 1/3, 1.00, 0.00>; #declare P51 = < 2/3, 1.00, 0.00>; #declare P52 = <1.00, 1.00, 0.00>; // uvw vector <1, 1, 0> // Second row #declare P53 = <0.00, 1.00, 1/3>; #declare P54 = < 1/3, 1.00, 1/3>; #declare P55 = < 2/3, 1.00, 1/3>; #declare P56 = <1.00, 1.00, 1/3>; // Third row #declare P57 = <0.00, 1.00, 2/3>; #declare P58 = < 1/3, 1.00, 2/3>; #declare P59 = < 2/3, 1.00, 2/3>; #declare P60 = <1.00, 1.00, 2/3>; // Rear row #declare P61 = <0.00, 1.00, 1.00>; // uvw vector <0, 1, 1> #declare P62 = < 1/3, 1.00, 1.00>; #declare P63 = < 2/3, 1.00, 1.00>; #declare P64 = <1.00, 1.00, 1.00>; // uvw vector <1, 1, 1> #macro MakeControlPointArray () // Create a 3D array of the control points #declare ControlPoints = array [4][4][4] { { {P1, P2, P3, P4 }, {P5, P6, P7, P8 }, {P9, P10, P11, P12}, {P13, P14, P15, P16} } { {P17, P18, P19, P20}, {P21, P22, P23, P24}, {P25, P26, P27, P28}, {P29, P30, P31, P32} } { {P33, P34, P35, P36}, {P37, P38, P39, P40}, {P41, P42, P43, P44}, {P45, P46, P47, P48} } { {P49, P50, P51, P52}, {P53, P54, P55, P56}, {P57, P58, P59, P60}, {P61, P62, P63, P64} } } #end // Simple, prototypical Bezier spline expressed as a Bernstein polynomial // #declare Bezier = function (A, B, C, D, T) {A*pow(1-T,3) + 3*B*T*pow(1-T,2) + 3*C*pow(T,2)*(1-T) + D*pow(T,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)} /* Original 2D version of macro to create 2D functions #macro F_Bezier (_Array, _Element) // Creates a function for one basis vector (axis) of a Bezier patch // (need to run 3 times for x, y, and z) #local Degree = 3; 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); Bernstein(Degree, i, _V) * Bernstein(Degree, j, _U) * _P + #end #end 0} // terminates equation for sum of all Bernstein polynomials in the patch #end */ // new 3D version of macro to create functions for 3D Bezier parallelpiped #macro F_Bezier (_Array, _Element) // Creates a function for one basis vector (axis) of a Bezier patch // (need to run 3 times for x, y, and z) #local Degree = 3; function (_U, _V, _W) { #for (k, 0, Degree) #for (j, 0, Degree) #for (i, 0, Degree) #local _P = select (_Element, _Array [i][j][k].x, _Array [i][j][k].y, _Array [i][j][k].z); Bernstein(Degree, i, _W) * Bernstein(Degree, j, _V) * Bernstein(Degree, k, _U) * _P + #end #end #end 0 } // terminates equation for sum of all Bernstein polynomials in the patch #end #macro MakeControlPoints () // gotta undef functions if running more than once #undef BezierX #undef BezierY #undef BezierZ // define the component functions for the parametric representation of the Bezier patch #declare BezierX = F_Bezier (ControlPoints, -1) // x #declare BezierY = F_Bezier (ControlPoints, 0) // y #declare BezierZ = F_Bezier (ControlPoints, 1) // z #end #macro MakeMeshArray () // this allows loops to use integers intead of fractional counters #declare UIStep = 1/UStep; #declare VIStep = 1/VStep; #declare WIStep = 1/WStep; #declare MeshArray = array [UStep+1] [VStep+1] [WStep+1]; #for (UI, 0, UStep) // integer steps for loop #for (VI, 0, VStep) // integer steps for loop #for (WI, 0, WStep) // integer steps for loop #local U = UI * UIStep; // convert parameter to 0-1 domain for use in Bezier functions #local V = VI * VIStep; // convert parameter to 0-1 domain for use in Bezier functions #local W = WI * WIStep; // convert parameter to 0-1 domain for use in Bezier functions #declare MeshArray [UI][VI][WI] = ; #end #end #end #end #macro BezierParallelpiped () //#local BPP = union { #if (Control_Points) #for (X, 0, 3) #for (Y, 0, 3) #for (Z, 0, 3) sphere {ControlPoints[X][Y][Z], 0.01 pigment {rgb CPcolor} no_shadow} #end #end #end #end // end if ControlPoints #if (ControlGrid) #for (Z, 0, 3) #for (Y, 0, 3) #for (X, 0, 3) sphere {ControlPoints[Z][Y][X], Sr pigment {rgb CGcolor} no_shadow} #if (X > 0) cylinder {ControlPoints[Z][Y][X], OldPoint, Cr pigment {rgb CGcolor} no_shadow} #end #declare OldPoint = ControlPoints[Z][Y][X]; #end // end for X #end // end for Y #end // end for Z #for (Z, 0, 3) #for (X, 0, 3) #for (Y, 0, 3) sphere {ControlPoints[Z][Y][X], Sr pigment {rgb CGcolor} no_shadow} #if (Y > 0) cylinder {ControlPoints[Z][Y][X], OldPoint, Cr pigment {rgb CGcolor} no_shadow} #end #declare OldPoint = ControlPoints[Z][Y][X]; #end // end for Y #end // end for X #end // end for Z // added for y direction #for (X, 0, 3) #for (Y, 0, 3) #for (Z, 0, 3) sphere {ControlPoints[Z][Y][X], Sr pigment {rgb CGcolor} no_shadow} #if (Z > 0) cylinder {ControlPoints[Z][Y][X], OldPoint, Cr pigment {rgb CGcolor} no_shadow} #end #declare OldPoint = ControlPoints[Z][Y][X]; #end // end for Y #end // end for X #end // end for Z #end // end if ControlGrid // Wireframe components #if (Vertices | Edges) #declare N1 = dimension_size (MeshArray, 1)-1; #declare N2 = dimension_size (MeshArray, 2)-1; #declare N3 = dimension_size (MeshArray, 3)-1; #for (U, 0, N1) #for (V, 0, N2) #for (W, 0, N3) #if (U != N1 & V != N2) #if (Vertices) sphere {MeshArray [U] [V] [W] Sr texture {T_S}} sphere {MeshArray [U+1][V+1][W] Sr texture {T_S}} sphere {MeshArray [U] [V+1][W] Sr texture {T_S}} //sphere {MeshArray [U][V+1] Sr texture {T_S}} sphere {MeshArray [U+1][V] [W] Sr texture {T_S}} //sphere {MeshArray [U+1][V] Sr texture {T_S}} #end // end if Vertices #if (Edges) cylinder {MeshArray [U] [V] [W] MeshArray [U+1][V+1][W] Cr texture {T_C}} cylinder {MeshArray [U+1][V+1][W] MeshArray [U] [V+1][W] Cr texture {T_C}} cylinder {MeshArray [U] [V+1][W] MeshArray [U] [V] [W] Cr texture {T_C}} cylinder {MeshArray [U] [V] [W] MeshArray [U+1][V] [W] Cr texture {T_C}} cylinder {MeshArray [U+1][V] [W] MeshArray [U+1][V+1][W] Cr texture {T_C}} cylinder {MeshArray [U+1][V+1][W] MeshArray [U] [V] [W] Cr texture {T_C}} #end #end // end if U, V #end // end for W #end // end for V #end // end for U #end // end if Vertices or Edges // Bezier bicubic patch surface #if (Faces) mesh { #declare N1 = dimension_size (MeshArray, 1)-1; #declare N2 = dimension_size (MeshArray, 2)-1; #declare N3 = dimension_size (MeshArray, 3)-1; #for (U, 0, N1) #for (V, 0, N2) #for (W, 0, N3) #if (U != N1 & V != N2) triangle { MeshArray [U][V][W], MeshArray [U+1][V+1][W], MeshArray [U] [V+1][W] } triangle { MeshArray [U][V][W], MeshArray [U+1][V] [W], MeshArray [U+1][V+1][W] } #end // end if U, V #end // end for W #end // end for V #end // end for U texture {T_F} interior_texture {T_FI} } #end // end if Faces //rotate x*60 } // end union #end // Now we need a test mesh: #declare MeshVertices = array [42] { <0.809017,0.309017,0.5>, <0.5,0.809017,0.309017>, <0.309017,0.5,0.809017>, <1,0,0>, <0.809017,-0.309017,0.5>, <-0.5,-0.809017,-0.309017>, <0,-1,0>, <-0.5,-0.809017,0.309017>, <-0.850651,0.525731,0>, <-0.5,0.809017,-0.309017>, <-0.809017,0.309017,-0.5>, <0.525731,0,0.850651>, <0.309017,-0.5,0.809017>, <-1,0,0>, <-0.809017,-0.309017,-0.5>, <0.850651,0.525731,0>, <0,-0.850651,0.525731>, <-0.525731,0,-0.850651>, <-0.309017,-0.5,-0.809017>, <0,0.850651,-0.525731>, <0.309017,0.5,-0.809017>, <0.5,0.809017,-0.309017>, <0.5,-0.809017,0.309017>, <-0.309017,-0.5,0.809017>, <-0.850651,-0.525731,0>, <-0.309017,0.5,-0.809017>, <-0.809017,-0.309017,0.5>, <0,-0.850651,-0.525731>, <-0.5,0.809017,0.309017>, <-0.809017,0.309017,0.5>, <0.5,-0.809017,-0.309017>, <0.309017,-0.5,-0.809017>, <0.850651,-0.525731,0>, <0.809017,-0.309017,-0.5>, <0,0.850651,0.525731>, <-0.309017,0.5,0.809017>, <0.525731,0,-0.850651>, <0.809017,0.309017,-0.5>, <0,1,0>, <-0.525731,0,0.850651>, <0,0,1>, <0,0,-1> } // declare the mesh using a macro, so the code is easily reused many times #macro MakeMesh (Distorted) mesh2 { vertex_vectors { 42, #for (V, 0, 41) #if (Distorted) DistortedVertex (MeshVertices [V]), #else MeshVertices [V], #end #end } face_indices { 80, <0,1,2>, <3,0,4>, <5,6,7>, <8,9,10>, <11,12,4>, <13,10,14>, <15,0,3>, <6,16,7>, <17,18,14>, <20,19,21>, <4,12,22>, <23,16,12>, <5,7,24>, <9,19,25>, <14,18,5>, <24,26,13>, <18,27,5>, <29,28,8>, <30,27,31>, <27,6,5>, <21,1,15>, <7,16,23>, <32,33,3>, <33,30,31>, <8,28,9>, <15,1,0>, <35,34,28>, <4,22,32>, <8,10,13>, <22,16,6>, <10,25,17>, <32,22,30>, <29,35,28>, <33,31,36>, <30,22,6>, <37,20,21>, <12,16,22>, <19,38,21>, <39,35,29>, <11,2,40>, <0,11,4>, <1,34,2>, <13,14,24>, <3,37,15>, <3,4,32>, <10,17,14>, <17,25,41>, <0,2,11>, <41,31,18>, <32,30,33>, <40,2,35>, <2,34,35>, <13,26,29>, <38,34,1>, <9,28,38>, <41,18,17>, <3,33,37>, <13,29,8>, <33,36,37>, <25,19,20>, <14,5,24>, <30,6,27>, <24,7,26>, <26,7,23>, <10,9,25>, <40,12,11>, <21,38,1>, <37,21,15>, <9,38,19>, <40,35,39>, <28,34,38>, <40,23,12>, <36,31,41>, <39,23,40>, <26,39,29>, <36,20,37>, <41,20,36>, <26,23,39>, <31,27,18>, <41,25,20> } inside_vector <0,1,0> } #end // We need to create the basis mesh so we can get extents and calculate u, v, & w parameters of the vertices // Basically, what percent toward the max_extent of the bounding box is each vertex, in all 3 directions MakeControlPointArray () MakeControlPoints () MakeMeshArray () #declare StartMesh = MakeMesh (0); #declare Min = min_extent (StartMesh); #declare Max = max_extent (StartMesh); #declare Range = Max-Min; // now that we have the fundamental mesh information, we can proceed with the rest of the mesh deformation code #macro DistortedVertex (Vert) //#debug concat ("V = ", vstr(3, Vert, ", ", 0, 3), " \n") #local VertexUVW = (Vert-Min)/Range; #local U = VertexUVW.x; //-Min.x; #local V = VertexUVW.y; //-Min.y; #local W = VertexUVW.z; //-Min.z; #local UVW = ; //#debug concat ("UVW = ", vstr(3, UVW, ", ", 0, 3), " \n") // this applies the Bezier "spline" deformation . . . . . and undoes the u,v,w scaling and translation to the origin * Range + Min #end // Simple macro to add random amounts of noise to each direction of the control points // this is a global change to the control point array // to reset to the original, just run MakeControlPointArray () instead to create the array from the basis control points #macro DistortControlPoints (RandSeed1, RandSeed2, RandSeed3, Amt1, Amt2, Amt3) #for (U, 0, 3) #for (V, 0, 3) #for (W, 0, 3) #declare ControlPoints[U][V][W] = ControlPoints[U][V][W] + <(rand(RandSeed1)*2 - 1)*Amt1, (rand(RandSeed2)*2 - 1)*Amt2, (rand(RandSeed3)*2 - 1)*Amt3>; #end // end for W #end // end for V #end // end for U #end //----------------------------------------------------------------------------- // Now we make the scene: // Original mesh2 {} union { object {BezierParallelpiped () translate x*Max.x*2} object {StartMesh texture {pigment {rgb z*0.2} finish {specular 0.1}}} translate -x*Max.x*3 } // Set seed values for random perturbation #declare Seed1 = seed (123); #declare Seed2 = seed (pi*100); #declare Seed3 = seed (1/3*1000); // apply distortion DistortControlPoints (Seed1, Seed2, Seed3, 0, 10, 5) // create functions MakeControlPoints () // create the Bezier parallelpiped MakeMeshArray () // apply Bezier deformation to original mesh2 {} union { object {BezierParallelpiped () translate x*Max.x*3} object {MakeMesh (1) texture {pigment {rgb y} finish {specular 0.1}} } translate x*Max.x }