/* * $Id: prune.c,v 1.5 2006/02/09 03:08:58 glynn Exp $ * **************************************************************************** * * MODULE: Vector library * * AUTHOR(S): Original author CERL, probably Dave Gerdes. * Update to GRASS 5.7 Radim Blazek. * Update to GRASS 6.2.1 Markus Metz * * PURPOSE: Lower level functions for reading/writing/manipulating vectors. * * COPYRIGHT: (C) 2001 by the GRASS Development Team * * This program is free software under the GNU General Public * License (>=v2). Read the file COPYING that comes with GRASS * for details. * *****************************************************************************/ /* @(#)prune.c 4.0 5/18/07 */ /* by Markus Metz for GRASS 6.2.1 - * This is a complete rewriting of the previous dig_prune subroutine. * The goal remains : it resamples a dense string of x,y coordinates to * produce a set of coordinates that approaches hand digitizing. * That is, the density of points is very low on straight lines, and * highest on tight curves. * * Here we use the Douglas-Peucker algorithm based on the paper * Douglas, D. H., and T. K. Peucker, Algorithms for the reduction * of the number of points required to represent a digitized line * of its caricature, Can. Cartogr., 10, 112-122, 1973. * * The implementation of this algorithm has been realised by * Dr. Gary J. Robinson, Environmental Systems Science Centre, * University of Reading, Reading, UK (gazza@mail.nerc-essc.ac.uk) * for the gshhs_dp tool * * The code implemented here is ported from the gshhs_dp tool version 1.5 * The primary site for gshhs is * ftp://ftp.soest.hawaii.edu/pwessel/gshhs * Source code and data are also available from * http://www.ngdc.noaa.gov/mgg/shorelines/data/gshhs/ * * The authors of the gshhs_dp tool are P. Wessel and W. H. F. Smith * and released it under the terms of the GNU General Public License as published by * the Free Software Foundation; version 2 of the License. * * To make it work in different projections, some simplifications were done. * The earth is now a flat disk, but we can use map units, i.e. both decimal degrees and meters. * * The first part of this code removes duplicate points (points with same coordinates) * and is always executed. Dig_prune can be called with a null thresh value. * In this case only cleaning is made. In fact, many vectors contain duplicate * points and should be pruned first with thresh=0 before doing an actual line or border * generalization/simplification. * * Another important note: Start with very small thresholds and work your way up. * As a rule of thumb, start with 1/10 of the resolution of the vector. * With a too large threshold you damage the vector, i.e. might * get duplicate centroids. * * The threshold: deviation of an intermediate point from a line connecting start and end point. * * * Input parameters : * points->x, ->y - double precision sets of coordinates. * points->n_points - the total number of points in the set. * thresh - the distance of a given point to the line connecting the two adjacent points. * * Value returned : - the final number of points in the pruned set. */ #include #include #include #include /* Douglas Peucker*/ #define VNULL (void *)NULL void *get_memory (void *prev_addr, int n, size_t size, char *progname); int dig_prune (struct line_pnts *points, double thresh) { double *ox, *oy, *nx, *ny; long o_num; /* original number of points */ long n_num; /* points left */ long at_num; /* current position */ /* Douglas Peuker */ int n_stack, start, end, sig, i, k; int *sig_start, *sig_end; /* indices of start&end of working section */ double band_sqr, max_dev_sqr, dev_sqr; /* threshold and current distances */ double x12, y12, d12, x13, y13, d13, x23, y23, d23; /* triangle lengths */ /* nothing to do if less than 3 points ! */ if (points->n_points <= 2) return (points->n_points); o_num = points->n_points; ox = points->x; /* old x coordinate */ oy = points->y; /* old y coordinate */ nx = points->x; /* new x coordinate */ ny = points->y; /* new y coordinate */ n_num = 0; /* Eliminate duplicate points */ at_num = 0; while (at_num < o_num) { *nx = *ox++; /* set new coor to current coor, increase old */ *ny = *oy++; n_num++; at_num++; while (*ox == *nx && *oy == *ny && at_num < o_num) { at_num++; /* duplicate, go to next */ ox++; oy++; } nx++; ny++; } if (thresh == 0.0) /* Thresh is null, nothing more to do */ return n_num; /* Douglas Peucker */ /* check for simple cases */ if (n_num <= 2) /* one or two points, leave */ return n_num; /* more complex case. initialize stack */ /* some (re)initialisations */ o_num = n_num; ox = points->x; oy = points->y; nx = points->x; ny = points->y; sig_start = (int *) get_memory (VNULL, o_num, sizeof (int), "Douglas-Peucker pruning"); sig_end = (int *) get_memory (VNULL, o_num, sizeof (int), "Douglas-Peucker pruning"); /* the actual pruning after Douglas Peucker... */ band_sqr = pow(thresh,2); n_num = 0; sig_start[0] = 0; sig_end[0] = o_num-1; n_stack = 1; /* while the stack is not empty ... */ while ( n_stack > 0 ) { /* ... pop the top-most entries off the stacks */ start = sig_start[n_stack-1]; end = sig_end[n_stack-1]; n_stack--; if ( end - start > 1 ) /* any intermediate points ? */ { /* ... yes, so find most deviant intermediate point to either side of line joining start & end points */ x12 = fabs(ox[end] - ox[start]); y12 = fabs(oy[end] - oy[start]); d12 = pow(x12,2) + pow(y12,2); for ( i = start + 1, sig = start, max_dev_sqr = -1.0; i < end; i++ ) { x13 = fabs(ox[i] - ox[start]); y13 = fabs(oy[i] - oy[start]); x23 = fabs(ox[i] - ox[end]); y23 = fabs(oy[i] - oy[end]); d13 = pow(x13,2) + pow(y13,2); d23 = pow(x23,2) + pow(y23,2); if ( d13 >= ( d12 + d23 ) ) dev_sqr = d23; else if ( d23 >= ( d12 + d13 ) ) dev_sqr = d13; else dev_sqr = pow(( x13 * y12 - y13 * x12 ),2) / d12; if ( dev_sqr > max_dev_sqr ) { sig = i; max_dev_sqr = dev_sqr; } } if ( max_dev_sqr < band_sqr ) /* is there a sig. intermediate point ? */ { /* ... no, so transfer current start point */ nx[n_num] = ox[start]; ny[n_num] = oy[start]; n_num++; } else { /* ... yes, so push two sub-sections on stack for further processing */ n_stack++; sig_start[n_stack-1] = sig; sig_end[n_stack-1] = end; n_stack++; sig_start[n_stack-1] = start; sig_end[n_stack-1] = sig; } } else { /* ... no intermediate points, so transfer current start point */ nx[n_num] = ox[start]; ny[n_num] = oy[start]; n_num++; } } /* transfer last point */ nx[n_num] = ox[o_num-1]; ny[n_num] = oy[o_num-1]; n_num++; free ((void *)sig_start); free ((void *)sig_end); return n_num; } void *get_memory (void *prev_addr, int n, size_t size, char *progname) { void *tmp; if (n == 0) return(VNULL); /* Take care of n = 0 */ if (prev_addr) { if ((tmp = G_realloc ((void *) prev_addr, (size_t) n * (size_t) size)) == VNULL) { fprintf (stderr, "Fatal Error: %s could not reallocate more memory, n = %d\n", progname, n); exit (EXIT_FAILURE); } } else { if ((tmp = G_malloc ((size_t) n * (size_t) size)) == VNULL) { fprintf (stderr, "Fatal Error: %s could not allocate memory, n = %d\n", progname, n); exit (EXIT_FAILURE); } } return (tmp); }