# Catmull-Clark and Doo-Sabin Subdivision Implementations

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.

• There is an implementation of Catmull-Clark on Rosetta Code Apr 18 '19 at 20:25

# Catmull-Clark Subdivision

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.

## Application

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:

## Implementation

getEdgesFromPolygons = Compile[{{f, _Integer, 1}},
Table[
{
Min[CompileGetElement[f, i], CompileGetElement[f, Mod[i + 1, Length[f], 1]]],
Max[CompileGetElement[f, i], CompileGetElement[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[
{
CompileGetElement[qq, i],
CompileGetElement[ee, i],
n,
CompileGetElement[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 += CompileGetElement[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 = CompileGetElement[rp, i] + 1;
c = CompileGetElement[rp, i + 1];
While[k < c + 1 && CompileGetElement[ci, k] != j, ++k];
If[k == c + 1, background, CompileGetElement[vals, k]]
],
RuntimeAttributes -> {Listable},
Parallelization -> True,
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];

ToPack = DeveloperToPackedArray;

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 = SparseArraySparseArraySort@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]]
]];
InternalPartitionRagged[
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]}}]
,
{}
]
];

If[Length[A] > 0,
With[{B = A.A\[Transpose]},
SparseArray[UnitStep[B - DiagonalMatrix[Diagonal[B]] - d]]
],
{}
];

d2_Integer] := If[(Length[Ad10] > 0) && (Length[A0d2] > 0),
If[d1 == d2,
UnitStep[B - DiagonalMatrix[Diagonal[B]] - d1],
UnitStep[B - (Min[d1, d2] + 1)]]
]
],
{}
];

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 _] :=
If[Length[A] > 0,
{}
]
];

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d_Integer, d_Integer] :=

polymesh /: MeshCellAdjacencyMatrix[M_polymesh, d1_Integer, d2_Integer] :=
Module[{r, m1, m2},
{m1, m2} = MinMax[{d1, d2}];
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];
A20 = Transpose[A02];
A10 = Transpose[A01];
A12 = getMeshCellAdjacencyMatrix[A10, A02, 1, 2];
valences = N[Total[A01, {2}]];
belist = RandomPrivatePositionsOf[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
];

• (+1) thanks for posting. One thing I do not get is the purpose of such sub divisions. What are they used for? Apr 18 '19 at 11:04
• Dear Henrik Schumacher, I appreciate your help since this is exactly what i was looking for Catmull-Clark . Anyone else with Doo-Sabin subdivision surface codes can also share them with us. Apr 18 '19 at 11:08
• @user21 Finite elements! ;) The limiting surfaces are proven to be C^2 apart from "exceptional" points where the initial mesh does not have the generic valence (for triangle meshes, the generic valence is 6 and for quad meshes, the generic valuence is 4). So, you get $C^1$-finite elements for surfaces. But the curvature aorund the exceptional vertices blows up so much that subdivision surface are actually not overly useful for minimizing curvature energies... (you have to add many quadrature points around the curvature singularieties. Apr 18 '19 at 11:08
• @user21 The fine thing about the limiting surfaces: You can evaluate also exact positions and normals with finite cost from the initial mesh (the so-called "control mesh"). Apr 18 '19 at 11:10
• @user21 So, subdivision surfaces are not that much used for finite element analysis, although there has been a small hype in the graphics processing community a few years ago. Mostly, subdivision surfaces are used as a flexible tool for surface design quite similar how splines and NURBS surfaces are used: The vertices of the initial mesh provide you with controls to deform the limiting surface in a surprisingly intuitive and predictable way. And you have also a hierarchy of mesh resolutions, so you can generate quick previews and high quality meshes with the same pipeline. Apr 18 '19 at 11:13

# Doo-Sabin Subdivision

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.

## Application

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]


## Implementation

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, CompileGetElement[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]],

bndelist = RandomPrivatePositionsOf[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])}};
{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}},
{
{CompileGetElement[polyidx, 1], CompileGetElement[edge, 1]},
{CompileGetElement[polyidx, 1], CompileGetElement[edge, 2]},
{CompileGetElement[polyidx, 2], CompileGetElement[edge, 2]},
{CompileGetElement[polyidx, 2], CompileGetElement[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]] = CompileGetElement[q, 1, 1];
p[[2]] = c = CompileGetElement[q, 1, 2];
i = 1;
k = 2;
While[k < n,
j = 1;
While[
And[j < n, Or[i == j, c != CompileGetElement[q, j, 1] && c != CompileGetElement[q, j, 2]]],
j++
];
k++;
i = j;
p[[k]] = c = If[c == CompileGetElement[q, j, 1],
CompileGetElement[q, j, 2],
CompileGetElement[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 = InternalPartitionRagged[Range[acc[[-1]]], lens];

edgevalences = Total[A12, {2}];
intelist = RandomPrivatePositionsOf[edgevalences, 2];
bndelist = RandomPrivatePositionsOf[edgevalences, 1];
intedges = MeshCells[M, 1][[intelist]];
intplist = Complement[Range[n0], Flatten[MeshCells[M, 1][[bndelist]]]];
{i, j} =
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}],
];
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
];