```# Douglas - Peucker algorithm # Author. John D. Coryat 01/2007 USNaviguide.com # Adapted from: http://mapserver.gis.umn.edu/community/scripts/thin.pl ## # Douglas-Peucker polyline simplification algorithm. First draws single line # from start to end. Then finds largest deviation from this straight line, and if # greater than tolerance, includes that point, splitting the original line into # two new lines. Repeats recursively for each new line created. ## package USNaviguide_Douglas_Peucker ; require 5.003 ; use strict ; BEGIN { use Exporter ; use vars qw ( \$VERSION @ISA @EXPORT) ; \$VERSION = 1.0 ; @ISA = qw ( Exporter ) ; @EXPORT = qw ( Douglas_Peucker perp_distance haversine_distance_meters angle3points ) ; } # Call as: @Opoints = &Douglas_Peucker( , ) ; # Returns: Array of points # Points Array Format: # ([lat1,lng1],[lat2,lng2],...[latn,lngn]) # 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 ; # 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 == \$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 ; } } if ( \$polygon ) # Check for 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( , , ) ; sub perp_distance # Perpendicular distance in meters { 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 * \$O ; my \$c = \$\$p2 * \$O ; my \$d = \$b - \$c ; my \$e = (\$\$p1 * \$O) - (\$\$p2 * \$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(\$_, sqrt(1 - \$_ * \$_)) ; } } # Calculate Angle in Radians between three points... sub angle3points # Angle between three points in radians { my \$p1 = shift ; my \$p2 = shift ; my \$p3 = shift ; my \$m1 = &slope( \$p2, \$p1 ) ; my \$m2 = &slope( \$p3, \$p1 ) ; return (\$m2 - \$m1) ; sub slope # Slope in radians { my \$p1 = shift ; my \$p2 = shift ; return( sprintf("%0.6f",atan2( (@\$p2 - @\$p1),( @\$p2 - @\$p1 ))) ) ; } } 1; __END__ =head1 SYNOPSIS # Test Program: # Douglas - Peucker Test Program # Author. John D. Coryat 01/2007 USNaviguide.com use strict; my \$infile = \$ARGV ; my \$outfile = \$ARGV ; my \$tolerance = \$ARGV ; my @Ipoints = ( ) ; my @Opoints = ( ) ; my \$data = '' ; if(!\$infile or !\$outfile or !\$tolerance) { print "Usage: douglas-peucker.pl \n"; print "Data format: (lat,lng)\n" ; exit ; } if ( \$tolerance <= 0 ) { print "Tolerance (meters) must be greater than zero.\n" ; exit } if ( !(-s \$infile) ) { print "Input File (\$infile) not found.\n" ; exit } if (-s \$outfile) { print "Output File (\$outfile) exists.\n" ; exit } open IN, \$infile ; while ( \$data = ) { if ( \$data =~ /\((.*),(.*)\)/ ) { push( @Ipoints, [\$1,\$2] ) ; } } close IN ; @Opoints = &Douglas_Peucker( \@Ipoints, \$tolerance ) ; open OUT, ">\$outfile" ; foreach \$data (@Opoints) { print OUT "(\$\$data,\$\$data)\n" ; } close OUT ; print "Input: " . \$#Ipoints . " Output: " . \$#Opoints . " Tolerance: \$tolerance\n" ; =cut ```