## Canvas geometrical calculations

TR - Geometrical calculations are important to many applications like GIS or drawing programs. Perhaps we can collect some good algorithms here. I had a look at http://www.geometryalgorithms.com , which is a great source of all kinds of calculations we might need for canvas items.

The area of a polygon

``` # compute the area of a polygon given its coordinates
#
# Argument: coords -> list of coordinates
#
# Returns: Area of polygon (taking care of holes inside the polygon)
#
proc polygonArea {coords}  {
# make sure we have a closed set of coords:
if {[lindex \$coords 0] != [lindex \$coords end-1] \
|| [lindex \$coords 1] != [lindex \$coords end]} {
lappend coords [lindex \$coords 0] [lindex \$coords 1]
}
# append another point for the calculation:
lappend coords [lindex \$coords 2] [lindex \$coords 3]
# area:
set area 0
# number of vertices:
set n [expr {([llength \$coords]-4)/2}]
# build lists with x and y coordinates only:
foreach {x y} \$coords {
lappend xList \$x
lappend yList \$y
}
for {set i 1; set j 2; set k 0} {\$i <= \$n} {incr i; incr j; incr k} {
set area [expr {\$area + ([lindex \$xList \$i] * ([lindex \$yList \$j] - [lindex \$yList \$k]))}]
}
return [expr {\$area/2.0}]
}```

This procedure may seem slooow, but trying it I got the area of a polygon with 2100 vertices in 0.0037 seconds on a 800MHz computer.

The algorithms takes care of holes. Try this polygon (.c is a canvas widget):

`  .c create polygon 0 0 100 0 100 100 50 100 50 75 75 75 75 25 25 25 25 75 50 75 50 100 0 100`

It has a hole inside (how sad, that independant holes are not suported by Tcl, we can only fake them here). The area is 7500. Without the hole

`  .c create polygon 0 0 100 0 100 100 0 100`

the procedure returns 10000. So the hole is treated correctly as not belonging to the inside of the polygon.

The distance between two points

This is trivial when you have a cartesian coordinate system. You just use Pythagoras' law. But in applications with geographical data, the distance must be calculated on a sphere in the simplest case. This is what the following procedure does:

``` # calculates the distance between two points on a sphere.
# The result is the distance on the great circle in Meters
# x1=longitude 1, y1=latitide 1, x2=longitude2, y2=latitude 2 - all given in decimal degrees
proc distance {x1 y1 x2 y2} {
set pi 3.1415926535
set x1 [expr {\$x1 *2*\$pi/360.0}]
set x2 [expr {\$x2 *2*\$pi/360.0}]
set y1 [expr {\$y1 *2*\$pi/360.0}]
set y2 [expr {\$y2 *2*\$pi/360.0}]
# calculate distance:
set d [expr {acos(sin(\$y1)*sin(\$y2)+cos(\$y1)*cos(\$y2)*cos(\$x1-\$x2))}]
# return distance in meters:
return [expr {20001600/\$pi*\$d}]
}```

Normally, d will be the distance measured in radians. The last expr is needed to translate the value to meters. It is the same as:

`    return [expr {180*60.0/\$pi*\$d*1852}]`

Here, 1852 is the number of meters that go into one nautical mile. By using other factors, you can get the result in other units.

This distance will be ok for general applications that don't have special requirements, but if you want to be more exact you have to explore the field of geodesy.

Intersection of two lines, Find Angle, Bisect Angle, and Trisect Angle

KPV 2003-07-10] -- In writing Triangle Madness I had to come up with a bunch of geometric functions that I thought I'd extract and put here. These routines all work with (x,y) duples and there are a few helper routines to manipulate duples.

NB. there's also an Intersect routine (called crosspoint) in Sun, moon, and stars.

``` #+##########################################################################
#
# Intersect -- find where two line intersect given two points on each line
# IntersectV -- same but w/ 2 point/vector pairs
#   we want K,J where p1+K(v1) = p3+J(v3)
#   convert into and solve matrix equation (a b / c d) ( K / J) = ( e / f )
#
proc Intersect {p1 p2 p3 p4} {
return [IntersectV \$p1 [VSub \$p2 \$p1] \$p3 [VSub \$p4 \$p3]]
}
proc IntersectV {p1 v1 p3 v3} {
foreach {x1 y1} \$p1 {vx1 vy1} \$v1 {x3 y3} \$p3 {vx3 vy3} \$v3 break

set a \$vx1
set b [expr {-1 * \$vx3}]
set c \$vy1
set d [expr {-1 * \$vy3}]
set e [expr {\$x3 - \$x1}]
set f [expr {\$y3 - \$y1}]

set det [expr {double(\$a*\$d - \$b*\$c)}]
if {\$det == 0} {error "Determinant is 0"}

set k [expr {(\$d*\$e - \$b*\$f) / \$det}]
#set j [expr {(\$a*\$f - \$c*\$e) / \$det}]
}
##+##########################################################################
#
# FindAngle -- returns the angle formed by p1,p2,p3
# Computes the dot product on vectors p1-p2, p3-p2
#
proc FindAngle {p1 p2 p3} {
foreach {x1 y1} [VSub \$p1 \$p2] {x2 y2} [VSub \$p3 \$p2] break

set m1 [expr {hypot(\$x1,\$y1)}]
set m2 [expr {hypot(\$x2,\$y2)}]
if {\$m1 == 0 || \$m2 == 0} { return 0 }      ;# Coincidental points
set dot [expr {\$x1 * \$x2 + \$y1 * \$y2}]

set theta [expr {acos(\$dot / \$m1 / \$m2)}]
if {\$theta < 1e-5} {set theta 0}
return \$theta
}
##+##########################################################################
#
# BisectAngle -- returns point on bisector of the angle formed by p1,p2,p3
#
proc BisectAngle {p1 p2 p3} {
foreach {x1 y1} [VSub \$p1 \$p2] {x2 y2} [VSub \$p3 \$p2] break

set s1 [expr {100.0 / hypot(\$x1, \$y1)}]
set s2 [expr {100.0 / hypot(\$x2, \$y2)}]
set v1 [VAdd \$p2 [list \$x1 \$y1] \$s1]        ;# 100*Unit vector from p2 to p1
set v2 [VAdd \$p2 [list \$x2 \$y2] \$s2]        ;# 100*Unit vector from p2 to p3
return [VAdd \$v1 [VSub \$v2 \$v1] .5]
}
##+##########################################################################
#
# TrisectAngle -- returns two points which are on the two lines trisecting
# the angle created by points p1,p2,p3. We use the cross product to tell
# us clockwise ordering.
#
proc TrisectAngle {p1 p2 p3} {
set cross [VCross [VSub \$p2 \$p1] [VSub \$p2 \$p3]]
if {\$cross < 0} {foreach {p1 p3} [list \$p3 \$p1] break}

set theta [FindAngle \$p1 \$p2 \$p3]           ;# What the angle is
set theta1 [expr {\$theta / 3.0}]            ;# 1/3 of that angle
set theta2 [expr {2 * \$theta1}]             ;# 2/3 of that angle

set v [VSub \$p3 \$p2]                        ;# We'll rotate this leg
set v1 [VRotate \$v \$theta1]                 ;# By 1/3
set v2 [VRotate \$v \$theta2]                 ;# By 2/3
set t1 [VAdd \$p2 \$v1]                       ;# Trisect point 1
set t2 [VAdd \$p2 \$v2]                       ;# Trisect point 2

if {\$cross < 0} {return [list \$t2 \$t1]}
return [list \$t1 \$t2]
}
##+##########################################################################
#
# HELPER FUNCTIONS:
#   VSub -- subtract two vectors
#   VCross -- returns the cross product for 2D vectors (z=0)
#   VRotate -- rotates vector counter-clockwise
#
proc VAdd {v1 v2 {scaling 1}} {
foreach {x1 y1} \$v1 {x2 y2} \$v2 break
return [list [expr {\$x1 + \$scaling*\$x2}] [expr {\$y1 + \$scaling*\$y2}]]
}
proc VSub {v1 v2} { return [VAdd \$v1 \$v2 -1] }
proc VCross {v1 v2} {
foreach {x1 y1} \$v1 {x2 y2} \$v2 break
return [expr {(\$x1*\$y2) - (\$y1*\$x2)}]
}
proc VRotate {v beta} {
foreach {x y} \$v break
set xx [expr {\$x * cos(-\$beta) - \$y * sin(-\$beta)}]
set yy [expr {\$x * sin(-\$beta) + \$y * cos(-\$beta)}]
return [list \$xx \$yy]
}```