Scope
This lab exemplifies simple geometric algorithms and formulas, dealing with points, lines and polygons.
All algorithms were prototyped into a software implementation, using the Perl language (http://www.perl.com/), and the various functions will be explained throughout this exercise. Although the implementation level may be language specific, the underlying algorithms are common geometric formulas, which can be applied to any software / programming integrated development environment.
courier
The nature of this software implementation was as follows:
geomath.pm
in this context.mathex1.pl
in this context.We have defined the following constants in our library:
use constant pi => 3.14159; use constant earthRadius => 6371;Generic Functions
The following generic functions were implemented as general purpose routines.
sub deg2rad { # degrees to radians my $degreeNum = $_[0]; my $radianNum = $degreeNum * ($pi / 180); return $radianNum; } sub rad2deg { # radians to degrees my $radianNum = $_[0]; my $degreeNum = $radianNum * (180 / $pi); return $degreeNum; } sub dd2dms { # decimal degrees to degrees minutes seconds my ($D, $d) = split /\./, $_[0]; $d = "." . $d; my $Mm = $d * 60; my @tmp = split /\./, $Mm; my $S = $1 if ($tmp[1] * 60) =~ /^(\d{2})/; return "$D $tmp[0] $S"; } sub dms2dd { # degrees minutes seconds to decimal degrees my ($d, $m, $s) = split / /, $_[0]; my $dd = $d + ($m / 60) + ($s / 3600); return $dd; } sub getMaxVal { # returns greater of two numbers my $val1 = $_[0]; my $val2 = $_[1]; if ($val1 > $val2) { return $val1; } if ($val1 < $val2) { return $val2; } if ($val1 == $val2) { return $val1; } } sub getMinVal { # returns lesser of two numbers my $val1 = $_[0]; my $val2 = $_[1]; if ($val1 < $val2) { return $val1; } if ($val1 > $val2) { return $val2; } if ($val1 == $val2) { return $val1; } }Circumference Calculations
We have implemented the following function to derive the circumference of the earth at a given latitude:
sub getEarthCircumference { my $latitude = $_[0]; my $circumference = 2 * $pi * $earthRadius * cos(deg2rad($latitude)); return $circumference; }
Below is how we call the function from our implementation script:
print " Circumference at 0: " . getEarthCircumference(0) . "\n";
The following code below calculates the circumference of the latitude circles at 15, 35, 49, 75, and 90 degrees:
foreach my $latitude(15, 35, 49, 75, 90) { print " Circumference at N$latitude: " . getEarthCircumference($latitude) . " km\n"; }..which produces the following results:
Circumference at N15: 38666.1481345194 km
Circumference at N35: 32790.7826809389 km
Circumference at N49: 26262.1564579107 km
Circumference at N75: 10360.6053048779 km
Circumference at N90: 0.0531117851733923 km
To convert from spherical (rho, phi, theta) coordinates to cartesian (x, y, z) coordinates
sub spherical2cartesian { my ($rho, $phi, $theta) = @_; my ($x, $y, $z); if ($theta < 0) { $x = $rho * sin($phi) * cos($theta); $y = $rho * sin($phi) * sin($theta); $z = $rho * cos($phi); } if ($theta > 0) { $x = $rho * cos($phi) * cos($theta); $y = $rho * cos($phi) * sin($theta); $z = $rho * sin($phi); } return "$x,$y,$z"; }
To implement this to derive the xyz coordinates of 45 25 00, -75 42 00 we call the function as follows:
print spherical2cartesian((getLineLength(dms2dd("-75 42 00"), dms2dd("45 25 00"), 0, 0)), deg2rad(dms2dd("45 25 00")), deg2rad(dms2dd("-75 42 00"))) . "\n";
16.7832301672891,-59.7079918581292,61.1264468502177
Line IntersectionBelow is the line intersection formula we implemented in our library:
sub lineIntersect { my $x1 = $_[0]; my $y1 = $_[1]; my $x2 = $_[2]; my $y2 = $_[3]; my $x3 = $_[4]; my $y3 = $_[5]; my $x4 = $_[6]; my $y4 = $_[7]; my $xm = getLineSlope($x1, $y1, $x2, $y2); my $ym = getLineSlope($x3, $y3, $x4, $y4); my $xb = getYIntercept($x1, $y1, $xm); my $yb = getYIntercept($x3, $y3, $ym); # check if bounding box of both lines intersect # if they don't, for sure they don't intersect if ( (getMaxVal($x1, $x2) >= getMinVal($x3, $x4)) && (getMaxVal($x3, $x4) >= getMinVal($x1, $x2)) && (getMaxVal($y1, $y2) >= getMinVal($y3, $y4)) && (getMaxVal($y3, $y4) >= getMinVal($y1, $y2))) { # print "Bounding boxes of both lines intersect\n"; } else { print " Lines do not physically intersect"; } my $xInt = ($xb - $yb) / ($ym - $xm); my $yInt = (((-$xm) * $yb) + ($ym * $xb)) / ($ym - $xm); return "$xInt,$yInt"; }
This function uses the following two additional functions to derive slope and y intercept values:
sub getLineSlope { my $x1 = $_[0]; my $y1 = $_[1]; my $x2 = $_[2]; my $y2 = $_[3]; my $slope = ($y2 - $y1) / ($x2 - $x1); return $slope; } sub getYIntercept { my $x = $_[0]; my $y = $_[1]; my $m = $_[2]; my $yIntercept = $y - $m * $x; return $yIntercept; }We use the lineIntersect function to derive the intersection point by passing the endpoints of two lines:
print " Intersection point: " . lineIntersect(1, 3, 10, 7, 2, 7, 8, 2) . "\n"; print " Intersection point: " . lineIntersect(1, 1, 4, 5, 7, 6, 10, 3) . "\n"; The intersection point of the two lines given by endpoints: Line 1: A(1, 3) B(10, 7) Line 2: C(2, 7) D(8, 2) is: 4.78260869565217, 4.68115942028986 (lines physically intersect) The intersection point of the two lines given by endpoints: Line 1: A(1, 1) B(4, 5) Line 2: C(7, 6) D(10, 3) is: 5.71428571428572, 7.28571428571429 (lines DO NOT physically intersect)The
getLineMidPoint
function
This function derives the xy midpoint of a given line defined by two endpoints. We can derive a test to find out whether the intersection xy points above are indeed the midpoints of the lines they were calculated on. The getLineMidPoint function is shown below:
sub getLineMidpoint { my $x1 = $_[0]; my $y1 = $_[1]; my $x2 = $_[2]; my $y2 = $_[3]; my $newX = ($x1 + $x2) / 2; my $newY = ($y1 + $y2) / 2; return "$newX,$newY"; }
Some pseudocode below describes this process
XYInt1 = Get x y intersection point of line AB XYInt2 = Get x y intersection point of line CD XYInt3 = Get midpoint of line AB XYInt4 = Get midpoint of line CD // Compare XYInt1 to XYint3 If XYint1 == XYint3 Intersection is at midpoints of line AB Else Intersection is not at midpoints of line AB // Compare XYInt2 to XYint4 If XYint2 == XYint4 Intersection is at midpoints of line CD Else Intersection is not at midpoints of line CD
The test shows that neither calculated intersection point in questions 4 and 5 lies between the endpoints of both lines.
Computing Intersection Points of Line EquationsFor the following two lines: y = 3x/2 + 7 x = 2y/3 - 5 ..normalize Y in equation 2 y = 3/2x + 15 / 2 3/2x + 7 = 3/2x + 15/2 7 = 15 / 2 FALSE
Therefore, the lines do not intersect, as they have the same slope
Line Length CalculationBelow is the code implemented to compute the length of a line defined by two endpoints:
sub getLineLength { my $x1 = $_[0]; my $y1 = $_[1]; my $x2 = $_[2]; my $y2 = $_[3]; my $theLength = sqrt((($x2 - $x1) * ($x2 - $x1)) + (($y2 - $y1) * ($y2 - $y1))); return $theLength; }
Using this as our base function, we call this function from our implementation script, with the specific coordinate details:
# define x coordinates array my @xList = (1, 6, 2, 3, 5, 2, 8); # define y coordinates array my @yList = (2, 9, 4, 7, 2, 9, 6); # loop through each coord set to derive line length # add to $total variable my $arrLength = scalar(@xList); for(my $i = 0; $i < ($arrLength - 1); $i++) { my $nextEl = $i + 1; $total += getLineLength($xList[$i], $yList[$i], $xList[$nextEl], $yList[$nextEl]); }
As a result, the length of the following line of coordinates
1. 1,2 2. 2,4 3. 5,2 4. 8,6 5. 6,9 6. 3,7 7. 2,9 is: 20.2887897813915
The entire implementation program is available for download from:
mathex1.pl
The library is available for download from:
geomath.pm