Catmull–Clark subdivision surface: Difference between revisions

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:point) =

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: point list) = 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

| [] -> failwith "edgeless face"

| a :: [] -> f a (List.hd pl)

| a :: b :: c -> f a b; next (b::c) in

next pl;;

module EdgeSet = Set.Make(struct type t = edge let compare = ecmp end)

let catmull_clark faces =

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 "facet") 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 c000 = { x=(-1.0); y=(-1.0); z=(-1.0) }

let c001 = { x=(-1.0); y=(-1.0); z= 1.0 }

let c010 = { x=(-1.0); y= 1.0 ; z=(-1.0) }

let c011 = { x=(-1.0); y= 1.0 ; z= 1.0 }

let c100 = { x= 1.0 ; y=(-1.0); z=(-1.0) }

let c101 = { x= 1.0 ; y=(-1.0); z= 1.0 }

let c110 = { x= 1.0 ; y= 1.0 ; z=(-1.0) }

let c111 = { x= 1.0 ; y= 1.0 ; z= 1.0 }

let cube = [

(Face [c000;c001;c011;c010]); (Face [c100;c101;c111;c110]);

(Face [c000;c100;c101;c001]); (Face [c010;c110;c111;c011]);

(Face [c000;c010;c110;c100]); (Face [c001;c011;c111;c101]) ];;

show_faces cube;;

show_faces (catmull_clark cube);;</lang>

with output:

<pre>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)

}</pre>