## 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]
}```

saito - 2023-07-20

Interesting work! I did a lot of similar calculations as a hobby project many years ago :-)

Do you have something similar for working in basic 2-D coordinate system? Let's assume that points are defined as (x,y), and lines as (p1, p2) where x and y are coordinates, and p1,p2 are points (i.e., p1: (x1,y1), p2: (x2,y2), line1: (p1, p2), etc.). The following functions may be useful to add to this list:

• Find distance between two points
• Find distance form a point to a line
• Find distance between two lines
• Get a line's slope
• Do two lines intersect?
• If two lines intersect, what is the intersection point?
• Given a line, define a new line that is perpendicular to it
• And something hard to explain: given a line as (p1, p2), get two perpendicular lines, one at p1 and the other at p2. Then, given a point, px, how do you figure out whether px lies within the region between the two perpendicular lines?

arjen - 2023-07-21 06:59:54

(Slightly off topic, but still relevant I guess) I have written a small script that allows me to use world coordinates to specify the things on the canvas. It helps me set up sketches. I will post it after my holiday, as it needs to be cleaned up, expanded and perhaps even documented. Here is a sample of what I can do with it:

```setup .c {-10 -10 10 10}
raster
axis -fill green -width 5

set coords {}
for {set x -10} {\$x < 10} {set x [expr {\$x+0.2}]} {
set y [expr {cos(\$x)}]
lappend coords \$x \$y
}
polyline \$coords -width 2

closedcurve [list 5.0 5.0 2.0 -2.0 -4.0 1.0] -outline green -fill {}
closedcurve [list 5.0 5.0 2.0 -2.0 -4.0 1.0] -fill green -smooth 1

circle -2 3 1.5 -fill red

square  6 0 1.5 -fill yellow

showtext -8 8 "This is simply text"```

 Category Mathematics Category Polygon