Arjen Markus (11 september 2002) I have had this idea for some time now: create a set of procedures that will allow me to draw mathematical shapes like a torus or even a knot, colour them in ingenious ways, in short, play with 3D space.The script below is just a start. It is not perfect, it is not flexible, but it does give you some sense of what I want.Notes:
- The sorting of the polygons must be improved
- The viewpoint is now hidden in the construction of the torus's polygons - this is a quick hack to avoid the work involved with the first note.
- The viewpoint is not easily changed, well, you can change the angle around the y-axis, but that is all.
. (Thanks to Reinhard Max for providing the space). Note that the full script saves each frame for the animation, so you may want to comment out the "exec" statement at the end.Martin Lemburg - 20.09.2002 - sorry, but does the animation only runs under UNIX or linux? It asks for xwd!AM Actually, you can comment out that line. xwd is used to store the individual frames. (This merits a separate page perhaps). This page will show the animation that resulted. Note that the script below offers more functionality.With this version you can look at three different shapes being wrung or rotated.Enjoy the show!
# solids3d.tcl --
#
# Package for displaying 3D solid bodies
# (sample Workbench module)
#
# Notes:
# This package is a quick hack to get started only
#
# Version information:
# version 0.3: added mapping algorithm for general transformations,
# september 2002
# version 0.2: improved sorting of polygons, september 2002
# (with help from Donal Fellows)
# version 0.1: initial implementation, september 2002
package provide Solids3D 0.1
namespace eval ::Solids3D {
namespace export func deriv
variable display_options
# normalToPlane --
# Return the normal vector to a plane given as three points
#
# Arguments:
# point1 First point
# point2 Second point
# point3 Third point
#
# Result:
# Vector {x,y,z}, a normal vector to the plane
#
proc normalToPlane {point1 point2 point3} {
foreach {dummy px1 py1 pz1} $point1 {break}
foreach {dummy px2 py2 pz2} $point2 {break}
foreach {dummy px3 py3 pz3} $point3 {break}
set vx1 [expr {$px2-$px1}]
set vy1 [expr {$py2-$py1}]
set vz1 [expr {$pz2-$pz1}]
set vx2 [expr {$px3-$px1}]
set vy2 [expr {$py3-$py1}]
set vz2 [expr {$pz3-$pz1}]
set nx [expr {$vy1*$vz2-$vz1*$vy2}]
set ny [expr {$vz1*$vx2-$vx1*$vz2}]
set nz [expr {$vx1*$vy2-$vy1*$vx2}]
return [list VECTOR $nx $ny $nz]
}
# pointOnTorus --
# Return the coordinates of a point on a torus, as given by
# two parameters (angles)
#
# Arguments:
# phi Angle with respect to x-axis
# theta Angle with respect to "inner" axis of torus
#
# Result:
# Point {x,y,z}, the coordinates of the point
#
proc pointOnTorus {phi theta} {
set cosphi [expr {cos($phi)}]
set sinphi [expr {sin($phi)}]
set costh4 [expr {0.25*cos($theta)}]
set sinth4 [expr {0.25*sin($theta)}]
set x [expr {$cosphi*(1.0+$costh4)}]
set y [expr {$sinphi*(1.0+$costh4)}]
set z [expr {1.0+$sinth4}]
return [list POINT $x $y $z]
}
# constructTorus --
# Return a list of polygons, together forming an approximation to a
# torus
#
# Arguments:
# nophi Number of steps along main perimeter
# notheta Number of steps along secondary perimeter
#
# Result:
# List of polygons
#
proc constructTorus {nophi notheta} {
variable angle
set pi 3.1415926
set dphi [expr {2.0*$pi/double($nophi)}]
set dtheta [expr {2.0*$pi/double($notheta)}]
set polygons {POLYGON-LIST}
for { set iphi 0 } { $iphi < $nophi } { incr iphi } {
set phi1 [expr {$dphi*double($iphi)}]
set phi2 [expr {$dphi*double($iphi+1)}]
for { set itheta 0 } { $itheta < $notheta } { incr itheta } {
set theta1 [expr {$dtheta*double($itheta)}]
set theta2 [expr {$dtheta*double($itheta+1)}]
set point1 [pointOnTorus $phi1 $theta1]
set point2 [pointOnTorus $phi1 $theta2]
set point3 [pointOnTorus $phi2 $theta2]
set point4 [pointOnTorus $phi2 $theta1]
set red [expr {int(125*(1.0+cos($phi1)))}]
set green [expr {int(125*(1.0+sin($phi1)))}]
set blue [expr {int(125*(1.0+sin($theta1)))}]
set colour [format "#%2.2X%2.2X%2.2X" $red $green $blue]
set normal [normalToPlane $point1 $point2 $point4]
lappend polygons \
[list POLYGON $colour $normal $point1 $point2 $point3 $point4]
}
}
return $polygons
}
# constructCube --
# Return a list of polygons, together forming a cube
#
# Arguments:
# None
#
# Result:
# List of polygons
#
proc constructCube {} {
set v {VECTOR 1 1 1}
set p1 {POINT -1 -1 -1}
set p2 {POINT 1 -1 -1}
set p3 {POINT 1 1 -1}
set p4 {POINT -1 1 -1}
set p5 {POINT -1 -1 1}
set p6 {POINT 1 -1 1}
set p7 {POINT 1 1 1}
set p8 {POINT -1 1 1}
set f1 [list POLYGON blue $v $p1 $p2 $p3 $p4]
set f2 [list POLYGON green $v $p1 $p2 $p6 $p5]
set f3 [list POLYGON yellow $v $p2 $p3 $p7 $p6]
set f4 [list POLYGON purple $v $p3 $p4 $p8 $p7]
set f5 [list POLYGON orange $v $p1 $p4 $p8 $p5]
set f6 [list POLYGON red $v $p5 $p6 $p7 $p8]
list POLYGON-LIST $f1 $f2 $f3 $f4 $f5 $f6
}
# constructSpiral --
# Return a list of polygons, together forming a spiral
#
# Arguments:
# nosegments Number of segments
#
# Result:
# List of polygons
#
proc constructSpiral {nosegments} {
set result {POLYGON-LIST}
set dz [expr {2.0/double($nosegments)}]
set dphi [expr {4.0*3.1415926/double($nosegments)}]
set r1 0.45
set r2 0.95
set v {VECTOR 1 1 1}
for {set i 0} {$i < $nosegments} {incr i} {
set z1 [expr {-1.0+$dz*$i}]
set z2 [expr {$z1+$dz}]
set phi1 [expr {$dphi*$i}]
set phi2 [expr {$phi1+$dphi}]
set x1 [expr {$r1*cos($phi1)}]
set x2 [expr {$r2*cos($phi1)}]
set y1 [expr {$r1*sin($phi1)}]
set y2 [expr {$r2*sin($phi1)}]
set x3 [expr {$r1*cos($phi2)}]
set x4 [expr {$r2*cos($phi2)}]
set y3 [expr {$r1*sin($phi2)}]
set y4 [expr {$r2*sin($phi2)}]
set p1 [list POINT $x1 $y1 $z1]
set p2 [list POINT $x2 $y2 $z1]
set p3 [list POINT $x3 $y3 $z2]
set p4 [list POINT $x4 $y4 $z2]
set red [expr {int(125*(1.0+cos($phi1)))}]
set green [expr {int(125*(1.0+sin($phi1)))}]
set blue 255
set colour [format "#%2.2X%2.2X%2.2X" $red $green $blue]
lappend result [list POLYGON $colour $v $p1 $p2 $p4 $p3]
}
return $result
}
# rotateX,Y,Z --
# Rotate around the x,y or z axis
#
# Arguments:
# angle Angle over which to rotate
# point Point to rotate
#
# Result:
# New coordinates
#
proc rotateX {angle point} {
foreach {dummy x y z} $point break
set xr $x
set yr [expr {$y*cos($angle)-$z*sin($angle)}]
set zr [expr {$y*sin($angle)+$z*cos($angle)}]
return [list POINT $xr $yr $zr]
}
proc rotateY {angle point} {
foreach {dummy x y z} $point break
set xr [expr {$x*cos($angle)-$z*sin($angle)}]
set yr $y
set zr [expr {$x*sin($angle)+$z*cos($angle)}]
return [list POINT $xr $yr $zr]
}
proc rotateZ {angle point} {
foreach {dummy x y z} $point break
set xr [expr {$x*cos($angle)-$y*sin($angle)}]
set yr [expr {$x*sin($angle)+$y*cos($angle)}]
set zr $z
return [list POINT $xr $yr $zr]
}
# projectOnYZ --
# Project a point on the YZ-plane (and scale as we do this)
#
# Arguments:
# point Point to rotate
#
# Result:
# XY coordinates (for the screen)
#
proc projectOnYZ {point} {
variable scale
variable yoffset
variable zoffset
foreach {dummy x y z} $point break
set xp [expr {int($scale*($y-$yoffset))}]
set yp [expr {int($scale*($zoffset-$z))}]
return [list $xp $yp]
}
# compareView --
# Compare the relative position of two polygons w.r.t. the viewpoint
#
# Arguments:
# polygon1 First polygon
# polygon2 Second polygon
#
# Result:
# 1, -1 if the first polygon lies before or after the second,
#
# Note:
# The comparison is made by looking at the x-coordinate. As the
# projection is in a plane perpendicular to the x-axis, this is
# all we need to do.
#
proc compareView {polygon1 polygon2} {
variable angle
variable viewpoint
#
# Determine the smallest z-coordinate for each point
#
set xmin1 {}
foreach point [lrange $polygon1 3 end] {
foreach {dummy px py pz} $point {break}
if { $xmin1 == {} || $xmin1 > $px } {
set xmin1 $px
}
}
set xmin2 {}
foreach point [lrange $polygon2 3 end] {
foreach {dummy px py pz} $point {break}
if { $xmin2 == {} || $xmin2 > $px } {
set xmin2 $px
}
}
if { $xmin1 < $xmin2 } {
return 1
} else {
return -1
}
}
# applyTorsion --
# Transformation modelling torsion around Y-axis
#
# Arguments:
# point Point to "tord"
#
# Result:
# New point
#
proc applyTorsion {point} {
variable torsion_angle
foreach {type x y z} $point break
set angle [expr {$torsion_angle*$y}]
set xr [expr {$x*cos($angle)-$z*sin($angle)}]
set yr $y
set zr [expr {$x*sin($angle)+$z*cos($angle)}]
return [list $type $xr $yr $zr]
}
# rotateCube --
# Rotate a cube over three axes
#
# Arguments:
# point Point to rotate
#
# Result:
# New point
#
proc rotateCube {point} {
variable angle
set anglex [expr {1.23*$angle}]
set angley [expr {$angle+0.2}]
set anglez [expr {-0.81*$angle+0.44}]
rotateX $anglex [rotateY $angley [rotateZ $anglez $point]]
}
# scaleZ --
# Scale a point in the Z-axis
#
# Arguments:
# point Point to "scale"
#
# Result:
# New point
#
proc scaleZ {point} {
variable scale_factor
foreach {type x y z} $point break
set zr [expr {$z*$scale_factor}]
return [list $type $x $y $zr]
}
# map --
# Map a list of polygons via a transformation of single points and
# vectors
#
# Arguments:
# mapPoint Transformation of a single point or vector
# polygons List of 3D polygons
#
# Result:
# List of transformed polygons
#
proc map {mapPoint polygons} {
#
# Check the type of the second argument (it should be a list whose
# first element gives information about the type):
# - For polygon-lists, call [map] on the individual polygons
# and assemble
# - For polygons, call [mapPoint] on the individual points
# and assemble
#
switch -- [lindex $polygons 0] {
"POLYGON-LIST" {
set result [lindex $polygons 0]
foreach polygon [lrange $polygons 1 end] {
lappend result [map $mapPoint $polygon]
}
}
"POLYGON" {
# $polygons is actually a single polygon now:
# element 0 - the keyword
# element 1 - the colour
# element 2 - the normal vector (should be re-evaluated,
# but right now it is not used)
set result [lrange $polygons 0 2]
foreach point [lrange $polygons 3 end] {
lappend result [$mapPoint $point]
}
}
default {
set result $polygons
}
}
return $result
}
# display --
# Quick and dirty implementation to display a set of polygons
#
# Arguments:
# polygons List of 3D polygons
#
# Result:
# None
#
# Side effect:
# Display of polygons in the canvas
#
proc display {polygons} {
variable angle
variable viewpoint
set viewpoint [list POINT [expr cos($angle)] 0.0 [expr sin($angle)]]
#
# Rotate the polygons first, so that the new coordinates can
# be used directly for sorting and plotting
#
set plane_polygons [lrange $polygons 1 end]
set rotated_polygons {}
foreach polygon $plane_polygons {
set points [lrange $polygon 3 end]
set coords [lrange $polygon 0 2]
foreach point $points {
lappend coords [rotateY $angle $point]
}
lappend rotated_polygons $coords
}
#
# Now sort the polygons w.r.t. the z-coordinate
#
set sorted_polygons [lsort -command compareView $rotated_polygons]
.cnv delete all
foreach polygon $sorted_polygons {
set colour [lindex $polygon 1]
set points [lrange $polygon 3 end]
set coords {}
foreach point $points {
set coords [concat $coords [projectOnYZ $point]]
}
.cnv create polygon $coords -fill $colour -outline black
}
}
#
# Initialise the variables
#
variable angle [expr {0.25*3.1415926}]
variable scale 70.0
variable yoffset -2.0
variable zoffset 2.0
variable torsion_angle 0.0
} ;# End of namespace
#
# Set up a simple UI
#
frame .f
radiobutton .f.radio1 -value torus -variable object -text "Contorted torus" -anchor w
#radiobutton .f.radio2 -value orange -variable object -text "Orange peals" -anchor w
radiobutton .f.radio3 -value cube -variable object -text "Colourful cube" -anchor w
radiobutton .f.radio4 -value spring -variable object -text "Spring" -anchor w
label .f.label -text " " -anchor w
button .f.show -text "Show" -command showAnimation
pack .f.radio1 -side top -fill x
#pack .f.radio2 -side top -fill x
pack .f.radio3 -side top -fill x
pack .f.radio4 -side top -fill x
pack .f.label -side top -fill x
pack .f.show -side top -fill x
pack .f -side left
canvas .cnv -width 300 -height 300 -background white
pack .cnv -fill both
button .exit -text "Exit" -command exit
pack .exit -side bottom
set ::object "torus"
#
# The animation commands showAnimation and next
#
set ::torus [::Solids3D::constructTorus 16 16]
set ::cube [::Solids3D::constructCube]
set ::spring [::Solids3D::constructSpiral 40]
proc showAnimation {} {
global object
global step
set step 0
next $step
}
proc next {step} {
global object
set maxstep 40
switch -- $object {
"torus" {
set ::Solids3D::torsion_angle [expr {1.57*$step/20.0}]
::Solids3D::display [::Solids3D::map applyTorsion $::torus]
}
"cube" {
set maxstep 100
set ::Solids3D::angle [expr {0.1*$step}]
::Solids3D::display [::Solids3D::map rotateCube $::cube]
}
"orange" {
set ::Solids3D::angle [expr {0.1*$step}]
::Solids3D::display [::Solids3D::map rotateCube $::cube]
}
"spring" {
set ::Solids3D::angle [expr {0.03*$step}]
set ::Solids3D::scale_factor [expr {($step+5)/30.0}]
::Solids3D::display \
[::Solids3D::map rotateCube \
[::Solids3D::map scaleZ $::spring]]
}
}
if { $step < $maxstep } {
incr step
after 100 next $step
}
}