I want to work on subdivision surfaces. Unfortunately, I don’t have any source code to start with. I need some Mathematica codes for applying Catmull-Clark and Doo-Sabin methods. I would like to request anyone on this platform who have such codes to kindly share them with me.
The codes I have from: http://hakenberg.de/subdivision/enclosed_volume.htm are giving the error code:
Part::partd: Part specification DooSabinListAlt[[1]] is longer than depth of object.
I will appreciate your help.
Indeed, I have some code for Catmull-Clark subdivision and I planned to post it here for quite some time. This seems to be a good opportunity.
The code is optimized for performance, so it involves a lot of CompiledFunction and SparseArray hacks. I am sorry if you find it somewhat unidiomatic.
CatmullClarkSubdivisionMatrix creates the, well, subdivision matrix while the actual subdivision is performed by CatmullClarkSubdivide. The code below assumes that the surface is a manifold-like polyhedral mesh in $\mathbb{R}^3$, possibly with boundary and not necessarily orientable.
Before we start, you might be also interested in Loop subdivision; that one is implemented here.
First, we have the load the code from the section "Implementation" below. Then we can employ the function CatmullClarkSubdivide:
pts = N[{{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, {1, 1, 1}, {-1, 1, 1}}];
polys = {{4, 3, 2, 1}, {1, 2, 6, 5}, {2, 3, 7, 6}, {8, 7, 3, 4}, {5, 8, 4, 1}, {5, 6, 7, 8}};
M = polymesh[pts, polys];
MList = NestList[CatmullClarkSubdivide, M, 4];
GraphicsRow[
GridMeshPlot /@ MList,
ImageSize -> Full
]
A somewhat more complex example:
R = ExampleData[{"Geometry3D", "Triceratops"}, "MeshRegion"];
M = polymesh[MeshCoordinates[R], MeshCells[R, 2, "Multicells" -> True][[1, 1]]];
MList = NestList[CatmullClarkSubdivide, M, 3];
GraphicsRow[
MeshPlot /@ MList,
ImageSize -> Full
]
Two ways to subdivide the boundary:
pts = N[{{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, {1, 1, 1}, {-1, 1, 1}}];
polys = {{4, 3, 2, 1},(*{1,2,6,5},*){2, 3, 7, 6}, {8, 7, 3, 4}, {5, 8, 4, 1}, {5, 6, 7, 8}};
M = polymesh[pts, polys];
With averaging, applying the standard $(1/8, 1/2, 1/8)$-subdivision along the boundary curves:
MList = NestList[
CatmullClarkSubdivide[#, "AverageBoundaryPoints" -> True] &,
M, 4];
GraphicsRow[GridMeshPlot /@ MList, ImageSize -> Full]
And without averaging:
getEdgesFromPolygons = Compile[{{f, _Integer, 1}},
Table[
{
Min[Compile`GetElement[f, i], Compile`GetElement[f, Mod[i + 1, Length[f], 1]]],
Max[Compile`GetElement[f, i], Compile`GetElement[f, Mod[i + 1, Length[f], 1]]]
},
{i, 1, Length[f]}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
getSubdividedPolygons =
Compile[{{qq, _Integer, 1}, {ee, _Integer, 1}, {n, _Integer}},
Table[
{
Compile`GetElement[qq, i],
Compile`GetElement[ee, i],
n,
Compile`GetElement[ee, Mod[i - 1, Length[qq], 1]]
},
{i, 1, Length[qq]}],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
AccumulateIntegerList = Compile[{{list, _Integer, 1}},
Block[{c = 0, r = 0},
Table[
If[i <= Length[list],
r = c; c += Compile`GetElement[list, i]; r,
c
]
, {i, 1, Length[list] + 1}]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
cExtractIntegerFromSparseMatrix = Compile[
{{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer,
1}, {background, _Integer},
{i, _Integer}, {j, _Integer}},
Block[{k, c},
k = Compile`GetElement[rp, i] + 1;
c = Compile`GetElement[rp, i + 1];
While[k < c + 1 && Compile`GetElement[ci, k] != j, ++k];
If[k == c + 1, background, Compile`GetElement[vals, k]]
],
RuntimeAttributes -> {Listable},
Parallelization -> True,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];
ToPack = Developer`ToPackedArray;
polymesh::usage = "";
polymesh /: polymesh[points0_, polygons0_] :=
Module[{polygons},
polygons = ToPack[polygons0];
polymesh[
Association[
"MeshCoordinates" -> ToPack[N[points0]],
"MeshCells" -> Association[
0 -> Partition[Range[Length[points0]], 1],
1 -> DeleteDuplicates[ToPack[Flatten[getEdgesFromPolygons[polygons], 1]]],
2 -> polygons
]
]
]
];
polymesh /: MeshCoordinates[M_polymesh] := M[[1]][["MeshCoordinates"]];
polymesh /: MeshCells[M_polymesh, d_Integer] := M[[1]][["MeshCells", Key[d]]];
polymesh /: MeshCellCount[M_polymesh, d_Integer] := Length[MeshCells[M, d]];
GridMeshPlot::usage = "";
polymesh /: GridMeshPlot[M_polymesh] := Graphics3D[{
ColorData[97][1], Specularity[White, 30], EdgeForm[{Thin, Black}];
GraphicsComplex[MeshCoordinates[M], Polygon[MeshCells[M, 2]]],
},
Lighting -> "Neutral",
Boxed -> False
];
MeshPlot::usage = "";
polymesh /: MeshPlot[M_polymesh] := Graphics3D[{
ColorData[97][1], Specularity[White, 30], EdgeForm[],
GraphicsComplex[MeshCoordinates[M], Polygon[MeshCells[M, 2]]],
},
Lighting -> "Neutral",
Boxed -> False
];
SignedPolygonsNeighEdges::usage = "";
polymesh /: SignedPolygonsNeighEdges[M_polymesh] :=
Module[{edges, n, A00, i, j},
edges = MeshCells[M, 1];
n = MeshCellCount[M, 0];
A00 = SparseArray`SparseArraySort@SparseArray[
Rule[
Join[edges, Transpose[Transpose[edges][[{2, 1}]]]],
Join[Range[1, Length[edges]], Range[-1, -Length[edges], -1]]
],
{n, n}
];
{i, j} = Transpose[Join @@ With[{cf = Compile[{{p, _Integer, 1}},
Transpose[{p, RotateLeft[p]}],
RuntimeAttributes -> {Listable},
Parallelization -> True
]},
cf[MeshCells[M, 2]]
]];
Internal`PartitionRagged[
cExtractIntegerFromSparseMatrix[
A00["NonzeroValues"], A00["RowPointers"],
Flatten[A00["ColumnIndices"]], 0, i, j
],
Length /@ MeshCells[M, 2]
]
];
SubdividedPolygons::usage = "";
polymesh /: SubdividedPolygons[M_polymesh] :=
With[{
n0 = MeshCellCount[M, 0],
n1 = MeshCellCount[M, 1],
n2 = MeshCellCount[M, 2]
},
Flatten[getSubdividedPolygons[
MeshCells[M, 2],
Abs[SignedPolygonsNeighEdges[M]] + n0,
Range[1 + n0 + n1, n0 + n1 + n2]
], 1]
];
getConnectivityMatrix::usage = "";
getConnectivityMatrix[n_Integer, cells_List] :=
With[{m = Length[cells]},
If[m > 0,
Module[{A, lens, nn, rp},
lens = Compile[{{cell, _Integer, 1}},
Length[cell],
CompilationTarget -> "WVM",
RuntimeAttributes -> {Listable},
Parallelization -> True
][cells];
rp = AccumulateIntegerList[lens];
nn = rp[[-1]];
A = SparseArray @@ {Automatic, {m, n}, 0, {1, {rp, Partition[Flatten[cells], 1]}, ConstantArray[1, nn]}}]
,
{}
]
];
getMeshCellAdjacencyMatrix::usage = "";
getMeshCellAdjacencyMatrix[A_?MatrixQ, d_Integer] :=
If[Length[A] > 0,
With[{B = A.A\[Transpose]},
SparseArray[UnitStep[B - DiagonalMatrix[Diagonal[B]] - d]]
],
{}
];
getMeshCellAdjacencyMatrix[Ad10_?MatrixQ, A0d2_?MatrixQ, d1_Integer,
d2_Integer] := If[(Length[Ad10] > 0) && (Length[A0d2] > 0),
With[{B = Ad10.A0d2}, SparseArray[
If[d1 == d2,
UnitStep[B - DiagonalMatrix[Diagonal[B]] - d1],
UnitStep[B - (Min[d1, d2] + 1)]]
]
],
{}
];
MeshCellAdjacencyMatrix::usage = "";
polymesh /: MeshCellAdjacencyMatrix[M_polymesh, 0, 0 _] := SparseArray[
Join[MeshCells[M, 1],
Transpose[Reverse[Transpose[MeshCells[M, 1]]]]] -> 1,
{1, 1} MeshCellCount[M, 0]
];
polymesh /: MeshCellAdjacencyMatrix[M_polymesh, 0 _, d_Integer] :=
With[{cells = MeshCells[M, d]},
If[Length[cells] > 0,
Transpose[getConnectivityMatrix[MeshCellCount[M, 0], MeshCells[M, d]]],
{}
]
];
polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d_Integer, 0 _] :=
With[{A = MeshCellAdjacencyMatrix[M, 0, d]},
If[Length[A] > 0,
Transpose[MeshCellAdjacencyMatrix[M, 0, d]],
{}
]
];
polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d_Integer, d_Integer] :=
getMeshCellAdjacencyMatrix[MeshCellAdjacencyMatrix[M, d, 0], d]
polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d_Integer] :=
MeshCellAdjacencyMatrix[M, d, d];
polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d1_Integer, d2_Integer] :=
Module[{r, m1, m2},
{m1, m2} = MinMax[{d1, d2}];
r = getMeshCellAdjacencyMatrix[
MeshCellAdjacencyMatrix[M, m1, 0],
MeshCellAdjacencyMatrix[M, 0, m2],
m1,
m2
];
If[d1 < d2, r, If[Length[r] > 0, Transpose[r], {}]]
];
CatmullClarkSubdivisionMatrix::usage = "";
polymesh /: CatmullClarkSubdivisionMatrix[M_polymesh, OptionsPattern[{"AverageBoundaryPoints" -> True}]] :=
Module[{avgbndpQ, A02, A01, A10, valences, bplist, edgevalencelist, χbndp, χbndpcomp, χbnde, χbndecomp, belist, A20, A12, vB, eB, pB, n0, n1},
avgbndpQ = OptionValue["AverageBoundaryPoints"];
n0 = MeshCellCount[M, 0];
n1 = MeshCellCount[M, 1];
A02 = MeshCellAdjacencyMatrix[M, 0, 2];
A20 = Transpose[A02];
A01 = MeshCellAdjacencyMatrix[M, 0, 1];
A10 = Transpose[A01];
A12 = getMeshCellAdjacencyMatrix[A10, A02, 1, 2];
valences = N[Total[A01, {2}]];
belist = Random`Private`PositionsOf[Total[A12, {2}], 1];
χbnde = SparseArray[Transpose[{belist}] -> 1, {n1}, 0];
χbndecomp = (1. - Normal[χbnde]);
bplist = Union @@ MeshCells[M, 1][[belist]];
χbndp = SparseArray[Partition[bplist, 1] -> 1, {n0}, 0];
χbndpcomp = (1. - Normal[χbndp]);
pB = A20 SparseArray[1./(Length /@ MeshCells[M, 2])];
eB = (0.5 χbnde) Transpose[χbndp A01] + SparseArray[0.25 χbndecomp] (A10 + A12.pB);
vB = Plus[
SparseArray[χbndpcomp/valences^2] (A02.pB + A01.A10),
DiagonalMatrix[SparseArray[χbndpcomp (1. - 3./valences) + If[avgbndpQ, 0.75, 1.] Normal[χbndp]]]
];
If[avgbndpQ,
vB += (0.125 χbndp) Transpose[χbndp MeshCellAdjacencyMatrix[M, 0, 0]]
];
Join[vB, eB, pB]
];
CatmullClarkSubdivide::usage = "";
polymesh /: CatmullClarkSubdivide[M0_polymesh, OptionsPattern[{
"Subdivisions" -> 1,
"AverageBoundaryPoints" -> True
}]
] :=
Module[{t, M, A},
M = M0;
If[OptionValue["Subdivisions"] > 0,
PrintTemporary["Subdividing..."];
t = AbsoluteTiming[
A = CatmullClarkSubdivisionMatrix[M, "AverageBoundaryPoints" -> OptionValue["AverageBoundaryPoints"]];
M = polymesh[A.MeshCoordinates[M], SubdividedPolygons[M]];
][[1]];
PrintTemporary["Subdivision done. Time elapsed: ", ToString[t]];
];
If[OptionValue["Subdivisions"] > 1,
M = CatmullClarkSubdivide[M,
"Subdivisions" -> OptionValue["Subdivisions"] - 1,
"AverageBoundaryPoints" -> OptionValue["AverageBoundaryPoints"]
]
];
M
];
To my own surprise, Doo-Sabin subdivision is in many ways much easier to implement than Catmull-Clark subdivision. The only real problem I met was to compute the faces created at vertices correctly. The method I use for this is reasonably fast for meshes with not-to-high vertex degrees, but it is not guaranteed that the subdivision of an oriented mesh is again oriented.
a = ConstantArray[1, {3, 3, 3}];
Do[a[[i, j, k]] = 0;, {i, {1, 3}}, {j, {1, 3}}, {k, {1, 3}}];
R = RegionBoundary[ArrayMesh[a]];
M = polymesh[MeshCoordinates[R], MeshCells[R, 2, "Multicells" -> True][[1, 1]]];
MList = NestList[DooSabinSubdivide, M, 4];
GraphicsRow[GridMeshPlot /@ MList, ImageSize -> Full]
An example with nontrivial boundary:
pts = N[{{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, {1, 1, 1}, {-1, 1, 1}}];
polys = {{4, 3, 2, 1},(*{1,2,6,5},*){2, 3, 7, 6}, {8, 7, 3, 4}, {5, 8, 4, 1}, {5, 6, 7, 8}};
M = polymesh[pts, polys];
MList = NestList[DooSabinSubdivide, M, 4];
GraphicsRow[GridMeshPlot /@ MList, ImageSize -> Full]
In addition to the code from my other post, one requires the following. This version also features correct(?) handling of boundaries.
getDooSabinSubdivisionMasks = Compile[{{n, _Real}},
Block[{ω, a},
ω = 2. Pi/n;
a = (n + 5.)/(4. n);
Flatten@Table[If[i == j, a, (3. + 2. Cos[ω (i - j)])/(4. n)], {i, 1, n}, {j, 1, n}]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
getDooSabinCombinatorics = Compile[{{face, _Integer, 1}, {idx, _Integer}},
Flatten[
Table[
{i, Compile`GetElement[face, j]},
{i, idx - Length[face] + 1, idx}, {j, 1, Length[face]}],
1
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
DooSabinSubdivisionMatrix::usage = "";
polymesh /: DooSabinSubdivisionMatrix[M_polymesh] :=
Module[{lens, acc, n0, n2, polys, L, A, A12, bndelist, edgesneighpolys, e1, e2, i, j, newbndplist, bndn1},
n0 = MeshCellCount[M, 0];
n2 = MeshCellCount[M, 2];
polys = MeshCells[M, 2];
lens = ToPack[Length /@ polys];
acc = AccumulateIntegerList[lens];
L = SparseArray[
Rule[
Join @@ getDooSabinCombinatorics[MeshCells[M, 2], Rest[acc]],
Join @@ getDooSabinSubdivisionMasks[Range[Max[lens]]][[lens]]], {acc[[-1]], n0}];
A12 = MeshCellAdjacencyMatrix[M, 1, 2];
bndelist = Random`Private`PositionsOf[Total[A12, {2}], 1];
If[Length[bndelist] > 0,
A = SparseArray @@ {Automatic, {n2, n0}, 0, {1, {acc, Partition[Join @@ polys, 1]}, Join @@ (Most[acc] + Range[lens])}};
edgesneighpolys = A12[[bndelist]]["AdjacencyLists"];
{e1, e2} = Transpose[MeshCells[M, 1][[bndelist]]];
{i, j} = Transpose[Riffle[Transpose[{Flatten[edgesneighpolys], e1}], Transpose[{Flatten[edgesneighpolys], e2}]]];
newbndplist = cExtractIntegerFromSparseMatrix[A["NonzeroValues"], A["RowPointers"], Flatten[A["ColumnIndices"]], 0, i, j];
bndn1 = Length[bndelist];
L[[newbndplist]] =
SparseArray[
Rule[
Transpose[{Join[Range[2 bndn1], Range[2 bndn1]], Join[Riffle[e1, e2], Riffle[e2, e1]]}], Join[ConstantArray[3./4., 2 bndn1],
ConstantArray[1./4., 2 bndn1]]], {2 bndn1, n0}];
];
L
];
getDooSabinEdgeQuads = Compile[{{edge, _Integer, 1}, {polyidx, _Integer, 1}},
{
{Compile`GetElement[polyidx, 1], Compile`GetElement[edge, 1]},
{Compile`GetElement[polyidx, 1], Compile`GetElement[edge, 2]},
{Compile`GetElement[polyidx, 2], Compile`GetElement[edge, 2]},
{Compile`GetElement[polyidx, 2], Compile`GetElement[edge, 1]}
},
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
getDooSabinVertexFaces = Compile[{{dualedges, _Integer, 2}, {idx, _Integer, 1}},
Block[{n, q, p, i, k, c, j},
n = Length[idx];
q = dualedges[[idx]];
p = Table[0, n];
p[[1]] = Compile`GetElement[q, 1, 1];
p[[2]] = c = Compile`GetElement[q, 1, 2];
i = 1;
k = 2;
While[k < n,
j = 1;
While[
And[j < n, Or[i == j, c != Compile`GetElement[q, j, 1] && c != Compile`GetElement[q, j, 2]]],
j++
];
k++;
i = j;
p[[k]] = c = If[c == Compile`GetElement[q, j, 1],
Compile`GetElement[q, j, 2],
Compile`GetElement[q, j, 1]
];
];
p
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
DooSabinSubdividedPolygons::usage = "";
polymesh /: DooSabinSubdividedPolygons[M_polymesh] :=
Module[{lens, acc, polys, edges, edgesneighpolys, n2, n0, A, B, A12, edgevalences, vertexfaces, edgefaces, facefaces, i, j, t, intelist, bndelist, intplist, intedges},
n0 = MeshCellCount[M, 0];
polys = MeshCells[M, 2];
n2 = Length[polys];
lens = ToPack[Length /@ polys];
acc = AccumulateIntegerList[lens];
facefaces = Internal`PartitionRagged[Range[acc[[-1]]], lens];
A12 = MeshCellAdjacencyMatrix[M, 1, 2];
edgevalences = Total[A12, {2}];
intelist = Random`Private`PositionsOf[edgevalences, 2];
bndelist = Random`Private`PositionsOf[edgevalences, 1];
intedges = MeshCells[M, 1][[intelist]];
intplist = Complement[Range[n0], Flatten[MeshCells[M, 1][[bndelist]]]];
edgesneighpolys = A12[[intelist]]["AdjacencyLists"];
{i, j} =
Transpose[Join @@ getDooSabinEdgeQuads[intedges, edgesneighpolys]];
A = SparseArray @@ {Automatic, {n2, n0}, 0, {1, {acc, Partition[Join @@ polys, 1]}, Join @@ (Most[acc] + Range[lens])}};
edgefaces = Partition[
cExtractIntegerFromSparseMatrix[
A["NonzeroValues"],
A["RowPointers"],
Flatten[A["ColumnIndices"]], 0, i, j
],
4
];
B = SparseArray[
Transpose[{Flatten[intedges], Range[2 Length[intelist]]}] -> 1,
{n0, 2 Length[intelist]}
];
vertexfaces = getDooSabinVertexFaces[
ArrayReshape[edgefaces[[All, {1, 4, 2, 3}]], {2 Length[edgefaces], 2}],
B[[intplist]]["AdjacencyLists"]
];
Join[facefaces, edgefaces, vertexfaces]
];
DooSabinSubdivide::usage = "";
polymesh /: DooSabinSubdivide[M0_polymesh, OptionsPattern[{"Subdivisions" -> 1}]
] := Module[{lens, acc, pat, polys, edges, edgesneighpolys, m, n, nn, A, B, A02, valences, vertexfaces, edgefaces, facefaces, i, j, t, M},
PrintTemporary["Subdividing..."];
t = AbsoluteTiming[
If[OptionValue["Subdivisions"] > 0,
M = M0;
M = polymesh[
DooSabinSubdivisionMatrix[M].MeshCoordinates[M],
DooSabinSubdividedPolygons[M]
];
,
M = M0;
]
][[1]];
PrintTemporary["Doo-Sabin subdivision done. Time elapsed: ", ToString[t]];
If[OptionValue["Subdivisions"] > 1,
M = DooSabinSubdivide[M, "Subdivisions" -> OptionValue["Subdivisions"] - 1]
];
M
];