Catmull-Clark and Doo-Sabin Subdivision Implementations

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[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
  ];

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 reasonable 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, 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
   ];