Catmull–Clark subdivision surface

Catmull–Clark subdivision surface
You are encouraged to solve this task according to the task description, using any language you may know.

Implement the Catmull-Clark surface subdivision (description on Wikipedia).

The process works as follow:

1. for each face, a face point is created which is the average of all the points of the face.
2. for each edge, an edge point is created which is the average between the center of the edge and the center of the segment made with the face points of the two adjacent faces.
3. for each vertex point, its coordinates are updated from (new_coords):
   1. the old coordinates (old_coords),
   2. the average of the face points of the faces the point belongs to (avg_face_points),
   3. the average of the centers of edges the point belongs to (avg_mid_edges),
   4. how many faces a point belongs to (n), then use this formula:

     m1 = (n - 3) / n
     m2 = 1 / n
     m3 = 2 / n
     new_coords = (m1 * old_vertex)
                + (m2 * avg_face_points)
                + (m3 * avg_mid_edges)

Then each face is replaced by new faces made with the new points,

• for a triangle face (a,b,c):

   (a, edge_point_ab, face_point, edge_point_ca)
   (b, edge_point_bc, face_point, edge_point_ab)
   (c, edge_point_ca, face_point, edge_point_bc)

• for a quad face (a,b,c,d):

   (a, edge_point_ab, face_point, edge_point_da)
   (b, edge_point_bc, face_point, edge_point_ab)
   (c, edge_point_cd, face_point, edge_point_bc)
   (d, edge_point_da, face_point, edge_point_cd)

This process is relevant when there are no holes in the surface.

When there is a hole, we can detect it as follow:

• an edge is the border of a hole if it belongs to only one face,
• a point is on the border of a hole if n1 != n2 with n1 the number of faces the point belongs to, and n2 the number of edges a point belongs to.

On the border of a hole the subdivision occurs as follow:

1. for the edges that are on the border of a hole, the edge point is just the middle of the edge.
2. for the vertex points that are on the border of a hole, the new coordinates are calculated as follow:
   1. in all the edges the point belongs to, only take in account the middles of the edges that are on the border of the hole then calculate the average between these points and the old coordinates.

Mathematica

This implementation supports tris, quads, and higher polys, as well as surfaces with holes.

<lang Mathematica>CatmullClark[{v_, i_}] := Block[{e, vc, fp, ep, vp},

 e = Function[a, {a, Select[Transpose[{i, Range@Length@i}], 
       Length@Intersection[#1, a] == 2 &]All, 2}] /@ 
   Union[Sort /@ Flatten[Partition[#, 2, 1, 1] & /@ i, 1]]; 
 vc = Table[{n, Select[Transpose[{i, Range@Length@i}], MemberQ[#1, n] &]All, 2, 
    Select[Transpose[{eAll, 1, Range@Length@eAll, 1}], MemberQ[#1, n] &]All, 2}, {n, Length@v}]; 
 fp = Mean[v#] & /@ i;
 ep = If[Length[#2] == 1, Mean[v[[#1]]], Mean@Join[v[[#1]], fp[[#2]]]] & /@ e; 
 vp = If[Length[#2] != Length[#3], Mean@Join[{v[[#1]]}, ep[[Select[#3, 
          Length[e#, 2] != 2 &]]]], ((Length@#2 - 3) v[[#1]] + Mean@fp[[#2]] +
           2 Mean@ep[[#3]])/Length@#2] & /@ vc;
 {Join[vp, ep, fp], Flatten[Function[a, Function[
       b, {a1, #1 + Length[vc], b + Length[vc] + Length[e], #2 + Length[vc]} &@
        Sort[Select[Transpose[{e, Range@Length@e}], MemberQ[#1, 1, a1] && MemberQ[#1, 2, b] &], 
          With[{f = i[[Intersection[#1, 2, #21, 2]1]], 
             n = Intersection[#1, 1, #21, 1]1}, 
            Xor[Abs[#] == 1, # < 0] &@(Position[f, Complement[#1, 1, {n}]1] - 
                Position[f, n])1, 1] &]All, 2] /@ a2] /@ vc, 1]}]

v = PolyhedronData["Cube", "VertexCoordinates"] // N i = PolyhedronData["Cube", "FaceIndices"] NestList[CatmullClark, {v, i}, 4]; Graphics3D[{FaceForm[{Opacity[0.3]}, {Opacity[0.1]}], GraphicsComplex[#1, Polygon[#2]]}] & /@ % Graphics3D[{EdgeForm[], FaceForm[White, Black],

   GraphicsComplex[#1, Polygon[#2], VertexNormals -> #1]}, Boxed -> False] & /@ %%

</lang>

The last few lines, after the function definition, do a test by using the built-in polyhedron data to generate the vertices and face indices. Then it repeatedly applies the method and graphs the results. Note that this was written in Mathematica 7, although it should be easy enough to port to maybe v5.2.

OCaml

The implementation below only supports quad faces, but it does handle surfaces with holes.

This code uses a module called Dynar (for dynamic array) because it needs a structure similar to arrays but with which we can push a new element at the end. (The source of this module is given in the sub-page.)

In the sub-page there is also a program in OCaml+OpenGL which displays a cube subdivided 2 times with this algorythm.

<lang ocaml>open Dynar

let add3 (x1, y1, z1) (x2, y2, z2) (x3, y3, z3) =

 ( (x1 +. x2 +. x3),
   (y1 +. y2 +. y3),
   (z1 +. z2 +. z3) )

let mul m (x,y,z) = (m *. x, m *. y, m *. z)

let avg pts =

 let n, (x,y,z) =
   List.fold_left
     (fun (n, (xt,yt,zt)) (xi,yi,zi) ->
        succ n, (xt +. xi, yt +. yi, zt +. zi))
     (1, List.hd pts) (List.tl pts)
 in
 let n = float_of_int n in
 (x /. n, y /. n, z /. n)

let catmull ~points ~faces =

 let da_points = Dynar.of_array points in
 let new_faces = Dynar.of_array [| |] in
 let push_face face = Dynar.push new_faces face in
 let h1 = Hashtbl.create 43 in
 let h2 = Hashtbl.create 43 in
 let h3 = Hashtbl.create 43 in
 let h4 = Hashtbl.create 43 in
 let blg = Array.make (Array.length points) 0 in (* how many faces a point belongs to *)
 let f_incr p = blg.(p) <- succ blg.(p) in
 let eblg = Array.make (Array.length points) 0 in (* how many edges a point belongs to *)
 let e_incr p = eblg.(p) <- succ eblg.(p) in
 let edge a b = (min a b, max a b) in  (* suitable for hash-table keys *)
 let mid_edge p1 p2 =
   let x1, y1, z1 = points.(p1)
   and x2, y2, z2 = points.(p2) in
   ( (x1 +. x2) /. 2.0,
     (y1 +. y2) /. 2.0,
     (z1 +. z2) /. 2.0 )
 in
 let mid_face p1 p2 p3 p4 =
   let x1, y1, z1 = points.(p1)
   and x2, y2, z2 = points.(p2)
   and x3, y3, z3 = points.(p3)
   and x4, y4, z4 = points.(p4) in
   ( (x1 +. x2 +. x3 +. x4) /. 4.0,
     (y1 +. y2 +. y3 +. y4) /. 4.0,
     (z1 +. z2 +. z3 +. z4) /. 4.0 )
 in
 Array.iteri (fun i (a,b,c,d) ->
   f_incr a; f_incr b; f_incr c; f_incr d;

   let face_point = mid_face a b c d in
   let face_pi = pushi da_points face_point in
   Hashtbl.add h3 a face_point;
   Hashtbl.add h3 b face_point;
   Hashtbl.add h3 c face_point;
   Hashtbl.add h3 d face_point;

   let process_edge a b =
     let ab = edge a b in
     if not(Hashtbl.mem h1 ab)
     then begin
       let mid_ab = mid_edge a b in
       let index = pushi da_points mid_ab in
       Hashtbl.add h1 ab (index, mid_ab, [face_point]);
       Hashtbl.add h2 a mid_ab;
       Hashtbl.add h2 b mid_ab;
       Hashtbl.add h4 mid_ab 1;
       (index)
     end
     else begin
       let index, mid_ab, fpl = Hashtbl.find h1 ab in
       Hashtbl.replace h1 ab (index, mid_ab, face_point::fpl);
       Hashtbl.add h4 mid_ab (succ(Hashtbl.find h4 mid_ab));
       (index)
     end
   in
   let mid_ab = process_edge a b
   and mid_bc = process_edge b c
   and mid_cd = process_edge c d
   and mid_da = process_edge d a in

   push_face (a, mid_ab, face_pi, mid_da);
   push_face (b, mid_bc, face_pi, mid_ab);
   push_face (c, mid_cd, face_pi, mid_bc);
   push_face (d, mid_da, face_pi, mid_cd);
 ) faces;

 Hashtbl.iter (fun (a,b) (index, mid_ab, fpl) ->
   e_incr a; e_incr b;
   if List.length fpl = 2 then
     da_points.ar.(index) <- avg (mid_ab::fpl)
 ) h1;

 Array.iteri (fun i old_vertex ->
   let n = blg.(i)
   and e_n = eblg.(i) in
   (* if the vertex doesn't belongs to as many faces than edges
      this means that this is a hole *)
   if n = e_n then
   begin
     let avg_face_points =
       let face_point_list = Hashtbl.find_all h3 i in
       (avg face_point_list)
     in
     let avg_mid_edges = 
       let mid_edge_list = Hashtbl.find_all h2 i in
       (avg mid_edge_list)
     in
     let n = float_of_int n in
     let m1 = (n -. 3.0) /. n
     and m2 = 1.0 /. n
     and m3 = 2.0 /. n in
     da_points.ar.(i) <-
         add3 (mul m1 old_vertex)
              (mul m2 avg_face_points)
              (mul m3 avg_mid_edges)
   end
   else begin
     let mid_edge_list = Hashtbl.find_all h2 i in
     let mid_edge_list =
       (* only average mid-edges near the hole *)
       List.fold_left (fun acc mid_edge ->
         match Hashtbl.find h4 mid_edge with
         | 1 -> mid_edge::acc
         | _ -> acc
       ) [] mid_edge_list
     in
     da_points.ar.(i) <- avg (old_vertex :: mid_edge_list)
   end
 ) points;

 (Dynar.to_array da_points,
  Dynar.to_array new_faces)

</lang>

Tcl

This code does not handle holes, but does handle arbitrary polygons in the input data. <lang tcl>package require Tcl 8.5

1. Use math functions and operators as commands (Lisp-like)

namespace path {tcl::mathfunc tcl::mathop}

1. Add 3 points

proc add3 {A B C} {

   lassign $A Ax Ay Az
   lassign $B Bx By Bz
   lassign $C Cx Cy Cz
   list [+ $Ax $Bx $Cx] [+ $Ay $By $Cy] [+ $Az $Bz $Cz]

}

1. Multiply a point by a constant

proc mulC {m A} {

   lassign $A x y z
   list [* $m $x] [* $m $y] [* $m $z]

}

1. Take the average/centroid of a set of points
2. Note that each of the arguments is a *list* of coordinate triples
3. This makes things easier later.

proc avg args {

   set x [set y [set z 0.0]]
   foreach plist $args {

incr n [llength $plist] foreach p $plist { lassign $p px py pz set x [+ $x $px] set y [+ $y $py] set z [+ $z $pz] }

   }
   set n [double $n]
   list [/ $x $n] [/ $y $n] [/ $z $n]

}

1. Select from the list the value from each of the indices in the *lists*
2. in the trailing arguments

proc selectFrom {list args} {

   foreach is $args {foreach i $is {lappend r [lindex $list $i]}}
   return $r

}

1. Rotate a list (only does one direction, only by small amounts)

proc lrot {list {n 1}} {

   list {*}[lrange $list $n end] {*}[lrange $list 0 [incr n -1]]

}

1. Generate an edge by putting the smaller coordinate index first.

proc edge {a b} {if {$a<$b} {list $a $b} {list $b $a}}

1. Perform one step of Catmull-Clark subdivision of a surface.

proc CatmullClark {points faces} {

   # Generate the face-points and list of edges, plus some lookup tables
   set edges {}
   set fnum -1
   foreach f $faces {

incr fnum lappend facepoints [avg [selectFrom $points $f]] foreach p [selectFrom $points $f] { lappend fp4p($p) $fnum } foreach p1 $f p2 [lrot $f] { set e [edge $p1 $p2] lappend edges $e lappend f4e($e) $fnum lappend e4f($fnum) $e }

   }

   # Generate the edge-points and let it also be looked up by point
   set enum -1
   foreach e $edges {

incr enum set ep [avg [selectFrom $points $e] [selectFrom $facepoints $f4e($e)]] lappend edgepoints $ep set en4e($e) $enum foreach p [selectFrom $points $e] { lappend ep4p($p) $enum }

   }

   # Generate the new vertex points
   foreach p $points {

set n [llength $fp4p($p)] set avgFace [avg [selectFrom $facepoints $fp4p($p)]] set avgEdge [avg [selectFrom $points {*}[selectFrom $edges $ep4p($p)]]] lappend newPoints [add3 [mulC [/ [- $n 3.0] $n] $p] \ [mulC [/ 1.0 $n] $avgFace] \ [mulC [/ 2.0 $n] $avgEdge]]

   }

   # Now compute the new set of quadrilateral faces
   set fStart [set eStart [llength $newPoints]]
   incr eStart [llength $facepoints]
   set fidx -1
   foreach f $faces {

incr fidx foreach a $f b [lrot $f] c [lrot $f 2] { lappend newFaces \ [list $a \ [+ $eStart $en4e([edge $a $b])] \ [+ $fStart $fidx] \ [+ $eStart $en4e([edge $b $c])]] }

   }

   return [list [concat $newPoints $facepoints $edgepoints] $newFaces]

}</lang>