// ------------------------------------------------ // Cube-based mesh generation and smoothing // // Created by Bill Pragnell, 2020 // This file is licensed under the terms of the CC-LGPL // // The CubeMeshShape macro traces points onto an object to make a mesh2 version of it, and // includes options for function-based surface perturbation and simple nearest-neighbour // smoothing. // Instead of the equirectangular approach taken by the meshrelief include, this version // uses a point distribution based on a cube surface. This avoids the 'pinching' at the // poles that the previous version suffered from. // // ------------------------------------------------ // Usage: // // CubeMeshShape(N, Obj, ObjSmooth, NormalSmooth, IncOpt, SaveName) // // N points per smallest edge of point cuboid (in proportion with target bounding box) // Obj target object // ObjSmooth number of passes for target object smoothing filter (0 = no smoothing) // NormalSmooth calculate normals for mesh (yes/no) // IncOpt include option (cmsNO_INC, cmsSAVE_INC, cmsLOAD_INC) // SaveName filename (no extension) for .inc file (can be empty if not used) // // N must be >= 4 because of the way the edges are joined together. // N gives 12*(N-1)^2 triangles. // Target object considerations: // - should be non-self-occluding from all angles (the mesh will still be created // correctly, but overhangs will be filled in). // - should be carefully bounded; the bounding box is used to locate and scale the // test points, so for best results should be as minimal as possible. // // Additionally, there are two optional perturbation stages, each with multiple options. // For ease of use, use #declare to globally declare the perturbation parameters before calling // the macro (substitute 'N' for '1' or '2' in the identifiers below). // // CMSPerturbFnN Function to be used (must return in range 0-1) // CMSPerturbDN Max distance of perturbation (can be negative) // CMSPerturbTexMapN Texture map to be applied to the function // CMSPerturbSmoothN Smoothing passes (0 = no smoothing) // // Perturbation functions should be float functions (they can be based on pigment functions by // using the . accessors). If a Fn is defined, the corresponding D must also be defined, because // they are used together (a separate D is necessary for ease of use by the texturing code). // // If perturbation is not required, simply omit the parameter declarations and the macro will // produce only the mesh of the target object. // // When loading a saved .inc file, the macro assumes that any #declared pertubation texture maps // are still available. When switching between 'save' and 'load', be sure not to change any of // them. // ------------------------------------------------ // Minimal usage example // Creates a mesh box with 2 smoothing passes on the base object, a single smoothed displacement // stage defined by a simple granite pattern, textured with a white-red map, smooth mesh normals // and no include file output. /* #declare TargetObject = box { -1, 1 } #declare T1 = texture { pigment { rgb 1 } finish { ambient 0 } } #declare T2 = texture { pigment { rgb <0.75, 0.25, 0> } finish { ambient 0 } } #declare PigmentFn = function { pigment { granite color_map { [0 rgb 0][1 rgb 1] } } } #declare CMSPerturbFn1 = function { PigmentFn(x,y,z).x } #declare CMSPerturbD1 = 0.1; #declare CMSPerturbSmooth1 = 2; #declare CMSPerturbTexMap1 = texture_map { [0 T1] [1 T2] }; CubeMeshShape(50, TargetObject, 2, yes, cmsNO_INC, "") */ // ------------------------------------------------ // ------------------------------------------------ // Helper macros // ------------------------------------------------ // Generate cube surface points. // // Cube points define the indexing; their order needs to be well-defined. // Top face, rastered starting from <-1,1,-1>. // Centre band of 4 faces, rastered in a wrapped loop starting from <-1,y,-1>. // Bottom face, rastered starting from <-1,-1,-1>. // Note that the bottom face has not been reversed, so triangulation of these // nodes must take this into account for normals. // #macro GenerateCubePositions(N, S, PointsOut) #local FaceN = N*N; #local StripN = (N-2)*(N-1)*4; #local TotalN = FaceN + FaceN + StripN; #local Points = array[TotalN]; #local PC = 0; #local D = 2/(N-1); // Top square #local Y = 1; #local Z = -1; #for (I, 0, N-1) #local X = -1; #for (J, 0, N-1) #local Points[PC] = S * ; #local PC = PC + 1; #local X = X + D; #end #local Z = Z + D; #end // Middle band #for (I, 0, N-3) #local Y = Y - D; #local X = -1; #local Z = -1; #local SubVal = array[8] { 1, 0, 0, 1, -1, 0, 0, -1 } #for (J, 0, 3) #local XM = SubVal[2*J]; #local ZM = SubVal[2*J+1]; #for (K, 0, N-2) #local Points[PC] = S * ; #local PC = PC + 1; #local X = X + XM*D; #local Z = Z + ZM*D; #end #end #end // Bottom square #local Y = -1; #local Z = -1; #for (I, 0, N-1) #local X = -1; #for (J, 0, N-1) #local Points[PC] = S * ; #local PC = PC + 1; #local X = X + D; #end #local Z = Z + D; #end // Return point list #declare PointsOut = Points; #end // A helper for the neighbour generation; make a list of the indices around // the edge of the top/bottom faces, to match the centre band ordering. // #macro CubeFaceLoopIndices(N, UberN, CornerI, IndicesOut) #local FaceN = N*N; #local RingN = 4*(N-1); #local Indices = array[RingN]; #local PC = 0; // Bottom row #for (I, 0, N-2) #local Indices[PC] = CornerI + I; #local PC = PC + 1; #end // Right row #for (I, 0, N-2) #local Indices[PC] = CornerI + (N-1) + I*UberN; #local PC = PC + 1; #end // Top row, in reverse #for (I, 0, N-2) #local Indices[PC] = CornerI + (N-1) + (N-1)*UberN - I; #local PC = PC + 1; #end // Left row #for (I, 0, N-2) #local Indices[PC] = CornerI + (N-1) + (N-1)*UberN - (N-1) - I*UberN; #local PC = PC + 1; #end #declare IndicesOut = Indices; #end // Generate cube point neighbour indices (L, R, D, U) // // These are stored in an array, a 4-neighbour entry per cube point, and define the // indices of the neighbouring point. This saves us the trouble of tricky index gymnastics // when triangulating or smoothing. // // Top face inner, L->R = +x, D->U = +z // Top face edge, L->R = ac around edge, D = centre band top row, U = top face inner // Centre band, L->R = ac around band, D->U = +y // #macro GenerateCubeNeighbours(N, NeighboursOut) #local FaceN = N * N; #local StripRowN = (N-1) * 4; #local StripN = (N-2) * StripRowN; #local TotalN = FaceN + FaceN + StripN; #local Neighbours = array[TotalN][4]; // ----- Top/bottom faces excluding edges // (well, including some edges, but these will be overwritten later) #local FirstTopCentre = N + 1; #local LastTopCentre = FaceN - N - 2; #for (TI, FirstTopCentre, LastTopCentre) // Top face #local Neighbours[TI][0] = TI-1; #local Neighbours[TI][1] = TI+1; #local Neighbours[TI][2] = TI-N; #local Neighbours[TI][3] = TI+N; // Bottom face #local BI = TI + FaceN + StripN; #local Neighbours[BI][0] = BI-1; #local Neighbours[BI][1] = BI+1; #local Neighbours[BI][2] = BI+N; // these are other way around so #local Neighbours[BI][3] = BI-N; // that normal is still outwards #end // ----- Top and bottom face edges, plus centre band top and bottom rows #local TopEdgePoints = array[1]; #local BotEdgePoints = array[1]; CubeFaceLoopIndices(N, N, 0, TopEdgePoints) CubeFaceLoopIndices(N, N, FaceN+StripN, BotEdgePoints) #local TopInnerEdgePoints = array[1]; #local BotInnerEdgePoints = array[1]; CubeFaceLoopIndices(N-2, N, N+1, TopInnerEdgePoints) CubeFaceLoopIndices(N-2, N, FaceN+StripN+N+1, BotInnerEdgePoints) // Edge loops have same counts top and bottom #local EdgeN = dimension_size(TopEdgePoints, 1); #local InnerEdgeN = dimension_size(TopInnerEdgePoints, 1); #local TopStripStart = FaceN; #local BotStripStart = FaceN + StripN - StripRowN; // Iterate over edge loop #local NextCornerN = 0; #local JJ = 0; #for (J, 0, EdgeN-1) #local TI = TopEdgePoints[J]; #local BI = BotEdgePoints[J]; // left #if (J = 0) // leftmost wrap #local Neighbours[TI][0] = TopEdgePoints[EdgeN-1]; // top edge #local Neighbours[BI][0] = BotEdgePoints[EdgeN-1]; // bottom edge #local Neighbours[TopStripStart+J][0] = TopStripStart + StripRowN - 1; // top strip first row #local Neighbours[BotStripStart+J][0] = BotStripStart + StripRowN - 1; // bottom strip last row #else // regular left #local Neighbours[TI][0] = TopEdgePoints[J-1]; // top edge #local Neighbours[BI][0] = BotEdgePoints[J-1]; // bottom edge #local Neighbours[TopStripStart+J][0] = TopStripStart + J - 1; // top strip first row #local Neighbours[BotStripStart+J][0] = BotStripStart + J - 1; // bottom strip last row #end // right #if (J = EdgeN-1) // rightmost wrap #local Neighbours[TI][1] = TopEdgePoints[0]; #local Neighbours[BI][1] = BotEdgePoints[0]; #local Neighbours[TopStripStart+J][1] = TopStripStart; #local Neighbours[BotStripStart+J][1] = BotStripStart; #else // regular right #local Neighbours[TI][1] = TopEdgePoints[J+1]; #local Neighbours[BI][1] = BotEdgePoints[J+1]; #local Neighbours[TopStripStart+J][1] = TopStripStart + J + 1; #local Neighbours[BotStripStart+J][1] = BotStripStart + J + 1; #end // towards centre band #local Neighbours[TI][2] = TopStripStart + J; // top edge down #local Neighbours[BI][3] = BotStripStart + J; // bottom edge up #local Neighbours[TopStripStart+J][2] = TopStripStart + StripRowN + J; // top strip first row down #local Neighbours[BotStripStart+J][3] = BotStripStart - StripRowN + J; // bottom strip last row up // towards faces // uses the inner loops for the face edge #if (J = NextCornerN) // special case for corners; don't increment inner loop index #local Neighbours[TI][3] = TopInnerEdgePoints[JJ]; // top edge #local Neighbours[BI][2] = BotInnerEdgePoints[JJ]; // bottom edge #local NextCornerN = NextCornerN + (N-1); #else // regular edge point inwards #local Neighbours[TI][3] = TopInnerEdgePoints[JJ]; // top edge #local Neighbours[BI][2] = BotInnerEdgePoints[JJ]; // bottom edge // only increment inner loop index if we're not yet at the next corner #if (J+1 < NextCornerN) #local JJ = JJ + 1; #end #if (JJ = InnerEdgeN) #local JJ = 0; #end #end // inwards for strip rows is face edges #local Neighbours[TopStripStart+J][3] = TopEdgePoints[J]; // top strip first row up #local Neighbours[BotStripStart+J][2] = BotEdgePoints[J]; // bottom strip first row down #end // ----- Centre band #local RowStart = FaceN + StripRowN; #local RowEnd = RowStart + StripRowN; #local NRows = N - 4; #for (J, 0, NRows-1) #for (I, RowStart, RowEnd-1) // left #if (I = RowStart) #local Neighbours[I][0] = RowEnd - 1; #else #local Neighbours[I][0] = I - 1; #end // right #if (I = RowEnd-1) #local Neighbours[I][1] = RowStart; #else #local Neighbours[I][1] = I + 1; #end // down #local Neighbours[I][2] = I + StripRowN; // up #local Neighbours[I][3] = I - StripRowN; #end #local RowStart = RowStart + StripRowN; #local RowEnd = RowEnd + StripRowN; #end // ----- Finished! Copy array to output #declare NeighboursOut = Neighbours; #end // Generate surface points given a target object and position list. // #macro GenerateSurfacePoints(Obj, CubePositions, PointsOut) #local Min = min_extent(Obj); #local Max = max_extent(Obj); #local R = vlength(Max-Min); #local P = (Min+Max)/2; #local Num = dimension_size(CubePositions, 1); #local Points = array[Num]; #for (I, 0, Num-1) #local Points[I] = trace(Obj, P + R*vnormalize(CubePositions[I]), -CubePositions[I]); #end #declare PointsOut = Points; #end // Perform Laplacian smoothing on the surface point list and neighbour point list. // #macro SmoothSurfacePoints(Points, Neighbours) #local Num = dimension_size(Points, 1); #local NewPoints = array[Num]; #for (I, 0, Num-1) #local AP = Points[Neighbours[I][0]]; #for (J, 1, 3) #local AP = AP + Points[Neighbours[I][J]]; #end #local NewPoints[I] = AP / 4; #end #declare Points = NewPoints; #end // Perturb the surface point list, also creating UV coords and corresponding // texture index if necessary. // #macro PerturbSurfacePoints(Points, Neighbours, PertFn, PertD, UV, TexIndex, TexI) #local Num = dimension_size(Points, 1); #local NewPoints = array[Num]; #for (I, 0, Num-1) #local PertVal = PertFn(Points[I].x,Points[I].y,Points[I].z); #local NewPoints[I] = Points[I] + Normal(I, Points, Neighbours) * PertD * PertVal; #if (TexI != -1) #ifdef (UV[I]) #local Def = true; #else #local Def = false; #end #if (!Def | PertVal > 0) #declare UV[I] = PertVal; #declare TexIndex[I] = TexI; #end #end #end #declare Points = NewPoints; #end // Make a single quad as two triangles in mesh2 format, with optional texture index. // Also write the triangles to the include file if necessary. // #macro Quad(I1, I2, I3, I4, P, T, Tex, Save) #local N = vcross(P[I4]-P[I1], P[I3]-P[I2]); #local V = (P[I1]+(P[I4]-P[I1])/2) - (P[I2]+(P[I3]-P[I2])/2); #if (vdot(N,V) < 0) , #if (Tex) , T[I1], T[I4], T[I3] #end , #if (Tex) , T[I1], T[I2], T[I4] #end #if (Save) #write (IncFile, ",\n<",Ix1,",",Ix4,",",Ix3,">") #if (Tex) #write (IncFile, ",",T[I1],",",T[I4],",",T[I3]) #end #write (IncFile, ",\n<",Ix1,",",Ix2,",",Ix4,">") #if (Tex) #write (IncFile, ",",T[I1],",",T[I2],",",T[I4]) #end #end #else , #if (Tex) , T[I1], T[I2], T[I3] #end , #if (Tex) , T[I2], T[I4], T[I3] #end #if (Save) #write (IncFile, ",\n<",Ix1,",",Ix2,",",Ix3,">") #if (Tex) #write (IncFile, ",",T[I1],",",T[I2],",",T[I3]) #end #write (IncFile, ",\n<",Ix2,",",Ix4,",",Ix3,">") #if (Tex) #write (IncFile, ",",T[I2],",",T[I4],",",T[I3]) #end #end #end #end // Calculate a point normal based on nearest neighbours. // #macro Normal(I, Points, Neighbours) (-vnormalize(vcross( Points[Neighbours[I][1]] - Points[Neighbours[I][0]], Points[Neighbours[I][3]] - Points[Neighbours[I][2]] ))) #end // Generate a list of textures from the input texture map list. // Also write the textures to the include file if necessary. // #macro GenerateTextures(TexMaps, TexListOut, Save) #local C = 0; #local PartTexList = array[2]; #for (I, 0, 1) #ifdef (TexMaps[I]) #local PartTexList[C] = texture { uv_mapping gradient x texture_map { TexMaps[I] } } #if (Save) #write (IncFile, "#local Texture",C+1," = texture {\n") #write (IncFile, " uv_mapping gradient x texture_map { CMSPerturbTexMap",I+1," }\n") #write (IncFile, "}\n") #end #local C = C + 1; #end #end #local TexList = array[C]; #for (I, 0, C-1) #local TexList[I] = PartTexList[I]; #end #declare TexListOut = TexList; #end // Generate a mesh2 given the point list, normals, UV coords and textures. // Also write the mesh2 definition to the include file if necessary. // #macro GenerateMesh(N, Points, Neighbours, UVs, TextureIs, Textures, Smooth, Save) #local Num = dimension_size(Points, 1); #local FaceN = N*N; #local BandN = (N-2) * (N-1) * 4; #local BandRowN = (N-1) * 4; #local FOff = FaceN + BandN; #local DoTextures = no; #ifdef (UVs[0]) // if any UVs have been created, need UV and texture lists #local DoTextures = yes; #end mesh2 { vertex_vectors { Num #if (Save) #write (IncFile, "mesh2 {\nvertex_vectors {\n",Num) #end #for (I, 0, Num-1) , Points[I] #if (Save) #write (IncFile, ",\n<",Points[I].x,",",Points[I].y,",",Points[I].z,">") #end #end } #if (Save) #write (IncFile, "\n}\n") #end #if (Smooth) normal_vectors { Num #if (Save) #write (IncFile, "normal_vectors {\n",Num) #end #for (I, 0, Num-1) , Normal(I, Points, Neighbours) #if (Save) #write (IncFile, ",\n",Normal(I, Points, Neighbours)) #end #end } #if (Save) #write (IncFile, "\n}\n") #end #end #if (DoTextures) uv_vectors { Num #if (Save) #write (IncFile, "uv_vectors {\n",Num) #end #for (I, 0, Num-1) , #if (Save) #write (IncFile, ",\n<",UVs[I],",0>") #end #end } #if (Save) #write (IncFile, "\n}\n") #end #local NumTex = dimension_size(Textures,1); texture_list { NumTex #if (Save) #write (IncFile, "texture_list {\n",NumTex) #end #for (I, 0, NumTex-1) , texture { Textures[I] } #if (Save) #write (IncFile, ",\ntexture { Texture",I+1," }") #end #end } #if (Save) #write (IncFile, "\n}\n") #end #end #local NumTris = 2*6*(N-1)*(N-1); face_indices { NumTris #if (Save) #write (IncFile, "face_indices {\n",NumTris) #end // top / bottom faces #for (I, 0, N-2) #for (J, 0, N-2) #local Ix1 = I*N + J; #local Ix2 = Ix1 + 1; #local Ix3 = (I+1)*N + J; #local Ix4 = Ix3 + 1; // top face tris Quad(Ix1, Ix2, Ix3, Ix4, Points, TextureIs, DoTextures, Save) // bottom face tris Quad(Ix1+FOff, Ix2+FOff, Ix3+FOff, Ix4+FOff, Points, TextureIs, DoTextures, Save) #end #end // centre band #for (I, 0, N-3) #for (J, 0, BandRowN-1) #local Ix1 = FaceN + I*BandRowN + J; #if (J = BandRowN-1) #local Ix2 = FaceN + I*BandRowN; #else #local Ix2 = FaceN + I*BandRowN + J + 1; #end #local Ix3 = Neighbours[Ix1][3]; #local Ix4 = Neighbours[Ix2][3]; Quad(Ix1, Ix2, Ix3, Ix4, Points, TextureIs, DoTextures, Save) // bottom row goes down too #if (I = N-3) #local Ix5 = Neighbours[Ix1][2]; #local Ix6 = Neighbours[Ix2][2]; Quad(Ix5, Ix6, Ix1, Ix2, Points, TextureIs, DoTextures, Save) #end #end #end } inside_vector y } #if (Save) #write (IncFile, "\n}\ninside_vector y\n}\n") #end #debug concat("Made cubemesh with ",str(NumTris,5,0)," triangles.\n") #end #declare cmsNO_INC = 0; #declare cmsSAVE_INC = 1; #declare cmsLOAD_INC = 2; // ------------------------------------------------ // Main user-facing macro // ------------------------------------------------ // Make a cube-based mesh version of a target shape. // // N points per smallest edge of point cuboid (in proportion with target bounding box) // Obj target object // ObjSmooth number of passes for target object smoothing filter (0 = no smoothing) // NormalSmooth calculate normals for mesh (yes/no) // IncOpt include option (cmsNO_INC, cmsSAVE_INC, cmsLOAD_INC) // SaveName filename (no extension) for .inc file (can be empty if not used) // // Additionally, there are two optional perturbation stages, each with multiple options. // For ease of use, use #declare to globally declare the perturbation parameters before calling // the macro (substitute N in the identifier for either 1 or 2). // If the perturbation is not required, simply omit the parameter declaration and the macro will // produce only the mesh of the target object. // // CMSPerturbFnN Function to be used (must return in range 0-1) // CMSPerturbDN Max distance of perturbation (can be negative) // CMSPerturbTexMapN Texture map to be applied to the function // CMSPerturbSmoothN Smoothing passes (0 = no smoothing) // // Perturbation functions should be float functions (they can be based on pigment functions by using // the .x or .gray accessors). // #macro CubeMeshShape(N, Obj, ObjSmooth, NormalSmooth, IncOpt, SaveName) #if (IncOpt = cmsLOAD_INC) // loading include - all other args irrelevant #include concat(SaveName,".inc") #else // not loading - build the mesh // if saving, open incfile here #if (IncOpt = cmsSAVE_INC) // file handle is global, so no need to pass it around #fopen IncFile concat(SaveName,".inc") write #local Save = yes; #else #local IncName = ""; #local Save = no; #end // set up arrays #local PositionList = array[1]; #local NeighbourList = array[1]; #local PointList = array[1]; #local TexList = array[1]; #local CurTex = 0; // calculate dimensions of cuboid #local ObjMin = min_extent(Obj); #local ObjMax = max_extent(Obj); #local OX = ObjMax.x - ObjMin.x; #local OY = ObjMax.y - ObjMin.y; #local OZ = ObjMax.z - ObjMin.z; #local OMax = max(OX, OY, OZ); #local Scale = ; // generate cube vectors and neighbour indices GenerateCubePositions(N, Scale, PositionList) GenerateCubeNeighbours(N, NeighbourList) // generate base object surface points GenerateSurfacePoints(Obj, PositionList, PointList) // smooth base object #for (I, 0, ObjSmooth-1) SmoothSurfacePoints(PointList, NeighbourList) #end // make UV and texture index arrays based on point list size #local NPoints = dimension_size(PointList, 1); #local UVList = array[NPoints]; #local TexIndexList = array[NPoints]; // first stage perturbation #ifdef (CMSPerturbFn1) // sort out texture index #ifdef (CMSPerturbTexMap1) #local TexI = CurTex; #local CurTex = CurTex + 1; #else #local TexI = -1; #end // apply perturbation to surface points PerturbSurfacePoints(PointList, NeighbourList, CMSPerturbFn1, CMSPerturbD1, UVList, TexIndexList, TexI) #end // smoothing of first perturbed points #ifdef (CMSPerturbSmooth1) #for (I, 0, CMSPerturbSmooth1-1) SmoothSurfacePoints(PointList, NeighbourList) #end #end // second stage perturbation #ifdef (CMSPerturbFn2) // sort out texture index #ifdef (CMSPerturbTexMap2) #local TexI = CurTex; #else #local TexI = -1; #end // apply perturbation to surface points PerturbSurfacePoints(PointList, NeighbourList, CMSPerturbFn2, CMSPerturbD2, UVList, TexIndexList, TexI) #end // smoothing of second perturbed points #ifdef (CMSPerturbSmooth2) #for (I, 0, CMSPerturbSmooth2-1) SmoothSurfacePoints(PointList, NeighbourList) #end #end // make textures for UV mapping #local NTexMaps = 0; #local TexMaps = array[2]; #ifdef (CMSPerturbTexMap1) #local NTexMaps = NTexMaps + 1; #local TexMaps[0] = CMSPerturbTexMap1; #end #ifdef (CMSPerturbTexMap2) #local NTexMaps = NTexMaps + 1; #local TexMaps[1] = CMSPerturbTexMap2; #end #if (NTexMaps > 0) GenerateTextures(TexMaps, TexList, Save) #end // make mesh GenerateMesh(N, PointList, NeighbourList, UVList, TexIndexList, TexList, NormalSmooth, Save) #end #end