package geoalglib;

# NOTE: some algorithms are taken from the following source:
#
# Authors   : Jon Orwant, Jarkko Hietaniemi, John Macdonald
# Title     : Mastering Algorithms with Perl
# Publisher : O'Reilly and Associates
# ISBN      : 1-56592-398-7
# URL       : http://www.oreilly.com/catalog/maperl/
#
# also, read O'Reilly's code policy at:
# http://www.oreilly.com/ask_tim/codepolicy_1101.html
#
require Exporter;
use vars       qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);

use constant epsilon => 1e-14;

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

@ISA         = qw(Exporter);
@EXPORT      = qw(
bounding_box_of_points
clockwise
closest_points
euclidean
great_circle_dist
great_circle_dist_err
line_intersection
manhattan
point_in_polygon
point_in_quadrangle
point_in_triangle
polygon_area
polygon_perimeter
triangle_area_heron
triangle_perimeter
spatial_extent

FileSizeSingleBandRaster
FileSizeSingleBandRasterHeader
FileSizeMultiBandRaster
FileSizeMultiBandRasterHeader
FileSizeMultiBandRasterMultiHeader

);

%EXPORT_TAGS = ( );     # eg: TAG => [ qw!name1 name2! ],

# your exported package globals go here,
# as well as any optionally exported functions
@EXPORT_OK   = qw($Var1 %Hashit &func3);

use vars qw($Var1 %Hashit);
# non-exported package globals go here
use vars      qw(@more $stuff);

# initialize package globals, first exported ones
$Var1   = '';
%Hashit = ();

# then the others (which are still accessible as $Some::Module::stuff)
$stuff  = '';
@more   = ();

# all file-scoped lexicals must be created before
# the functions below that use them.

# file-private lexicals go here
my $priv_var    = '';
my %secret_hash = ();

# here's a file-private function as a closure,
# callable as &$priv_func.
my $priv_func = sub {
    # stuff goes here.
};

# make all your functions, whether exported or not;
# remember to put something interesting in the {} stubs

# bounding_box($d, @p [,@b])
#   Return the bounding box of the points @p in $d dimensions.
#   The @b is an optional initial bounding box: we can use this
#   to create a cumulative bounding box that includes boxes found
#   by earlier runs of the subroutine (this feature is used by
#   bounding_box_of_points()).
#
#   The bounding box is returned as a list.  The first $d elements
#   are the minimum coordinates, the last $d elements are the
#   maximum coordinates.

sub bounding_box {
    my ( $d, @bb ) = @_; # $d is the number of dimensions.
    # Extract the points, leave the bounding box.
    my @p = splice( @bb, 0, @bb - 2 * $d );

    @bb = ( @p, @p ) unless @bb;

    # Scan each coordinate and remember the extrema.
    for ( my $i = 0; $i < $d; $i++ ) {
        for ( my $j = 0; $j < @p; $j += $d ) {
            my $ij = $i + $j;
            # The minima.
            $bb[ $i      ] = $p[ $ij ] if $p[ $ij ] < $bb[ $i      ];
            # The maxima.
            $bb[ $i + $d ] = $p[ $ij ] if $p[ $ij ] > $bb[ $i + $d ];
        }
    }

    return @bb;
}

# bounding_box_intersect($d, @a, @b)
#   Return true if the given bounding boxes @a and @b intersect
#   in $d dimensions.  Used by line_intersection().

sub bounding_box_intersect {
    my ( $d, @bb ) = @_; # Number of dimensions and box coordinates.
    my @aa = splice( @bb, 0, 2 * $d ); # The first box.
    # (@bb is the second one.)

    # Must intersect in all dimensions.
    for ( my $i_min = 0; $i_min < $d; $i_min++ ) {
        my $i_max = $i_min + $d; # The index for the maximum.
        return 0 if ( $aa[ $i_max ] + epsilon ) < $bb[ $i_min ];
        return 0 if ( $bb[ $i_max ] + epsilon ) < $aa[ $i_min ];
    }

    return 1;
}

sub bounding_box_of_points {
    my ($d, @points) = @_;

    my @bb;

    while (my @p = splice @points, 0, $d) {
        @bb = bounding_box($d, @p, @bb); # Defined below.
    }

    return @bb;
}


sub clockwise {
  my ($x0, $y0, $x1, $y1, $x2, $y2) = @_;
  return ($x2 - $x0) * ($y1 - $y0) - ( $x1 - $x0) * ($y2 - $y0);
}
sub closest_points {
    my ( @p ) = @_;

    return () unless @p and @p % 2 == 0;

    my $unsorted_x = [ map { $p[ 2 * $_     ] } 0..$#p/2 ];
    my $unsorted_y = [ map { $p[ 2 * $_ + 1 ] } 0..$#p/2 ];

    # Compute the permutation and ordinal indices.

    # X Permutation Index.
    #
    # If @$unsorted_x is (18, 7, 25, 11), @xpi will be (1, 3, 0, 2),
    # e.g., $xpi[0] == 1 meaning that the $sorted_x[0] is in
    # $unsorted_x->[1].
    #
    # We do this because we may now sort @$unsorted_x to @sorted_x
    # and can still restore the original ordering as @sorted_x[@xpi].
    # This is needed because we will want to sort the points by x and y
    # but might also want to identify the result by the original point
    # indices: "the 12th point and the 45th point are the closest pair".

    my @xpi = sort { $unsorted_x->[ $a ] <=> $unsorted_x->[ $b ] }
                   0..$#$unsorted_x;

    # Y Permutation Index.
    #
    my @ypi = sort { $unsorted_y->[ $a ] <=> $unsorted_y->[ $b ] }
                   0..$#$unsorted_y;

    # Y Ordinal Index.
    #
    # The ordinal index is the inverse of the permutation index: If
    # @$unsorted_y is (16, 3, 42, 10) and @ypi is (1, 3, 0, 2), @yoi
    # will be (2, 0, 3, 1), e.g. $yoi[0] == 1 meaning that
    # $unsorted_y->[0] is the $sorted_y[1].

    my @yoi;
    @yoi[ @ypi ] = 0..$#ypi;

    # Recurse to find the closest points.
    my ( $p, $q, $d ) = _closest_points_recurse( [ @$unsorted_x[@xpi] ],
                                                  [ @$unsorted_y[@xpi] ],
                                                  \@xpi, \@yoi, 0, $#xpi
                                                );

    my $pi = $xpi[ $p ];                           # Permute back.
    my $qi = $xpi[ $q ];

    ( $pi, $qi ) = ( $qi, $pi ) if $pi > $qi;      # Smaller id first.
    return ( $pi, $qi, $d );
}

sub _closest_points_recurse {
    my ( $x, $y, $xpi, $yoi, $x_l, $x_r ) = @_;

    # $x, $y:  array references to the x- and y-coordinates of the points
    # $xpi:    x permutation indices: computed by closest_points_recurse()
    # $yoi:    y ordering indices: computed by closest_points_recurse()
    # $x_l:    the left  bound of the currently interesting point set
    # $x_r:    the right bound of the currently interesting point set
    #          That is, only points $x->[$x_l..$x_r] and $y->[$x_l..$x_r]
    #          will be inspected.

    my $d;     # The minimum distance found.
    my $p;     # The index of the other end of the minimum distance.
    my $q;     # Ditto.

    my $N = $x_r - $x_l + 1;      # Number of interesting points.

    if ( $N > 3 ) {               # We have lots of points.  Recurse!
        my $x_lr = int( ( $x_l + $x_r ) / 2 ); # Right bound of left half.
        my $x_rl = $x_lr + 1;                  # Left bound of right half.

        # First recurse to find out how the halves do.

        my ( $p1, $q1, $d1 ) =
            _closest_points_recurse( $x, $y, $xpi, $yoi, $x_l, $x_lr );
        my ( $p2, $q2, $d2 ) =
            _closest_points_recurse( $x, $y, $xpi, $yoi, $x_rl, $x_r );

        # Then merge the halves' results.

        # Update the $d, $p, $q to be the closest distance
        # and the indices of the closest pair of points so far.

        if ( $d1 < $d2 ) { $d = $d1;  $p = $p1;  $q = $q1 }
        else             { $d = $d2;  $p = $p2;  $q = $q2 }

        # Then check the straddling area.

        # The x-coordinate halfway between the left and right halves.
        my $x_d = ( $x->[ $x_lr ] + $x->[ $x_rl ] ) / 2;

        # The indices of the "potential" points: those point pairs
        # that straddle the area and have the potential to be closer
        # to each other than the closest pair so far.
        #
        my @xi;

        # Find the potential points from the left half.

        # The left bound of the left segment with potential points.
        my $x_ll;

        if ( $x_lr == $x_l ) { $x_ll = $x_l }
        else {                                     # Binary search.
            my $x_ll_lo = $x_l;
            my $x_ll_hi = $x_lr;
            do { $x_ll = int( ( $x_ll_lo + $x_ll_hi ) / 2 );
                 if ( $x_d - $x->[ $x_ll ] > $d ) {
                    $x_ll_lo = $x_ll + 1;
                 } elsif ( $x_d - $x->[ $x_ll ] < $d ) {
                    $x_ll_hi = $x_ll - 1;
                 }
            } until $x_ll_lo > $x_ll_hi
                or ( $x_d - $x->[ $x_ll ] < $d
                     and ( $x_ll == 0 or
                           $x_d - $x->[ $x_ll - 1 ] > $d ) );
        }
        push @xi, $x_ll..$x_lr;

        # Find the potential points from the right half.

        # The right bound of the right segment with potential points.
        my $x_rr;

        if ( $x_rl == $x_r ) { $x_rr = $x_r }
        else {                                     # Binary search.
            my $x_rr_lo = $x_rl;
            my $x_rr_hi = $x_r;
            do { $x_rr = int( ( $x_rr_lo + $x_rr_hi ) / 2 );
                 if ( $x->[ $x_rr ] - $x_d > $d ) {
                    $x_rr_hi = $x_rr - 1;
                 } elsif ( $x->[ $x_rr ] - $x_d < $d ) {
                    $x_rr_lo = $x_rr + 1;
                 }
            } until $x_rr_hi < $x_rr_lo
                or ( $x->[ $x_rr ] - $x_d < $d
                     and ( $x_rr == $x_r or
                           $x->[ $x_rr + 1 ] - $x_d > $d ) );
        }
        push @xi, $x_rl..$x_rr;

        # Now we know the potential points.  Are they any good?
        # This gets kind of intense.

        # First sort the points by their original indices.

        my @x_by_y   = @$yoi[ @$xpi[ @xi ] ];
        my @i_x_by_y = sort { $x_by_y[ $a ] <=> $x_by_y[ $b ] }
                       0..$#x_by_y;
        my @xi_by_yi;
        @xi_by_yi[ 0..$#xi ] = @xi[ @i_x_by_y ];

        my @xi_by_y = @$yoi[ @$xpi[ @xi_by_yi ] ];
        my @x_by_yi = @$x[ @xi_by_yi ];
        my @y_by_yi = @$y[ @xi_by_yi ];

        # Inspect each potential pair of points (the first point
        # from the left half, the second point from the right).

        for ( my $i = 0; $i <= $#xi_by_yi; $i++ ) {
            my $i_i = $xi_by_y[ $i ];
            my $x_i = $x_by_yi[ $i ];
            my $y_i = $y_by_yi[ $i ];
            for ( my $j = $i + 1; $j <= $#xi_by_yi; $j++ ) {
                # Skip over points that can't be closer
                # to each other than the current best pair.
                last if $xi_by_y[ $j ] - $i_i > 7; # Too far?
                my $y_j = $y_by_yi[ $j ];
                my $dy = $y_j - $y_i;
                last if $dy > $d;                  # Too tall?
                my $x_j = $x_by_yi[ $j ];
                my $dx = $x_j - $x_i;
                next if abs( $dx ) > $d;           # Too wide?
                # Still here?  We may have a winner.
                # Check the distance and update if so.
                my $d3 = sqrt( $dx**2 + $dy**2 );
                if ( $d3 < $d ) {
                    $d = $d3;
                    $p = $xi_by_yi[ $i ];
                    $q = $xi_by_yi[ $j ];
                }
            }
        }
    } elsif ( $N == 3 ) {      # Just three points?  No need to recurse.
        my $x_m = $x_l + 1;
        # Compare the square sums and leave the sqrt for later.
        my $s1 = ($x->[ $x_l ]-$x->[ $x_m ])**2 +
                 ($y->[ $x_l ]-$y->[ $x_m ])**2;
        my $s2 = ($x->[ $x_m ]-$x->[ $x_r ])**2 +
                 ($y->[ $x_m ]-$y->[ $x_r ])**2;
        my $s3 = ($x->[ $x_l ]-$x->[ $x_r ])**2 +
                 ($y->[ $x_l ]-$y->[ $x_r ])**2;
        if ( $s1 < $s2 ) {
            if ( $s1 < $s3 )  { $d = $s1;  $p = $x_l;  $q = $x_m }
            else              { $d = $s3;  $p = $x_l;  $q = $x_r }
        } elsif ( $s2 < $s3 ) { $d = $s2;  $p = $x_m;  $q = $x_r }
          else                { $d = $s3;  $p = $x_l;  $q = $x_r }

        $d = sqrt $d;
    } elsif ( $N == 2 ) {         # Just two points?  No need to recurse.
        $d = sqrt(($x->[ $x_l ]-$x->[ $x_r ])**2 +
                  ($y->[ $x_l ]-$y->[ $x_r ])**2);
        $p = $x_l;
        $q = $x_r;
    } else {                      # Less than two points?  Strange.
        return ( );
    }

    return ( $p, $q, $d );
}

sub determinant {
  $_[0] * $_[3] - $_[1] * $_[2];
}

sub euclidean {
  my @p = @_;
  my $d = @p / 2;

  return sqrt( ($_[0] - $_[2]) ** 2 + ($_[1] - $_[3]) ** 2) if $d == 2;

  my $S = 0;

  my @p0 = splice @p, 0, $d;

  for(my $i = 0; $i < $d; $i++) {
    my $di = $p0[$i] - $p[$i];
    $S += $di * $di;
  }
  return sqrt($S);
}

sub great_circle_dist {
# dist -- find great-circle distance between two points on earth's surface
#                                                              -*- perl -*-
#
# This code was written in 1998 by Darrell Kindred <dkindred@cmu.edu>.
# I have released it into the public domain.  Do what you like with it,
# but if you make any modifications, please either list your name and
# describe the changes here, or remove my name and address above.
#
# This code is distributed without any warranty, express or implied.
#
# Calculate the great-circle distance and initial heading from one point on
# the Earth to another, given latitude & longitude.
#   usage: dist <location1> <location2>
#
# Calculations assume a spherical Earth with radius 6367 km.
# (I think this should cause no more than 0.2% error.)
#
# Here are some examples of acceptable location formats:
#
#   40:26:46N,79:56:55W
#   40:26:46.302N 79:56:55.903W
#   40026'21"N 79d58'36"W
#   40d 26' 21" N 79d 58' 36" W
#   40.446195N 79.948862W
#   40.446195N -79.948862E
#
# Be sure to quote arguments from the shell as necessary.
#
# For a good discussion of the formula used here for calculating distances,
# as well as several more and less accurate techniques, see
#    http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
# and
#    http://www.best.com/~williams/avform.htm

  use strict;

  if (scalar(@_) == 2) {
    my($lat1,$long1) = &parse_location(shift);
    my($lat2,$long2) = &parse_location(shift);

    my $meters_per_mile = 1609.344;
    my $nautical_miles_per_mile = (5280 / 6076);

    my $dist = &great_circle_distance($lat1,$long1,$lat2,$long2);

    my $heading = &radians_to_degrees(&initial_heading($lat1,$long1, $lat2,$long2));

    print "great-circle distance from ", &loc_to_string($lat1,$long1), " to ", &loc_to_string($lat2,$long2), " is\n";

    # round to 3 significant digits since that's all we are justified
    # in printing given the spherical earth assumption.
    printf("%8g km\n", &round_to_3($dist / 1000));
    printf("%8g miles\n", &round_to_3($dist / $meters_per_mile));
    printf("%8g nautical miles\n", &round_to_3($dist / $meters_per_mile * $nautical_miles_per_mile));

    printf("initial heading: %.0f degrees (%s)\n", $heading, &heading_string($heading));
  }

  else {
    great_circle_dist_err();
  }

  sub great_circle_dist_err {
    print STDERR "$0: two arguments required, ", scalar(@ARGV), " found\n";
    print STDERR "usage: $0 <loc1> <loc2>\n";
    print STDERR "       allowed loc formats: \n";
    print STDERR <<\EndFormats;
          40:26:46N,79:56:55W
          40:26:46.302N 79:56:55.903W
          40026'21"N 79d58'36"W
          40d 26' 21" N 79d 58' 36" W
          40.446195N 79.948862W
          40.446195N -79.948862E
EndFormats
  exit(1);
  }

  # given coordinates of two places in radians, compute distance in meters
  sub great_circle_distance {
    my ($lat1,$long1,$lat2,$long2) = @_;

    # approx radius of Earth in meters.  True radius varies from
    # 6357km (polar) to 6378km (equatorial).
    my $earth_radius = 6367000;

    my $dlon = $long2 - $long1;
    my $dlat = $lat2 - $lat1;
    my $a = (sin($dlat / 2)) ** 2 + cos($lat1) * cos($lat2) * (sin($dlon / 2)) ** 2;
    my $d = 2 * atan2(sqrt($a), sqrt(1 - $a));

    # This is a simpler formula, but it's subject to rounding errors
    # for small distances.  See http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
    # my $d = &acos(sin($lat1) * sin($lat2) + cos($lat1) * cos($lat2) * cos($long1-$long2));

    return $earth_radius * $d;
  }

  # compute the initial bearing (in radians) to get from lat1/long1 to lat2/long2
  sub initial_heading {
    my ($lat1,$long1,$lat2,$long2) = @_;

    BEGIN {
      $::pi = 4 * atan2(1,1);
    }

    # note that this is the same $d calculation as above.
    # duplicated for clarity.
    my $dlon = $long2 - $long1;
    my $dlat = $lat2 - $lat1;
    my $a = (sin($dlat / 2)) ** 2 + cos($lat1) * cos($lat2) * (sin($dlon / 2)) ** 2;
    my $d = 2 * atan2(sqrt($a), sqrt(1 - $a));

    my $heading = acos((sin($lat2) - sin($lat1) * cos($d)) / (sin($d) * cos($lat1)));
    if (sin($long2 - $long1) < 0) {
      $heading = 2 * $::pi - $heading;
    }
    return $heading;
  }

  # return an angle in radians, between 0 and pi, whose cosine is x
  sub acos {
    my($x) = @_;
    die "bad acos argument ($x)\n" if (abs($x) > 1.0);
    return atan2(sqrt(1 - $x * $x), $x);
  }

  # round to 3 significant digits
  sub round_to_3 {
    my ($num) = @_;
    my ($lg,$round);
    if ($num == 0) {
      return 0;
    }
    $lg = int(log(abs($num)) / log(10.0));      # log base 10 of num
    $round = 10 ** ($lg - 2);
    return int($num / $round + 0.5) * $round;
  }

  # round to nearest integer
  sub round_to_int {
    return int(abs($_[0]) + 0.5) * ($_[0] < 0 ? -1 : 1);
  }

  # print a location in a canonical form
  sub loc_to_string {
    my($lat,$long) = @_;

    $lat  = &radians_to_degrees($lat);
    $long  = &radians_to_degrees($long);

    my $ns = "N";
    my $ew = "E";

    if ($lat < 0) {
      $lat = -$lat;
      $ns = "S";
    }
    if ($long < 0) {
      $long = -$long;
      $ew = "W";
    }
    $lat = int($lat * 3600 + 0.5);
    my $lat_string = sprintf("%d:%02d:%02d%s",
                             int($lat/3600),
                             int(($lat - 3600*int($lat / 3600))/60),
                             $lat - 60*int($lat / 60),
                             $ns);
    $long = int($long * 3600 + 0.5);
    my $long_string = sprintf("%d:%02d:%02d%s",
                              int($long/3600),
                              int(($long - 3600*int($long / 3600))/60),
                              $long - 60*int($long / 60),
                              $ew);

    return "$lat_string $long_string";
  }

  # convert a string which looks like "34:45:12N,15:34:10W" into a pair
  # of degrees.  Also accepts "34.233N,90.134E" etc.
  sub parse_location {
    my($str) = @_;
    my($lat,$long);

    if ($str =~ /^([0-9:.\260'"d -]*)([NS])[, ]+([0-9:.\260'"d -]*)([EW])$/i) {
      return undef if (!defined($lat = &parse_degrees($1)));
      $lat *= (($2 eq "N" || $2 eq "n") ? 1.0 : -1.0);
      return undef if (!defined($long = &parse_degrees($3)));
      $long *= (($4 eq "E" || $4 eq "e") ? 1.0 : -1.0);
      return(&degrees_to_radians($lat), &degrees_to_radians($long));
    }
    else {
      return undef;
    }
  }

  # given a bearing in degrees, return a string "north", "southwest",
  # "east-southeast", etc.
  sub heading_string {
    my($deg) = @_;
    my($rounded,$s);
    my(@dirs) = ("north","east","south","west");
    $rounded = &round_to_int($deg / 22.5) % 16;
    if (($rounded % 4) == 0) {
      $s = $dirs[$rounded/4];
    }
    else {
      $s = $dirs[2 * int(((int($rounded / 4) + 1) % 4) / 2)];
      $s .= $dirs[1 + 2 * int($rounded / 8)];
      if ($rounded % 2 == 1) {
        $s = $dirs[&round_to_int($rounded/4) % 4] . "-" . $s;
      }
    }
    return $s;
  }

  BEGIN { $::pi = 4 * atan2(1,1); }
  sub degrees_to_radians {
    return $_[0] * $::pi / 180.0;
  }

  sub radians_to_degrees {
    return $_[0] * 180.0 / $::pi;
  }

  # convert a string like 34:45:12.34 or 38:40 or 34.124 or 34d45'12.34"
  # or 250 02' 30" to degrees
  # (also handles a leading `-')
  sub parse_degrees {
    my($str) = @_;
    my($d,$m,$s,$sign);

    # yeah, this could probably be done with one regexp.
    if ($str =~ /^\s*(-?)([\d.]+)\s*(:|d|\260)\s*([\d.]+)\s*(:|\')\s*([\d.]+)\s*\"?\s*$/) {
      $sign = ($1 eq "-") ? -1.0 : 1.0;
      $d = $2 + 0.0;
      $m = $4 + 0.0;
      $s = $6 + 0.0;
    }
    elsif ($str =~ /^\s*(-?)([\d.]+)\s*(:|d|\260)\s*([\d.]+)(\')?\s*$/) {
      $sign = ($1 eq "-") ? -1.0 : 1.0;
      $d = $2 + 0.0;
      $m = $4 + 0.0;
      $s = 0.0;
    }
    elsif ($str =~ /^\s*(-?)([\d.]+)(d|\260)?\s*$/) {
      $sign = ($1 eq "-") ? -1.0 : 1.0;
      $d = $2 + 0.0;
      $m = 0.0;
      $s = 0.0;
    }
    else {
      die "parse_degrees: can't parse $str\n";
    }
    return ($sign * ($d + ($m / 60.0) + ($s / 3600.0)));
  }
}

sub line_intersection {
  my ( $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3 );

  if ( @_ == 8 ) {
    ( $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3 ) = @_;

    # The bounding boxes chop the lines into line segments.
    # bounding_box() is defined later in this chapter.
    my @box_a = bounding_box( 2, $x0, $y0, $x1, $y1 );
    my @box_b = bounding_box( 2, $x2, $y2, $x3, $y3 );

    # Take this test away and the line segments are
    # turned into lines going from infinite to another.
    # bounding_box_intersect() defined later in this chapter.
    return "out of bounding box" unless bounding_box_intersect( 2, @box_a, @box_b );
  }
  elsif ( @_ == 4 ) { # The parametric form.
    $x0 = $x2 = 0;
    ( $y0, $y2 ) = @_[ 1, 3 ];
    # Need to multiply by 'enough' to get 'far enough'.
    my $abs_y0 = abs $y0;
    my $abs_y2 = abs $y2;
    my $enough = 10 * ( $abs_y0 > $abs_y2 ? $abs_y0 : $abs_y2 );
    $x1 = $x3 = $enough;
    $y1 = $_[0] * $x1 + $y0;
    $y3 = $_[2] * $x2 + $y2;
  }

  my ($x, $y);  # The as-yet-undetermined intersection point.

  my $dy10 = $y1 - $y0; # dyPQ, dxPQ are the coordinate differences
  my $dx10 = $x1 - $x0; # between the points P and Q.
  my $dy32 = $y3 - $y2;
  my $dx32 = $x3 - $x2;

  my $dy10z = abs( $dy10 ) < epsilon; # Is the difference $dy10 "zero"?
  my $dx10z = abs( $dx10 ) < epsilon;
  my $dy32z = abs( $dy32 ) < epsilon;
  my $dx32z = abs( $dx32 ) < epsilon;

  my $dyx10;                            # The slopes.
  my $dyx32;


  $dyx10 = $dy10 / $dx10 unless $dx10z;
  $dyx32 = $dy32 / $dx32 unless $dx32z;

  # Now we know all differences and the slopes;
  # we can detect horizontal/vertical special cases.
  # E.g., slope = 0 means a horizontal line.

  unless ( defined $dyx10 or defined $dyx32 ) {
    return "parallel vertical";
  }
  elsif ( $dy10z and not $dy32z ) { # First line horizontal.
    $y = $y0;
    $x = $x2 + ( $y - $y2 ) * $dx32 / $dy32;
  }
  elsif ( not $dy10z and $dy32z ) { # Second line horizontal.
    $y = $y2;
    $x = $x0 + ( $y - $y0 ) * $dx10 / $dy10;
  }
  elsif ( $dx10z and not $dx32z ) { # First line vertical.
    $x = $x0;
    $y = $y2 + $dyx32 * ( $x - $x2 );
  }
  elsif ( not $dx10z and $dx32z ) { # Second line vertical.
    $x = $x2;
    $y = $y0 + $dyx10 * ( $x - $x0 );
  }
  elsif ( abs( $dyx10 - $dyx32 ) < epsilon ) {
    # The slopes are suspiciously close to each other.
    # Either we have parallel collinear or just parallel lines.

    # The bounding box checks have already weeded the cases
    # "parallel horizontal" and "parallel vertical" away.

    my $ya = $y0 - $dyx10 * $x0;
    my $yb = $y2 - $dyx32 * $x2;

    return "parallel collinear" if abs( $ya - $yb ) < epsilon;
    return "parallel";
  }
  else {
    # None of the special cases matched.
    # We have a "honest" line intersection.

    $x = ($y2 - $y0 + $dyx10*$x0 - $dyx32*$x2)/($dyx10 - $dyx32);
    $y = $y0 + $dyx10 * ($x - $x0);
  }

  my $h10 = $dx10 ? ($x - $x0) / $dx10 : ($dy10 ? ($y - $y0) / $dy10 : 1);
  my $h32 = $dx32 ? ($x - $x2) / $dx32 : ($dy32 ? ($y - $y2) / $dy32 : 1);

  return ($x, $y, $h10 >= 0 && $h10 <= 1 && $h32 >= 0 && $h32 <= 1);
}

sub manhattan {
  my @p = @_;
  my $d = @p / 2;

  my $S = 0;
  my @p0 = splice @p, 0, $d;

  for(my $i = 0; $i < $d; $i++) {
    my $di = $p0[$i] - $p[$i];
    $S += abs $di;
  }
  return $S;
}

sub point_in_polygon {
  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 point_in_quadrangle {
  my ($x, $y, $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3) = @_;

  return point_in_triangle($x,$y,$x0,$y0,$x1,$y1,$x2,$y2)||
	 point_in_triangle($x,$y,$x0,$y0,$x2,$y2,$x3,$y3);
}

sub point_in_triangle {
  my ( $x, $y, $x0, $y0, $x1, $y1, $x2, $y2 ) = @_;

  # clockwise() from earlier in the chapter.
  my $cw0 = clockwise( $x0, $y0, $x1, $y1, $x, $y );
  return 1 if abs( $cw0 ) < epsilon; # On 1st edge.

  my $cw1 = clockwise( $x1, $y1, $x2, $y2, $x, $y );
  return 1 if abs( $cw1 ) < epsilon; # On 2nd edge.

  # Fail if the sign changed.
  return 0 if ( $cw0 < 0 and $cw1 > 0 ) or ( $cw0 > 0 and $cw1 < 0 );

  my $cw2 = clockwise( $x2, $y2, $x0, $y0, $x, $y );
  return 1 if abs( $cw2 ) < epsilon; # On 3rd edge.

  # Fail if the sign changed.
  return 0 if ( $cw0 < 0 and $cw2 > 0 ) or ( $cw0 > 0 and $cw2 < 0 );

  # Jubilate!
  return 1;
}

sub polygon_area {
  my @xy = @_;
  my $A = 0;

  for(my ($xa, $ya) = @xy[-2,-1];my ($xb, $yb) = splice @xy, 0, 2;($xa, $ya) = ($xb, $yb))
  { $A += determinant($xa, $ya, $xb, $yb) }
  return abs $A / 2;
}

sub polygon_perimeter {
  my @xy = @_;
  my  $P = 0;

  for(my ($xa, $ya) = @xy[-2,-1];my ($xb, $yb) = splice @xy, 0, 2;($xa, $ya) = ($xb, $yb))
  { $P += euclidean($xa, $ya, $xb, $yb) }
  return $P;
}

sub triangle_area_heron {
  my ($a, $b, $c);
  if (@_ == 3)
  { ($a, $b, $c) = @_ }
  elsif (@_ == 6) {
    ($a, $b, $c) = (euclidean($_[0], $_[1], $_[2], $_[3]),
		                euclidean($_[2], $_[3], $_[4], $_[5]),
		                euclidean($_[4], $_[5], $_[0], $_[1]), );
  }
  my $s = ($a + $b + $c) / 2;  # semiperimeter
  return sqrt( $s * ($s - $a) * ($s - $b) * ($s - $c));
}

sub triangle_perimeter {
  my ($a, $b, $c);
  if (@_ == 3)
  { ($a, $b, $c) = @_ }
  elsif (@_ == 6) {
    ($a, $b, $c) = (euclidean($_[0], $_[1], $_[2], $_[3]),
		                euclidean($_[2], $_[3], $_[4], $_[5]),
		                euclidean($_[4], $_[5], $_[0], $_[1]), );
  }
  else
  {  print "Illegal number of inputs" }
  return $a + $b + $c if (@_ == 3) or (@_ == 6);
}

sub spatial_extent {
  my ($minx,$miny,$maxx,$maxy,$numrows,$numcols) = @_;
  my $xres = ($maxx - $minx) / $numcols;
  my $yres = ($maxy - $miny) / $numrows;
  return $xres,$yres;
}

sub FileSizeSingleBandRaster {
  my ($numcols, $numrows,$numbytes) = @_;
  return ($numcols*$numrows*$numbytes);
}

sub FileSizeSingleBandRasterHeader {
  my ($numcols, $numrows, $numbytes, $headerbytes) = @_;
  return (($numcols*$numrows*$numbytes) + $headerbytes);
}

sub FileSizeMultiBandRaster {
  my ($numcols, $numrows, $numbytes, $numbands) = @_;
  return ($numcols*$numrows*$numbytes*$numbands);
}

sub FileSizeMultiBandRasterHeader {
  my ($numcols, $numrows, $numbytes, $numbands, $headerbytes) = @_;
  return (($numcols*$numrows*$numbytes*$numbands) + $headerbytes);
}

sub FileSizeMultiBandRasterMultiHeader {
  my ($numcols, $numrows, $numbytes, $numbands, $headerbytes) = @_;
  return (($numcols * $numbands * $numrows * $numbytes) + ($headerbytes * $numbands));
}
1;

__END__


=head1 NAME

  geoalglib.pm - Geometric Algorithim Library / Perl module.

=head1 VERSION

  1.0

=head1 SYNOPSIS

  ./geoalglib -s 45.03N -75.45W 46.09N -80.32W
  great-circle distance from 45:01:48N 75:27:00E to 46:05:24N 80:19:12E is
     397 km
     247 miles
     214 nautical miles
  initial heading: 71 degrees (east-northeast)

=head1 DESCRIPTION

  geoalglib.pm provides numerous geometric algorithms for everyday use in a Perl environment.
  The accompanying geoalglib Perl script gives a useful API into the module.

=head1 PREREQUISITES / MODULES

  None.

=head1 OPTIONS / FLAGS / ARGUMENTS

  A full listing below:

Usage: geoalglib [flag] [inputs]

  Geometric Algorithms

  flag function                     inputs
  -e   Euclidean (direct) distance  x1 y1 x2 y2
  -m   Manhattan (grid) distance    x1 y1 x2 y2
  -s   Great Circle distance        x1 y1 x2 y2
  -c   Closest pair of points       x1 y1 x2 y2 x3 y3 ...
  -at  Area of triangle             len1 len2 len3 OR x1 y1 x2 y2 x3 y3
  -pt  Perimeter of triangle        len1 len2 len3 OR x1 y1 x2 y2 x3 y3
  -ap  Area of polygon              x1 y1 x2 y2 x3 y3 ...
  -pp  Perimeter of polygon         x1 y1 x2 y2 x3 y3 ...
  -pip Point in polygon             point x y x1 y1 x2 y2 x3 y3 ...
  -pit Point in triangle            point x y x1 y1 x2 y2 x3 y3
  -vd  Vector direction             x1 y1 x2 y2 x3 y3
  -piq Point in a quadrangle        point x y x1 y1 x2 y2 x3 y3 ...
  -bb  Bounding box of object       x1 y1 x2 y2 x3 y3 ...
  -li  Line intersection            x1 y1 x2 y2 x3 y3 x4 y4
  -se  Compute spatial extent       minx miny maxx maxy numrows numcols

  Computing File Sizes

  -fssbr   Single band raster                          numcols numrows numbytes
  -fssbrh  Single band raster with header              numcols numrows numbytes headerbytes
  -fsmbr   Multi-band raster                           numcols numrows numbytes numbands
  -fsmbrh  Multi-band raster with header               numcols numrows numbytes numbands headerbytes
  -fsmbrmh Multi-band raster with multi headers        numcols numrows numbytes numbands headerbytes

=head1 AUTHOR

  Tom Kralidis
  http://www.kralidis.ca/gis/

=cut