#!/bin/env tclsh

# polyhedra.tcl 
# Author:      Gerard Sookahet
# Date:        30 Mai 2005 
# Description: Rotating polyhedra using a 'standard' tk canvas.
#              Flat shading and wireframe mode.

package require Tcl 8.5
package require Tk 8.4

bind all <Escape> {exit}

proc Barycenter {lcoords} {
    set X 0
    set Y 0
    set n [llength $lcoords]
  
    foreach vtx $lcoords {
        foreach {x y} $vtx {
            set X [expr {$X + $x}]
            set Y [expr {$Y + $y}]
        }
    }
    return [list [expr {$X/$n}] [expr {$Y/$n}]]
}

proc CrossProduct {x1 y1 z1 x2 y2 z2} {
    return [list [expr {$y1*$z2 - $y2*$z1}] \
                     [expr {$z1*$x2 - $z2*$x1}] \
                     [expr {$x1*$y2 - $x2*$y1}]]
}

proc DotProduct {x1 y1 z1 x2 y2 z2} {
     return [expr {$x1*$x2 + $y1*$y2 + $z1*$z2}]
}

proc MatrixVectorProduct {M V} {
    set x [lindex $V 0]
    set y [lindex $V 1]
    set z [lindex $V 2]
    return [list [expr {[lindex $M 0 0]*$x+[lindex $M 1 0]*$y+[lindex $M 2 0]*$z}] \
                     [expr {[lindex $M 0 1]*$x+[lindex $M 1 1]*$y+[lindex $M 2 1]*$z}] \
                     [expr {[lindex $M 0 2]*$x+[lindex $M 1 2]*$y+[lindex $M 2 2]*$z}]]
}

proc MatrixProduct {M1 M2} {
    set M {{0 0 0 0} {0 0 0 0} {0 0 0 0} {0 0 0 0}}
  
    for {set i 0} {$i<4} {incr i} {
        for {set j 0} {$j<4} {incr j} {
            lset M $i $j 0
            for {set k 0} {$k<4} {incr k} {
                lset M $i $j [expr {[lindex $M $i $j]+[lindex $M1 $i $k]*[lindex $M2 $k $j]}]
            }
        }
    }
    return $M
}

proc MatrixRotation { ax ay az } {
    set sax [expr {sin($ax)}]
    set cax [expr {cos($ax)}]
    set say [expr {sin($ay)}]
    set cay [expr {cos($ay)}]
    set saz [expr {sin($az)}]
    set caz [expr {cos($az)}]
    set Mx {{1 0 0 0} {0 0 0 0} {0 0 0 0} {0 0 0 1}}
    set My {{0 0 0 0} {0 1 0 0} {0 0 0 0} {0 0 0 1}}
    set Mz {{0 0 0 0} {0 0 0 0} {0 0 1 0} {0 0 0 1}}
    
    # Rotation matrix around X axis with angle ax
    lset Mx 1 1 $cax
    lset Mx 1 2 $sax
    lset Mx 2 1 [expr {-1*$sax}]
    lset Mx 2 2 $cax
    # Rotation matrix around Y axis with angle ay
    lset My 0 0 $cay
    lset My 0 2 [expr {-1*$say}]
    lset My 2 0 $say
    lset My 2 2 $cay
    # Rotation matrix around Z axis with angle az
    lset Mz 0 0 $caz
    lset Mz 0 1 $saz
    lset Mz 1 0 [expr {-1*$saz}]
    lset Mz 1 1 $caz
    return [MatrixProduct [MatrixProduct $Mx $My] $Mz]
}

# Compute normal vector and norm for each face
# -------------------------------------------------------------------
proc NormalVector {lvtx lcnx} {
    set lnv {} 
    set lmv {} 
  
    foreach face $lcnx {
        foreach {nx ny nz} [CrossProduct \
        [expr {[lindex $lvtx [lindex $face 1] 0] - [lindex $lvtx [lindex $face 0] 0]}] \
        [expr {[lindex $lvtx [lindex $face 1] 1] - [lindex $lvtx [lindex $face 0] 1]}] \
        [expr {[lindex $lvtx [lindex $face 1] 2] - [lindex $lvtx [lindex $face 0] 2]}] \
        [expr {[lindex $lvtx [lindex $face 2] 0] - [lindex $lvtx [lindex $face 1] 0]}] \
        [expr {[lindex $lvtx [lindex $face 2] 1] - [lindex $lvtx [lindex $face 1] 1]}] \
        [expr {[lindex $lvtx [lindex $face 2] 2] - [lindex $lvtx [lindex $face 1] 2]}]] {}
        lappend lnv [list $nx $ny $nz]
        lappend lmv [DotProduct $nx $ny $nz $nx $ny $nz]
    }
    return [list $lnv $lmv]
}

# 2D projection
# -------------------------------------------------------------------
proc Projection {x y z M} {
    global scx scy vdist
 
    set nx [expr  {[lindex $M 0 0]*$x+[lindex $M 1 0]*$y+[lindex $M 2 0]*$z}]
    set ny [expr  {[lindex $M 0 1]*$x+[lindex $M 1 1]*$y+[lindex $M 2 1]*$z}]
    set nz [expr {([lindex $M 0 2]*$x+[lindex $M 1 2]*$y+[lindex $M 2 2]*$z+10)/$vdist}]
    return [list [expr {$nx/$nz+$scx/2.0}] [expr {$ny/$nz+$scy/2.0}]]
}

# Apply transformations to vertex coordinates
# -------------------------------------------------------------------
proc Transformations {lvtx lnv} {
    global ax ay az 
    set lnew {}
    set lvn {}
    # Compute matrix rotation
    set M [MatrixRotation $ax $ay $az]
    set i 0
    # Apply projection
    foreach vtx $lvtx {
        lappend lnew [Projection [lindex $vtx 0] [lindex $vtx 1] [lindex $vtx 2] $M]
        incr i
    }
    # Normal vector rotation
    foreach v $lnv {lappend lvn [MatrixVectorProduct $M $v]}
  
    return [list $M $lnew $lvn]
}

# Compute color entensity for each face
# -------------------------------------------------------------------
proc Intensity {lnv lmv lvv} {
    set lclr {}
    set v [DotProduct [lindex $lvv 0] [lindex $lvv 1] [lindex $lvv 2] \
                            [lindex $lvv 0] [lindex $lvv 1] [lindex $lvv 2]]
    set i 0
    foreach nv $lnv {
        set clr 31
        set a [DotProduct [lindex $nv 0]  [lindex $nv 1]  [lindex $nv 2] \
                                [lindex $lvv 0] [lindex $lvv 1] [lindex $lvv 2]]
        set b [expr {sqrt([lindex $lmv $i]*$v)}]
        set clr [expr {round(31*($a/$b))}] 
        lappend lclr [expr {$clr < 0 ? 31 : [expr {32 - $clr}]}]
        incr i
    }
    return $lclr
}

# Start the display and rotation loop
# -------------------------------------------------------------------
proc DisplayModel {w s} {
    global stop
    global display
    global ax ay az tx ty tz
    global form
  
    $w.c delete all
    set stop 0
    global iterations
    set ax 0.2
    set ay 0.1
    set az 0.3
    set tx 0
    set ty 0
    set tz 0
  
    set d $display
    foreach {t lvtx lcnx lclr} [ReadData $s] {}
    $w.c create text 10 10 -anchor w -fill white -text $t
    foreach {lnv lmv} [NormalVector $lvtx $lcnx] {}
  
    set lpoly [DisplayInit $w $d $lcnx $lclr]

    after cancel $::run
    set ::run [after 0 [list Display$d $w $lpoly $lvtx $lcnx $lnv $lmv]]
}

# Data structure for models with vertices and connectivity
# -------------------------------------------------------------------
proc ReadData { n } {
  
    set lvtx {}
    set lcnx {}
    set lclr {}
    set txt ""
    
    switch $n {
        tetrahedron {
            set txt "tetrahedron: 4 faces 4 vertices 5 edges" 
            set a [expr {1.0/sqrt(3.0)}]
            set lvtx [list [list $a $a $a] [list $a -$a -$a] \
                                [list -$a $a -$a] [list -$a -$a $a]]
            set lcnx {{0 3 1} {2 0 1} {3 0 2} {1 3 2}}
        }
      
        cube {
            set txt "cube: 6 faces 8 vertices 12 edges"
            set lvtx {{0.7 0.7 0.7} {-0.7 0.7 0.7} {-0.7 -0.7 0.7} {0.7 -0.7 0.7} 
                         {0.7 0.7 -0.7} {-0.7 0.7 -0.7} {-0.7 -0.7 -0.7} {0.7 -0.7 -0.7}}
            set lcnx {{4 7 6 5} {0 1 2 3} {3 2 6 7} {4 5 1 0} {0 3 7 4} {5 6 2 1}}
        }
      
        octahedron {
            set txt "octahedron 8 faces  6 vertices 16 edges"
            set lvtx {{1 0 0} {0 1 0} {-1 0 0} {0 -1 0} {0 0 1} {0 0 -1}}
            set lcnx {{0 1 4} {1 2 4} {2 3 4} {3 0 4}
                         {1 0 5} {2 1 5} {3 2 5} {0 3 5}}
        }
      
        dodecahedron {
            set txt "dodecahedron 12 faces 20 vertices 30 edges"
            set s3 [expr sqrt(3)]
            set s5 [expr sqrt(5)]
            set alpha [expr {sqrt(2.0/(3 + $s5))/$s3}]
            set beta [expr {(1.0 + sqrt(6.0/(3 + $s5) - 2 + 2*sqrt(2.0/(3.0 + $s5))))/$s3}]
            set gamma [expr {1.0/$s3}]
            set lvtx [list \
                [list -$alpha 0 $beta] \
                [list  $alpha 0 $beta] \
                [list -$gamma -$gamma -$gamma] \
                [list -$gamma -$gamma  $gamma] \
                [list -$gamma  $gamma -$gamma] \
                [list -$gamma  $gamma  $gamma] \
                [list  $gamma -$gamma -$gamma] \
                [list  $gamma -$gamma  $gamma] \
                [list  $gamma  $gamma -$gamma] \
                [list  $gamma  $gamma  $gamma] \
                [list  $beta $alpha 0] \
                [list  $beta -$alpha 0] \
                [list -$beta $alpha  0] \
                [list -$beta -$alpha 0] \
                [list -$alpha 0 -$beta] \
                [list  $alpha 0 -$beta] \
                [list  0 $beta $alpha] \
                [list  0 $beta -$alpha] \
                [list  0 -$beta $alpha] \
                [list  0 -$beta -$alpha]]
            set lcnx {{0 1 9 16 5} {1 0 3 18 7} {1 7 11 10 9} {11 7 18 19 6} 
                         {8 17 16 9 10} {2 14 15 6 19} {2 13 12 4 14} {2 19 18 3 13} 
                         {3 0 5 12 13} {6 15 8 10 11} {4 17 8 15 14} {4 12 5 16 17}} 
        } 
      
        icosahedron {
            set txt "icosahedron: 20 faces 12 vertices 30 edges"
            set X 0.525731112119133606
            set Z 0.850650808352039932
            set lvtx [list [list -$X 0.0 $Z] [list $X 0.0 $Z] [list -$X 0.0 -$Z] \
                                [list $X 0.0 -$Z] [list 0.0 $Z $X] [list 0.0 $Z -$X] \
                                [list 0.0 -$Z $X] [list 0.0 -$Z -$X] [list $Z $X 0.0] \
                                [list -$Z $X 0.0] [list $Z -$X 0.0] [list -$Z -$X 0.0]]
            set lcnx {{4 0 1} {9 0 4} {5 9 4} {5 4 8}
                         {8 4 1} {10 8 1} {3 8 10} {3 5 8}
                         {2 5 3} {7 2 3} {10 7 3} {6 7 10}
                         {11 7 6} {0 11 6} {1 0 6} {1 6 10}
                         {0 9 11} {11 9 2} {2 9 5} {2 7 11}}
        }
    }
  
    for {set i 0} {$i <= [llength $lcnx]} {incr i} {
        lappend lclr "0000[format %2.2x 255]"
    }
    return [list $txt $lvtx $lcnx $lclr]
}

# Initialization of canvas with polygonal objects filled or not
# -------------------------------------------------------------------
proc DisplayInit {w d lcnx lclr} {
    set lpoly {}
    set i 0
    if {$d == "Shaded"} then {
        foreach cnx $lcnx {
            lappend lpoly [$w.c create polygon \
                              [string repeat " 0" [expr {2*[llength $cnx]}]] \
                              -fill "#[lindex $lclr $i]"]
            incr i
        }
    } else {
        foreach cnx $lcnx {
            lappend lpoly [$w.c create polygon \
                              [string repeat " 0" [expr {2*[llength $cnx]}]] \
                              -fill black -outline blue]
        }
    }
    return $lpoly
}

# Flat shaded display with gradient color
# -------------------------------------------------------------------
proc DisplayShaded {w lpoly lvtx lcnx lnv lmv} {
    if {$::stop} return

    global iterations
    global ax ay az

    set ax [expr {$ax-0.02}]
    set az [expr {$az+0.02}]
    set ay [expr {$ay+0.025}]
    
    set lgradB {}
    foreach {M lnew lvn} [Transformations $lvtx $lnv] {}
    # Light vector is set to <1 1 -1> 
    foreach i [Intensity $lvn $lmv [list 1 1 -1]] {
        lappend lgradB [format %2.2x [expr {100+154*$i/32}]]
    }
    set i 0
    foreach cnx $lcnx {
        set lcoords {}
        foreach j $cnx {lappend lcoords [lindex $lnew $j]}
        # Backface culing for hidden face. Not removed but only reduced to a point
        if {[lindex $lvn $i 2] < 0} {
            eval $w.c coords [lindex $lpoly $i] [join $lcoords]
            $w.c itemconfigure [lindex $lpoly $i] -fill "#0000[lindex $lgradB $i]"
        } else {
            $w.c coords [lindex $lpoly $i] [string repeat " [join [Barycenter $lcoords]]" [llength $cnx]]
        }
        incr i
    }
    if {[incr ::iterations]} {
        set ::run [after $::speed [list DisplayShaded $w $lpoly $lvtx $lcnx $lnv $lmv]]
    } else {
        return 
    }
}

# Wireframe display 
# -------------------------------------------------------------------
proc DisplayWireframe {w lpoly lvtx lcnx lnv lmv} {
    if {$::stop} return
    global ax az ay

    set ax [expr {$ax-0.02}]
    set az [expr {$az+0.02}]
    set ay [expr {$ay+0.025}]
    foreach {M lnew lvn} [Transformations $lvtx $lnv] {}
    set i 0
    foreach cnx $lcnx {
        set lcoords {}
        foreach j $cnx {lappend lcoords [lindex $lnew $j]}
        # Backface culing for hidden face. Not removed but only reduced to a point
        if {[lindex $lvn $i 2] < 0} {
            eval $w.c coords [lindex $lpoly $i] [join $lcoords]
        } else {
            $w.c coords [lindex $lpoly $i] [string repeat " [join [Barycenter $lcoords]]" [llength $cnx]]
        }
        incr i
    }
    if {[incr ::iterations] } {
        set ::run [after $::speed [list DisplayWireframe $w $lpoly $lvtx $lcnx $lnv $lmv]]
    } else {
        return 
    }
}

# -------------------------------------------------------------------
proc Main {} {
    global stop
    global display
    global scx scy vdist speed
    set ::run {}

    set w .tdc
    catch {destroy $w}
    toplevel $w
    wm withdraw .
    wm title $w "Rotating polyhedra in Tk canvas "

    set display Shaded
    set scx 420
    set scy 420
    set vdist 1200
    set ::scaspeed 40
    set ::speed 40

    pack [canvas $w.c -width $scx -height $scy -bg white -bg black -bd 0]
    $w.c delete all
    bind $w.c <Destroy> {
        after cancel $::run
    }

    set f1 [frame $w.f1 -relief sunken -borderwidth 2]
    pack $f1 -fill x
    button $f1.brun -text Stop -command {set stop 1}
    button $f1.bq   -text Quit -command exit
    label  $f1.l1   -text "  "
    radiobutton $f1.rbs -text "Shaded" -variable display -value Shaded 
    radiobutton $f1.rbw -text "Wireframe" -variable display -value Wireframe
    pack {*}[winfo children $f1] -side left
    set f2 [frame $w.f2 -relief sunken -borderwidth 2]
    pack $f2 -fill x
    foreach i {tetrahedron cube octahedron dodecahedron icosahedron} {
        button $f2.b$i -text $i -command "DisplayModel $w $i"
    }
    pack {*}[winfo children $f2] -side left
    set f3 [frame $w.f3 -relief sunken -borderwidth 2]
    pack $f3 -fill x
    label $f3.l1 -text "View distance " -width 12 
    scale $f3.sca -from 300 -to 1600 -length 300 \
                      -orient horiz -bd 1 -showvalue true -variable vdist 
    pack {*}[winfo children $f3] -side left

    set f4 [frame $w.f4 -relief sunken -borderwidth 2]
    pack $f4 -fill x
    label $f4.l1 -text "Speed " -width 12 
    scale $f4.speed -from 1 -to 99 -length 300 \
                      -orient horiz -bd 1 -showvalue true -variable scaspeed \
                      -command {set speed [expr {100-$scaspeed}];#}
    pack {*}[winfo children $f4] -side left
}

Main