package geomath;

# geomath.pm
# geometric function library
# Author : Tom Kralidis
# email  : tomkralidis@hotmail.com
# website: http://www.kralidis.ca/
#

require Exporter;
use strict;

use POSIX;

use vars       qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);

# set the version for version checking
$VERSION     = 1.0;

@ISA         = qw(Exporter);

@EXPORT      = qw(
deg2rad
rad2deg
scaleXY
translateXY
rotateXY
getLineLength
lineIntersect
getLineSlope
getYIntercept
getMaxVal
getMinVal
getEarthCircumference
getLineMidpoint
getEuclideanDistance
spherical2cartesian
cartesian2spherical
dms2dd
dd2dms
getPointToLineDistance
getPolygonAreaTrapezoidal
getPolygonAreaTriangular
getPolygonPerimeter
getPolygonCentroidTriangular
getPolygonCentroidByCoordAvg
getMBR
getPolygonFromMBR
isEvenXYPairs
isPolygonClosed
isPointInPolygon
);

my $pi = 3.14159;
my $earthRadius = 6371;

# functions

sub deg2rad {
  my $degreeNum = $_[0];
  my $radianNum = $degreeNum * ($pi / 180);
  return $radianNum;
}

sub rad2deg {
  my $radianNum = $_[0];
  my $degreeNum = $radianNum * (180 / $pi);
  return $degreeNum;
}

sub scaleXY {
  my $theX        = $_[0];
  my $theY        = $_[1];
  my $scaleFactor = $_[2];

  my $newX = $theX * $scaleFactor;
  my $newY = $theY * $scaleFactor;
  return "$newX,$newY";
}

sub translateXY {
  my $theX   = $_[0];
  my $theY   = $_[1];
  my $shiftX = $_[2];
  my $shiftY = $_[3];

  my $newX = $theX + $shiftX;
  my $newY = $theY + $shiftY;
  return "$newX,$newY";
}

sub rotateXY {
  my $theX            = $_[0];
  my $theY            = $_[1];
  my $rotationTheta   = $_[2];
  my $rotationOriginX = $_[3];
  my $rotationOriginY = $_[4];

  my $newX;
  my $newY;

  if ($rotationOriginX == 0 && $rotationOriginY == 0) {
    $newX = $theX * cos(deg2rad($rotationTheta)) - sin(deg2rad($rotationTheta)) * $theY;
    $newY = $theX * sin(deg2rad($rotationTheta)) + cos(deg2rad($rotationTheta)) * $theY;
  }
  else {
    $newX = ($theX - $rotationOriginX) * cos(deg2rad($rotationTheta)) - ($theY - $rotationOriginY) * sin($rotationTheta) + $rotationOriginX;
    $newY = ($theX - $rotationOriginX) * sin(deg2rad($rotationTheta)) + ($theY - $rotationOriginY) * cos($rotationTheta) + $rotationOriginY;
  }
  return "$newX,$newY";
}

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

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 {
  #  return " Lines do not intersect";
    print " Lines do not physically intersect";
  }

  my $xInt = ($xb - $yb) / ($ym - $xm);
  my $yInt = ((($xm * -1) * $yb) + ($ym * $xb)) / ($ym - $xm);
  return "$xInt,$yInt";
}

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

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

sub getMaxVal {
  my $val1 = $_[0];
  my $val2 = $_[1];

  if ($val1 > $val2) {
    return $val1;
  }
  if ($val1 < $val2) {
    return $val2;
  }
  if ($val1 == $val2) {
   return $val1;
  }
}

sub getMinVal {
  my $val1 = $_[0];
  my $val2 = $_[1];

  if ($val1 < $val2) {
    return $val1;
  }
  if ($val1 > $val2) {
    return $val2;
  }
  if ($val1 == $val2) {
    return $val1;
  }
}

sub getEarthCircumference {
  my $latitude = $_[0];
  my $circumference = 2 * $pi * $earthRadius * cos(deg2rad($latitude));
  return $circumference;
}

sub dd2dms {
  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 {
  my ($d, $m, $s) = split / /, $_[0];
  my $dd = $d + ($m / 60) + ($s / 3600);
  return $dd;
}

sub getEuclideanDistance{
  my $x1 = $_[0];
  my $y1 = $_[1];
  my $x2 = $_[2];
  my $y2 = $_[3];

  my $dist = sqrt((($x2 - $x1) ** 2) + (($y2 - $y1) ** 2));
  return $dist;
}

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

sub cartesian2spherical {
  my ($x, $y, $z) = @_;
  my $rho   = sqrt (($x*$x) + ($y*$y) + ($z*$z));
  my $theta = atan($y/$x);
  my $phi   = acos($z/sqrt (($x**2) + ($y**2) + ($z**2)));
  return "$rho,$phi,$theta";
}

sub getPointToLineDistance {
  my $x1 = $_[0];
  my $y1 = $_[1];
  my $x2 = $_[2];
  my $y2 = $_[3];
  my $x3 = $_[4]; # x of lone point
  my $y3 = $_[5]; # y of lone point

  my $A = $x3 - $x1;
  my $B = $y3 - $y1;
  my $m = getLineSlope($x1, $y1, $x2, $y2); 

  my $dist = ($B - $m * $A) / sqrt(($m ** 2) + 1); 
  return $dist;
}

sub getPolygonPerimeter {
  my @polygonVertices = @_;
  my $total           = 0;

  if (isPolygonClosed(@polygonVertices) == 1)  {
    return "Corrupt inputs, not a closed polygon";
  }

  if (isEvenXYPairs(@polygonVertices) == 1) {
    return "Corrupt inputs, must be even numbered xy vertices";
  }

  for (my $i = 0; $i < (scalar(@polygonVertices) - 2); $i = $i + 2) {
    $total += sqrt(
     (($polygonVertices[$i] - $polygonVertices[$i + 2]) ** 2) + 
     (($polygonVertices[$i + 1] - $polygonVertices[$i + 3]) ** 2));  
  }
  return $total;
}

sub isPolygonClosed {
  my @polygonVertices = @_;
  if ( ($polygonVertices[0] != $polygonVertices[-2]) ||
       ($polygonVertices[1] != $polygonVertices[-1])) {
    return 1;
  }
  else {
    return 0;
  }
}

sub isEvenXYPairs {
  my @coordArray = @_;
  if (! (scalar(@coordArray) % 2) == 0) {
    return 1;
  }
  else {
    return 0;
  }
}

sub getPolygonAreaTrapezoidal {
  my $total = 0;
  my @polygonVertices = @_;

  if (isEvenXYPairs(@polygonVertices) == 1) {
    return "Corrupt inputs, must be even numbered xy vertices";
  }

  if (isPolygonClosed(@polygonVertices) == 1)  {
    return "Corrupt inputs, not a closed polygon";
  }

  for(my $i = 0; $i < (scalar(@polygonVertices) - 2); $i = $i + 2) {
    $total += ($polygonVertices[$i+2] - $polygonVertices[$i]) * ($polygonVertices[$i+1] + $polygonVertices[$i+3]);
  }
  $total *= 0.5;
  return $total;
}

sub getPolygonAreaTriangular {
  my $total = 0;
  my @polygonVertices = @_;

  if (isEvenXYPairs(@polygonVertices) == 1) {
    return "Corrupt inputs, must be even numbered xy vertices";
  }

  if (isPolygonClosed(@polygonVertices) == 1)  {
    return "Corrupt inputs, not a closed polygon";
  }

  for(my $i = 0; $i < (scalar(@polygonVertices) - 2); $i = $i + 2) {
    $total += ( ($polygonVertices[$i] * $polygonVertices[$i + 3]) -
                ($polygonVertices[$i + 2] * $polygonVertices[$i + 1]));
  }
  $total *= -0.5;
  return $total;
}

sub getPolygonCentroidTriangular {
  my @polygonVertices = @_;

  my $totalX = 0;
  my $totalY = 0;
  my $centX  = 0;
  my $centY  = 0;

  my $polygonArea = 0;

  if (isEvenXYPairs(@polygonVertices) == 1) {
    return "Corrupt inputs, must be even numbered xy vertices";
  }

  if (isPolygonClosed(@polygonVertices) == 1)  {
    return "Corrupt inputs, not a closed polygon";
  }

  $polygonArea = getPolygonAreaTriangular(@_);

  my $arrLength = (scalar(@polygonVertices) / 2) - 1;

  for (my $i = 0; $i < ($arrLength + 3); $i = $i + 2) {
    $totalX += ($polygonVertices[$i] + $polygonVertices[$i+2]) * ($polygonVertices[$i+2] * $polygonVertices[$i+1] - $polygonVertices[$i] * $polygonVertices[$i+3]);

    $totalY += ($polygonVertices[$i+1] + $polygonVertices[$i+3]) * ($polygonVertices[$i+2] * $polygonVertices[$i+1] - $polygonVertices[$i] * $polygonVertices[$i+3]);
  }

  $centX = $totalX / (6 * $polygonArea);
  $centY = $totalY / (6 * $polygonArea);
  return "$centX,$centY";
}

sub getPolygonCentroidByCoordAvg {
  my @polygonVertices = @_;

  my $totalX = 0;
  my $totalY = 0;
  my $centX  = 0;
  my $centY  = 0;

  if (isEvenXYPairs(@polygonVertices) == 1) {
    return "Corrupt inputs, must be even numbered xy vertices";
  }

  if (isPolygonClosed(@polygonVertices) == 1)  {
    return "Corrupt inputs, not a closed polygon";
  }

  my $arrLength = (scalar(@polygonVertices) / 2) - 1;

  for (my $i = 0; $i < ($arrLength + 3); $i = $i + 2) {
    $totalX += $polygonVertices[$i]; 
    $totalY += $polygonVertices[$i + 1]; 
  }

  $centX = $totalX / $arrLength;
  $centY = $totalY / $arrLength;

  return "$centX,$centY";
}


sub isPointInPolygon {
  my ( $x, $y, @xy ) = @_;

  my $n = @xy / 2;                      # Number of points in polygon.
  my @i = map { 2 * $_ } 0 .. (@xy/2);  # The even indices of @xy.
  my @x = map { $xy[ $_ ]     } @i;     # Even indices: x-coordinates.
  my @y = map { $xy[ $_ + 1 ] } @i;     # Odd indices: y-coordinates.

  my ( $i, $j );                        # Indices.

  my $side = 0;                         # 0 = outside, 1 = inside.

  for ( $i = 0, $j = $n - 1 ; $i < $n; $j = $i++ ) {
    if (
         (
           # If the y is between the (y-) borders ...
           ( ( $y[ $i ] <= $y ) && ( $y < $y[ $j ] ) ) ||
           ( ( $y[ $j ] <= $y ) && ( $y < $y[ $i ] ) )
         )
         and
         # ...the (x,y) to infinity line crosses the edge
         # from the ith point to the jth point...
         ($x
         <
         ( $x[ $j ] - $x[ $i ] ) *
         ( $y - $y[ $i ] ) / ( $y[ $j ] - $y[ $i ] ) + $x[ $i ] ))
     {
       $side = not $side; # Jump the fence.
     }
  }
  return $side ? 1 : 0;
}

sub getMBR {
  my @polygonVertices = @_;
  my $minx = 0;
  my $miny = 0;
  my $maxx = 0;
  my $maxy = 0;

  if (isEvenXYPairs(@polygonVertices) == 1) {
    return "Corrupt inputs, must be even numbered xy vertices";
  }

  if (isPolygonClosed(@polygonVertices) == 1)  {
    return "Corrupt inputs, not a closed polygon";
  }

  for(my $i = 0; $i < (scalar(@polygonVertices) / 2 + 2); $i = $i + 2) {
    # initialize bounding box vals
    $minx = $polygonVertices[$i] if $i == 0;
    $miny = $polygonVertices[$i+1] if $i == 0;
    $maxx = $polygonVertices[$i] if $i == 0;
    $maxy = $polygonVertices[$i+1] if $i == 0;

    # test and set MBR
    $maxx = $polygonVertices[$i] if $polygonVertices[$i] > $maxx;
    $minx = $polygonVertices[$i] if $polygonVertices[$i] < $minx;
    $maxy = $polygonVertices[$i+1] if $polygonVertices[$i+1] > $maxy;
    $miny = $polygonVertices[$i+1] if $polygonVertices[$i+1] < $miny;
  }
  return "$minx,$miny,$maxx,$maxy";
}

sub getPolygonFromMBR {
  my ($minx, $miny, $maxx, $maxy) = split /,/, $_[0];
  my @polyOut = ($minx, $miny, $minx, $maxy, $maxx, $maxy, $maxx, $miny, $minx, $miny);
  my $str  = "$minx,$miny,$minx,$maxy,$maxx,$maxy,$maxx,$miny,$minx,$miny";
  return $str;
}

1;