User:Mga geo/Poly.pm

package Poly;
# auteur: Marc GAUTHIER
#
# manipulation de lignes brisées/polygone
use strict;
use warnings;
use Carp;
use Data::Dumper;
our @ISA = qw(Exporter);
our @EXPORT_OK = qw(&Douglas_Peucker);
# http://www.usnaviguide.com/douglas-peucker.htm
#http://fr.wikipedia.org/wiki/Algorithme_de_Douglas-Peuker
# http://softsurfer.com/Archive/algorithm_0205/algorithm_0205.htm
sub Douglas_Peucker {
  my $href	= shift ;
  my $tolerance	= shift ;
  my @Ipoints	= @$href ;
  my @Opoints	= ( ) ;
  my @stack	= ( ) ;
  my $fIndex	= 0 ;
  my $fPoint	= '' ;
  my $aIndex	= 0 ;
  my $anchor	= '' ;
  my $max		= 0 ;
  my $maxIndex	= 0 ;
  my $point	= '' ;
  my $dist	= 0 ;
  my $polygon	= 0 ;					# Line Type
  $anchor = $Ipoints[0] ; 				# save first point
  push( @Opoints, $anchor ) ;
  $aIndex = 0 ;						# Anchor Index
# Check for a polygon: At least 4 points and the first point == last point...
  if ( $#Ipoints >= 4 and $Ipoints[0] == $Ipoints[$#Ipoints] ) {
    $fIndex = $#Ipoints - 1 ;				# Start from the next to last point
   $polygon = 1 ;						# It's a polygon
  } else {
    $fIndex = $#Ipoints ;					# It's a path (open polygon)
  }
  push( @stack, $fIndex ) ;

# Douglas - Peucker algorithm...
  while(@stack) {
    $fIndex = $stack[$#stack] ;
    $fPoint = $Ipoints[$fIndex] ;
    $max = $tolerance ;		 			# comparison values
    $maxIndex = 0 ;
 # Process middle points...
    for (($aIndex+1) .. ($fIndex-1)) {
      $point = $Ipoints[$_] ;
      $dist = &perp_distance($anchor, $fPoint, $point);
      if( $dist >= $max ) {
        $max = $dist ;
        $maxIndex = $_;
      }
    }
    if( $maxIndex > 0 ) {
      push( @stack, $maxIndex ) ;
    } else {
      push( @Opoints, $fPoint ) ;
      $anchor = $Ipoints[(pop @stack)] ;
      $aIndex = $fIndex ;
    }
  }
# Check for Polygon
  if ( $polygon ) {
    push( @Opoints, $Ipoints[$#Ipoints] ) ;		# Add the last point
 # Check for collapsed polygons, use original data in that case...

    if( $#Opoints < 4 ) {
      @Opoints = @Ipoints ;
    }
  }
  return ( @Opoints ) ;
}

# Calculate Perpendicular Distance in meters between a line (two points) and a point...
# my $dist = &perp_distance( <line point 1>, <line point 2>, <point> ) ;
sub perp_distance {
  my $lp1	= shift ;
  my $lp2	= shift ;
  my $p		= shift ;
  my $dist	= &haversine_distance_meters( $lp1, $p ) ;
  my $angle	= &angle3points( $lp1, $lp2, $p ) ;
  return ( sprintf("%0.6f", abs($dist * sin($angle)) ) ) ;
}

# Calculate Distance in meters between two points...
sub haversine_distance_meters {
  my $p1	= shift ;
  my $p2	= shift ;
  my $O = 3.141592654/180 ;
  my $b = $$p1[0] * $O ;
  my $c = $$p2[0] * $O ;
  my $d = $b - $c ;
  my $e = ($$p1[1] * $O) - ($$p2[1] * $O) ;
  my $f = 2 * &asin( sqrt( (sin($d/2) ** 2) + cos($b) * cos($c) * (sin($e/2) ** 2)));
  return sprintf("%0.4f",$f * 6378137) ; 		# Return meters
  sub asin {
   atan2($_[0], sqrt(1 - $_[0] * $_[0])) ;
  }
}

# Calculate Angle in Radians between three points...
sub angle3points {
  my $p1	= shift ;
  my $p2	= shift ;
  my $p3 = shift ;
  my $m1 = &slope( $p2, $p1 ) ;
  my $m2 = &slope( $p3, $p1 ) ;

  return ($m2 - $m1) ;
  sub slope {
    my $p1	= shift ;
    my $p2	= shift ;
    return( sprintf("%0.6f",atan2( (@$p2[1] - @$p1[1]),( @$p2[0] - @$p1[0] ))) ) ;
  }
}
1;