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), which is an algorithm that maps from a surface (described as a set of points and a set of polygons with vertices at those points) to another more refined surface. The resulting surface will always consist of a mesh of quadrilaterals.
The process for computing the new locations of the points works as follows when the surface is free of holes:
- for each face, a face point is created which is the average of all the points of the face.
- 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.
- for each vertex point, its coordinates are updated from (new_coords):
- the old coordinates (old_coords),
- the average of the face points of the faces the point belongs to (avg_face_points),
- the average of the centers of edges the point belongs to (avg_mid_edges),
- 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_coords) + (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_pointab, face_pointabc, edge_pointca) (b, edge_pointbc, face_pointabc, edge_pointab) (c, edge_pointca, face_pointabc, edge_pointbc)
- for a quad face (a,b,c,d):
(a, edge_pointab, face_pointabcd, edge_pointda) (b, edge_pointbc, face_pointabcd, edge_pointab) (c, edge_pointcd, face_pointabcd, edge_pointbc) (d, edge_pointda, face_pointabcd, edge_pointcd)
When there is a hole, we can detect it as follows:
- 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 nfaces != nedges with nfaces the number of faces the point belongs to, and nedges the number of edges a point belongs to.
On the border of a hole the subdivision occurs as follows:
- for the edges that are on the border of a hole, the edge point is just the middle of the edge.
- for the vertex points that are on the border of a hole, the new coordinates are calculated as follows:
- 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
- calculate the average between these points (on the hole boundary) and the old coordinates (also on the hole boundary).
For edges and vertices not next to a hole, the standard algorithm from above is used.
C
Only the subdivision part. The full source is way too long to be shown here. Lots of macros, you'll have to see the full code to know what's what. <lang c>vertex face_point(face f) { int i; vertex v;
if (!f->avg) { f->avg = vertex_new(); foreach(i, v, f->v) if (!i) f->avg->pos = v->pos; else vadd(f->avg->pos, v->pos);
vdiv(f->avg->pos, len(f->v)); } return f->avg; }
- define hole_edge(e) (len(e->f)==1)
vertex edge_point(edge e) { int i; face f;
if (!e->e_pt) { e->e_pt = vertex_new(); e->avg = e->v[0]->pos; vadd(e->avg, e->v[1]->pos); e->e_pt->pos = e->avg;
if (!hole_edge(e)) { foreach (i, f, e->f) vadd(e->e_pt->pos, face_point(f)->pos); vdiv(e->e_pt->pos, 4); } else vdiv(e->e_pt->pos, 2);
vdiv(e->avg, 2); }
return e->e_pt; }
- define hole_vertex(v) (len((v)->f) != len((v)->e))
vertex updated_point(vertex v) { int i, n = 0; edge e; face f; coord_t sum = {0, 0, 0};
if (v->v_new) return v->v_new;
v->v_new = vertex_new(); if (hole_vertex(v)) { v->v_new->pos = v->pos; foreach(i, e, v->e) { if (!hole_edge(e)) continue; vadd(v->v_new->pos, edge_point(e)->pos); n++; } vdiv(v->v_new->pos, n + 1); } else { n = len(v->f); foreach(i, f, v->f) vadd(sum, face_point(f)->pos); foreach(i, e, v->e) vmadd(sum, edge_point(e)->pos, 2, sum); vdiv(sum, n); vmadd(sum, v->pos, n - 3, sum); vdiv(sum, n); v->v_new->pos = sum; }
return v->v_new; }
model catmull(model m) { int i, j, a, b, c, d; face f; vertex v, x;
model nm = model_new(); foreach (i, f, m->f) { foreach(j, v, f->v) { _get_idx(a, updated_point(v)); _get_idx(b, edge_point(elem(f->e, (j + 1) % len(f->e)))); _get_idx(c, face_point(f)); _get_idx(d, edge_point(elem(f->e, j))); model_add_face(nm, 4, a, b, c, d); } } return nm; }</lang>
J
<lang j>havePoints=: e."1/~ i.@# catmullclark=:3 :0
'mesh points'=. y face_point=. avg"2 mesh{points point_face=. |: mesh havePoints points avg_face_points=. point_face avg@#"1 2 face_point edges=. ~.,/ meshEdges=. mesh /:~@,"+1|."1 mesh edge_face=. *./"2 edges e."0 1/ mesh edge_center=. avg"2 edges{points edge_point=. (0.5*edge_center) + 0.25 * edge_face +/ .* face_point point_edge=. |: edges havePoints points avg_mid_edges=. point_edge avg@#"1 2 edge_center n=. +/"1 point_edge 'm3 m2 m1'=. (2,1,:n-3)%"1 n new_coords=. (m1 * points) + (m2 * avg_face_points) + (m3 * avg_mid_edges) pts=. face_point,edge_point,new_coords c0=. (#edge_point)+ e0=. #face_point msh=. (,c0+mesh),.(,e0+edge i. meshEdges),.((#i.)~/$mesh),.,e0+_1|."1 edge i. meshEdges msh;pts
)</lang>
Example use:
<lang j>NB.cube points=: _1+2*#:i.8 mesh=: 1 A."1 I.(,1-|.)8&$@#&0 1">4 2 1
catmullclark mesh;points
┌──────────┬─────────────────────────────┐ │22 6 0 9│ 1 0 0│ │23 7 0 6│ 0 1 0│ │25 8 0 7│ 0 0 1│ │24 9 0 8│ 0 0 _1│ │20 10 1 12│ 0 _1 0│ │21 11 1 10│ _1 0 0│ │25 8 1 11│ 0.75 _0.75 0│ │24 12 1 8│ 0.75 0 0.75│ │19 13 2 14│ 0.75 0.75 0│ │21 11 2 13│ 0.75 0 _0.75│ │25 7 2 11│ _0.75 0.75 0│ │23 14 2 7│ 0 0.75 0.75│ │18 15 3 16│ 0 0.75 _0.75│ │20 12 3 15│ _0.75 0 0.75│ │24 9 3 12│ 0 _0.75 0.75│ │22 16 3 9│ _0.75 0 _0.75│ │18 17 4 16│ 0 _0.75 _0.75│ │19 14 4 17│ _0.75 _0.75 0│ │23 6 4 14│_0.555556 _0.555556 _0.555556│ │22 16 4 6│_0.555556 _0.555556 0.555556│ │18 17 5 15│_0.555556 0.555556 _0.555556│ │19 13 5 17│_0.555556 0.555556 0.555556│ │21 10 5 13│ 0.555556 _0.555556 _0.555556│ │20 15 5 10│ 0.555556 _0.555556 0.555556│ │ │ 0.555556 0.555556 _0.555556│ │ │ 0.555556 0.555556 0.555556│ └──────────┴─────────────────────────────┘</lang>
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 algorithm.
<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>
Another implementation
Another implementation which should work with holes, but has only been tested on a cube <lang OCaml>type point = { x: float; y : float; z : float } let add a b = { x = a.x+.b.x; y = a.y+.b.y; z = a.z+.b.z } let mul a k = { x = a.x*.k; y = a.y*.k; z= a.z*.k } let div p k = mul p (1.0/.k) let fsgn x y = if x < y then -1 else if x > y then 1 else 0 let cmp a b = if a.x=b.x then if a.y=b.y then fsgn b.z a.z else fsgn b.y a.y else fsgn b.x a.x
type face = Face of point list type edge = Edge of point*point let ecmp (Edge (p1,p2)) (Edge (p3,p4)) = let sgn = cmp p1 p3 in if sgn = 0 then cmp p2 p4 else sgn
let make_edge a b = if cmp a b < 0 then Edge (b,a) else Edge (a,b) let make_face a b c d = Face [a;b;c;d]
let centroid plist = div (List.fold_left add {x=0.0;y=0.0;z=0.0} plist) (float (List.length plist)) let mid_edge (Edge (p1,p2)) = div (add p1 p2) 2.0 let face_point (Face pl) = centroid pl let point_in_face p (Face pl) = List.mem p pl let point_in_edge p (Edge (p1,p2)) = (p = p1 || p = p2) let edge_in_face (Edge (p1,p2)) (Face pl) = (List.mem p1 pl && List.mem p2 pl)
let border_edge faces e =
List.length (List.filter (edge_in_face e) faces) < 2
let edge_point faces e =
if border_edge faces e then mid_edge e else let adjacent = List.filter (edge_in_face e) faces in let fps = List.map face_point adjacent in centroid [mid_edge e; centroid fps]
let mod_vertex faces edges p =
let v_edges = List.filter (point_in_edge p) edges in let v_faces = List.filter (point_in_face p) faces in let n = List.length v_faces in let is_border = n != (List.length v_edges) in if is_border then let border_mids = List.map mid_edge (List.filter (border_edge faces) v_edges) in (* description ambiguity: average (border+p) or average(average(border),p) ?? *) centroid (p :: border_mids) else let avg_face = centroid (List.map face_point v_faces) in let avg_mid = centroid (List.map mid_edge v_edges) in div (add (add (mul p (float(n-3))) avg_face) (mul avg_mid 2.0)) (float n)
let iter_edges f (Face pl) =
let rec next = function | [] -> () | a :: [] -> f a (List.hd pl) | a :: b :: c -> f a b; next (b::c) in next pl;;
let catmull_clark faces =
let module EdgeSet = Set.Make(struct type t = edge let compare = ecmp end) in let eset = ref EdgeSet.empty in let add_edge a b = eset := EdgeSet.add (make_edge a b) !eset in let edges = (List.iter (iter_edges add_edge) faces; EdgeSet.elements !eset) in let new_faces = ref [] in let mod_face ((Face pl) as face) = let fp = face_point face in let ep = ref [] in ( iter_edges (fun a b -> ep := (edge_point faces (make_edge a b)):: !ep) face; let e_tl = List.hd (List.rev !ep) in let v' = List.map (mod_vertex faces edges) pl in let rec add_facet e vl el = (match (vl,el) with | (h1::t1),(h2::t2) -> new_faces := (make_face e h1 h2 fp) :: !new_faces; add_facet h2 t1 t2 | ([],[]) -> () | _ -> failwith "vertex/edge mismatch") in add_facet e_tl v' !ep) in (List.iter mod_face faces; !new_faces)
let show_faces fl =
let pr_point p = Printf.printf " (%.4f, %.4f, %.4f)" p.x p.y p.z in let pr_face (Face(pl)) = print_string "Face:"; List.iter pr_point pl; print_string "\n" in (print_string "surface {\n"; List.iter pr_face fl; print_string "}\n")
let c p q r = let s i = if i = 0 then -1.0 else 1.0 in { x = s p; y = s q; z = s r } ;; let cube = [
Face [c 0 0 0; c 0 0 1; c 0 1 1; c 0 1 0]; Face [c 1 0 0; c 1 0 1; c 1 1 1; c 1 1 0]; Face [c 0 0 0; c 1 0 0; c 1 0 1; c 0 0 1]; Face [c 0 1 0; c 1 1 0; c 1 1 1; c 0 1 1]; Face [c 0 0 0; c 0 1 0; c 1 1 0; c 1 0 0]; Face [c 0 0 1; c 0 1 1; c 1 1 1; c 1 0 1] ] in
show_faces cube; show_faces (catmull_clark cube)</lang> with output:
surface { Face: (-1.0000, -1.0000, -1.0000) (-1.0000, -1.0000, 1.0000) (-1.0000, 1.0000, 1.0000) (-1.0000, 1.0000, -1.0000) Face: (1.0000, -1.0000, -1.0000) (1.0000, -1.0000, 1.0000) (1.0000, 1.0000, 1.0000) (1.0000, 1.0000, -1.0000) Face: (-1.0000, -1.0000, -1.0000) (1.0000, -1.0000, -1.0000) (1.0000, -1.0000, 1.0000) (-1.0000, -1.0000, 1.0000) Face: (-1.0000, 1.0000, -1.0000) (1.0000, 1.0000, -1.0000) (1.0000, 1.0000, 1.0000) (-1.0000, 1.0000, 1.0000) Face: (-1.0000, -1.0000, -1.0000) (-1.0000, 1.0000, -1.0000) (1.0000, 1.0000, -1.0000) (1.0000, -1.0000, -1.0000) Face: (-1.0000, -1.0000, 1.0000) (-1.0000, 1.0000, 1.0000) (1.0000, 1.0000, 1.0000) (1.0000, -1.0000, 1.0000) } surface { Face: (0.0000, 0.7500, 0.7500) (0.5556, -0.5556, 0.5556) (-0.7500, 0.0000, 0.7500) (0.0000, 0.0000, 1.0000) Face: (0.7500, 0.0000, 0.7500) (0.5556, 0.5556, 0.5556) (0.0000, 0.7500, 0.7500) (0.0000, 0.0000, 1.0000) Face: (0.0000, -0.7500, 0.7500) (-0.5556, 0.5556, 0.5556) (0.7500, 0.0000, 0.7500) (0.0000, 0.0000, 1.0000) Face: (-0.7500, 0.0000, 0.7500) (-0.5556, -0.5556, 0.5556) (0.0000, -0.7500, 0.7500) (0.0000, 0.0000, 1.0000) Face: (0.0000, 0.7500, -0.7500) (0.5556, -0.5556, -0.5556) (-0.7500, 0.0000, -0.7500) (0.0000, 0.0000, -1.0000) Face: (0.7500, 0.0000, -0.7500) (0.5556, 0.5556, -0.5556) (0.0000, 0.7500, -0.7500) (0.0000, 0.0000, -1.0000) Face: (0.0000, -0.7500, -0.7500) (-0.5556, 0.5556, -0.5556) (0.7500, 0.0000, -0.7500) (0.0000, 0.0000, -1.0000) Face: (-0.7500, 0.0000, -0.7500) (-0.5556, -0.5556, -0.5556) (0.0000, -0.7500, -0.7500) (0.0000, 0.0000, -1.0000) Face: (0.7500, 0.7500, 0.0000) (-0.5556, 0.5556, 0.5556) (0.0000, 0.7500, -0.7500) (0.0000, 1.0000, 0.0000) Face: (0.0000, 0.7500, 0.7500) (0.5556, 0.5556, 0.5556) (0.7500, 0.7500, 0.0000) (0.0000, 1.0000, 0.0000) Face: (-0.7500, 0.7500, 0.0000) (0.5556, 0.5556, -0.5556) (0.0000, 0.7500, 0.7500) (0.0000, 1.0000, 0.0000) Face: (0.0000, 0.7500, -0.7500) (-0.5556, 0.5556, -0.5556) (-0.7500, 0.7500, 0.0000) (0.0000, 1.0000, 0.0000) Face: (0.7500, -0.7500, 0.0000) (-0.5556, -0.5556, 0.5556) (0.0000, -0.7500, -0.7500) (0.0000, -1.0000, 0.0000) Face: (0.0000, -0.7500, 0.7500) (0.5556, -0.5556, 0.5556) (0.7500, -0.7500, 0.0000) (0.0000, -1.0000, 0.0000) Face: (-0.7500, -0.7500, 0.0000) (0.5556, -0.5556, -0.5556) (0.0000, -0.7500, 0.7500) (0.0000, -1.0000, 0.0000) Face: (0.0000, -0.7500, -0.7500) (-0.5556, -0.5556, -0.5556) (-0.7500, -0.7500, 0.0000) (0.0000, -1.0000, 0.0000) Face: (0.7500, 0.0000, 0.7500) (0.5556, 0.5556, -0.5556) (0.7500, -0.7500, 0.0000) (1.0000, 0.0000, 0.0000) Face: (0.7500, 0.7500, 0.0000) (0.5556, 0.5556, 0.5556) (0.7500, 0.0000, 0.7500) (1.0000, 0.0000, 0.0000) Face: (0.7500, 0.0000, -0.7500) (0.5556, -0.5556, 0.5556) (0.7500, 0.7500, 0.0000) (1.0000, 0.0000, 0.0000) Face: (0.7500, -0.7500, 0.0000) (0.5556, -0.5556, -0.5556) (0.7500, 0.0000, -0.7500) (1.0000, 0.0000, 0.0000) Face: (-0.7500, 0.0000, 0.7500) (-0.5556, 0.5556, -0.5556) (-0.7500, -0.7500, 0.0000) (-1.0000, 0.0000, 0.0000) Face: (-0.7500, 0.7500, 0.0000) (-0.5556, 0.5556, 0.5556) (-0.7500, 0.0000, 0.7500) (-1.0000, 0.0000, 0.0000) Face: (-0.7500, 0.0000, -0.7500) (-0.5556, -0.5556, 0.5556) (-0.7500, 0.7500, 0.0000) (-1.0000, 0.0000, 0.0000) Face: (-0.7500, -0.7500, 0.0000) (-0.5556, -0.5556, -0.5556) (-0.7500, 0.0000, -0.7500) (-1.0000, 0.0000, 0.0000) }
Tcl
This code handles both holes and arbitrary polygons in the input data. <lang tcl>package require Tcl 8.5
- Use math functions and operators as commands (Lisp-like).
namespace path {tcl::mathfunc tcl::mathop}
- 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]
}
- Multiply a point by a constant.
proc mulC {m A} {
lassign $A x y z list [* $m $x] [* $m $y] [* $m $z]
}
- Take the centroid of a set of points.
- Note that each of the arguments is a *list* of coordinate triples
- This makes things easier later.
proc centroid 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]
}
- Select from the list the value from each of the indices in the *lists*
- in the trailing arguments.
proc selectFrom {list args} {
foreach is $args {foreach i $is {lappend r [lindex $list $i]}} return $r
}
- Rotate a list.
proc lrot {list {n 1}} {
set n [% $n [llength $list]] list {*}[lrange $list $n end] {*}[lrange $list 0 [incr n -1]]
}
- Generate an edge by putting the smaller coordinate index first.
proc edge {a b} {
list [min $a $b] [max $a $b]
}
- Perform one step of Catmull-Clark subdivision of a surface.
proc CatmullClark {points faces} {
# Generate the new face-points and list of edges, plus some lookup tables. set edges {} foreach f $faces {
set ps [selectFrom $points $f] set fp [centroid $ps] lappend facepoints $fp foreach p $ps { lappend fp4p($p) $fp } foreach p1 $f p2 [lrot $f] { set e [edge $p1 $p2] if {$e ni $edges} { lappend edges $e } lappend fp4e($e) $fp }
}
# Generate the new edge-points and mid-points of edges, and a few more # lookup tables. set i [+ [llength $points] [llength $faces]] foreach e $edges {
set ep [selectFrom $points $e] if {[llength $fp4e($e)] > 1} { set mid [centroid $ep $fp4e($e)] } else { set mid [centroid $ep] foreach p $ep { lappend ep_heavy($p) $mid } } lappend edgepoints $mid set en4e($e) $i foreach p $ep { lappend ep4p($p) $mid } incr i
}
# Generate the new vertex points with our lookup tables. foreach p $points {
if {[llength $fp4p($p)] >= 4} { set n [llength $fp4p($p)] lappend newPoints [add3 [mulC [/ [- $n 3.0] $n] $p] \ [mulC [/ 1.0 $n] [centroid $fp4p($p)]] \ [mulC [/ 2.0 $n] [centroid $ep4p($p)]]] } else { # Update a point on the edge of a hole. This formula is not # described on the WP page, but produces a nice result. lappend newPoints [centroid $ep_heavy($p) [list $p $p]] }
}
# Now compute the new set of quadrilateral faces. set i [llength $points] foreach f $faces {
foreach a $f b [lrot $f] c [lrot $f -1] { lappend newFaces [list \ $a $en4e([edge $a $b]) $i $en4e([edge $c $a])] } incr i
}
list [concat $newPoints $facepoints $edgepoints] $newFaces
}</lang>
The test code for this solution is available as well. The example there produces the following partial toroid output image: