Catmull-Clark - full up program with test harness

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]
	set mid [centroid $ep]
	if {[llength $fp4e($e)] > 1} {
	    lappend edgepoints [centroid $ep $fp4e($e)]
	} else {
	    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 {
	set n [llength $fp4p($p)]
	if {$n == [llength $ep4p($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 {
	    lappend newPoints [centroid [list $p] $ep4p($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
}

package require Tk
 
# A simple-minded ordering function for faces
proc orderf {points face1 face2} {
    set d1 [set d2 0.0]
    foreach p [selectFrom $points $face1] {
	lassign $p x y z
	set d1 [expr {$d1 + sqrt($x*$x + $y*$y + $z*$z)}]
    }
    foreach p [selectFrom $points $face2] {
	lassign $p x y z
	set d2 [expr {$d2 + sqrt($x*$x + $y*$y + $z*$z)}]
    }
    expr {$d1<$d2 ? -1 : $d1>$d2 ? 1 : 0}
}
 
# Plots a net defined in points-and-faces fashion
proc visualizeNet {w points faces args} {
    foreach face [lsort -command [list orderf $points] $faces] {
	set c {}
	set polyCoords [selectFrom $points $face]
	set sum {[list 0. 0. 0.]}
	set centroid [centroid $polyCoords]
	foreach coord $polyCoords {
	    lassign $coord x y z
	    lappend c \
		[expr {200. + 190. * (0.867 * $x - 0.9396 * $y)}] \
		[expr {200 + 190. * (0.5 * $x + 0.3402 * $y - $z)}]
	}
	lassign $centroid x y z
	set depth [expr {int(255*sqrt($x*$x + $y*$y + $z*$z) / sqrt(3.))}]
	set grey [format #%02x%02x%02x $depth $depth $depth]
	$w create polygon $c -fill $grey {*}$args
    }
}
# Make a display surface
pack [canvas .c -width 400 -height 400 -background #7f7f7f]
 
# Points to define the unit cube
set points {
    {0.0 0.0 0.0}
    {1.0 0.0 0.0}
    {1.0 1.0 0.0}
    {0.0 1.0 0.0}
    {0.0 0.0 1.0}
    {1.0 0.0 1.0}
    {1.0 1.0 1.0}
    {0.0 1.0 1.0}
}
foreach pt $points {
    lassign $pt x y z
    lappend points [list [expr {0.25 + 0.5*$x}] [expr {0.25 + 0.5*$y}] $z]
}
 
# Try removing {1 2 6 5} to demonstrate holes.
set faces {
    {0 8 9 1}
    {1 9 10 2}
    {2 10 11 3}
    {3 11 8 0}
    {0 1 5 4}
    {1 2 6 5}
    {2 3 7 6}
    {3 0 4 7}
    {4 5 13 12}
    {5 6 14 13}
    {6 7 15 14}
    {7 4 12 15}
    {8 9 13 12}
    {9 10 14 13}
    {10 11 15 14}
    {11 8 12 15}
}
 
# Show the initial layout
visualizeNet .c $points $faces -outline white -fill {}
 
# Apply the Catmull-Clark algorithm to generate a new surface
lassign [CatmullClark $points $faces] points2 faces2
 
## Uncomment the next line to get the second level of subdivision
lassign [CatmullClark $points2 $faces2] points2 faces2
lassign [CatmullClark $points2 $faces2] points2 faces2
 
# Visualize the new surface
visualizeNet .c $points2 $faces2 -outline #0000cc