Geometric Calculations of Points and Lines


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.

Software Setup

The nature of this software implementation was as follows:

Constants

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 circumference of the earth at the equator is 40030.13978 km

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

Spherical to Cartesian Coordinates Transformation

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";

..which results in the following result:

16.7832301672891,-59.7079918581292,61.1264468502177

Line Intersection

Below 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 Equations
For 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 Calculation

Below 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



Analytical and Computer Cartography Home