[GS] - This code displays 3d polyhedra in shaded or wireframe mode. 
It uses only the tk canvas. The illumination model is a simple flat shading model 
[http://online.caup.washington.edu/courses/arch481/00HomePage/1.Tapestry%20Reader/2.Rendering/4.Shading/2.Flat%20Shading.html].
The color intensity of a face is proportional to the angle between its normal and 
a light direction.

[http://gersoo.free.fr/wiki/w14283/polyhedra.jpg] [http://gersoo.free.fr/wiki/w14283/polyhedra1.jpg]

- A starkit version with more demos is available at [http://gersoo.free.fr/wiki/w14283/3display.kit]

- A lightweight tclet version can be seen at [http://gersoo.free.fr/wiki/w14283/3display-tclet.html] 
(sources [http://gersoo.free.fr/wiki/w14283/3display-tclet.txt])

The hidden face removal algorithm works well with convex objects but is very limited for 
the others. See for instance the torus or the shuttle as bad examples.

[http://gersoo.free.fr/wiki/w14283/3display.jpg] [http://gersoo.free.fr/wiki/w14283/3display1.jpg]

[PWQ] '''2 May 05''', moved bind command to after package require Tk!

 # polyhedra.tcl 
 # Author:      Gerard Sookahet
 # Date:        30 Mai 2005 
 # Description: Rotating polyhedra using a 'standard' tk canvas.
 #              Flat shading and wireframe mode.
 
 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 
 
  update 
  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
 
  $w.c delete all
  set stop 0
  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]
   
  if {$d == "shaded"} then {
   for {set i 1} {$i<=820} {incr i} {
      if $stop break
      set ax [expr {$ax-0.02}]
      set az [expr {$az+0.02}]
      set ay [expr {$ay+0.025}]
      after 40
      DisplayShaded $w $lpoly $lvtx $lcnx $lnv $lmv
   }
  } else {
   for {set i 1} {$i<=820} {incr i} {
      if $stop break
      set ax [expr {$ax-0.02}]
      set az [expr {$az+0.02}]
      set ay [expr {$ay+0.025}]
      after 40
      Display $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} {
  update 
  
  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
  }
 }
 
 # Wireframe display 
 # -------------------------------------------------------------------
 proc Display {w lpoly lvtx lcnx lnv lmv} {
  update 
  
  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
  }
 }
 
 # -------------------------------------------------------------------
 proc Main {} {
     global stop
     global display
     global scx scy vdist
 
  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
 
  pack [canvas $w.c -width $scx -height $scy -bg white -bg black -bd 0]
  $w.c delete all
 
  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
  eval 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"
  }
  eval 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 "
  scale $f3.sca -from 300 -to 1600 -length 300 \
                -orient horiz -bd 1 -showvalue true -variable vdist 
  eval pack [winfo children $f3] -side left
 }
 
 Main

----
[Peter Newman] 2 June 2005: This is without a shadow of doubt, the best piece of software 
I have ever seen/used (Charles Moore's Forth interpreter being its only possible competition). 
Quake and Pov-Ray are all very powerful. But so complicated. And don't even talk about the rest. 
But you've got useful 3D, in a script language, on a Tk canvas, in a few A4 pages of code. Unbelievable!

But please make all my Christmases come true at once. How does one create the models? 
Please tell me that there's an equally easy to use and powerful modelling tool as well.

[GS] Thanks for all the kudos :-) The models were converted from .obj files since the data 
are almost similar (vertex coordinates and conectivity).

----
[AM] (2 june 2005) I have a few comments on this delightful starkit:

   * Pressing the little x to close the application causes it to display messages about .tdc.c not being a valid command rather than closing the program.
   * The hidden-surface algorithm is not perfect: it fails for instance with the torus.

[GS] Unfortunately, the hidden surface algorithm works only with convex objects. To go further, 
we must use a Z-buffer algorithm, which is more time-consuming since it works at the pixel level.

[AM] In [Three-dimensional shapes] I used a simple method: estimate the distance to the viewpoint and sort the polygons with respect to that distance. It seems to work well as long as the polygons are indeed flat. Since some transformations I apply there are non-linear that is at some point no longer the case. 
It may also be slower than the method applied here (I must admit I have not studied it yet). 

[Lars H]: I was just going to suggest the same thing...  It wouldn't be foolproof, because if the polygons vary in size you could have e.g. a large polygon coming close to the viewpoint whose most distant end is obscured by a smaller polygon in between:
   \
    \  A
     \
      \
 E     \
        \
         \
       B  \
           \
E is the viewpoint, the diagonal line (D) is a large polygon. It's hard to tell from 
distance calculations alone that D should obscure A but not B.
However in real life it would probably work quite well (and fast).

[GS] (050604) Thanks for these suggestions. I have tested a version with polygon depth sorting. 
It is available as starkit here [http://gersoo.free.fr/wiki/w14283/3display_lsort.kit]. 
It sorts a list mean distances. The mean distance of a polygon is the distance between 
the viewpoint and its barycenter. The torus example always failled.

[PWQ] ''23 Aug 05'', Elron the Elf says ''Insert the mystical code below and all will '''NOT''' be revealed''.

 set p 0
 foreach j [lsort_index $lbcz] {
      ...
      $w.c coords [lindex $lpoly $p] [join $lcoords]
      $w.c itemconfigure [lindex $lpoly $p] -fill "#0000[lindex $lgradB $j]"
      ...
    incr p
 }
----
See also [Advanced Optimisation - TestCase 3d polyhedra]
----
[[
[Category Toys] | 
[Category Graphics]
[Category 3D Graphics]
]]