Posted to tcl by kbk at Thu Dec 16 02:26:28 GMT 2010view pretty
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