prune.c

Go to the documentation of this file.
00001 /*
00002  ****************************************************************************
00003  *
00004  * MODULE:       Vector library 
00005  *              
00006  * AUTHOR(S):    Original author CERL, probably Dave Gerdes.
00007  *               Update to GRASS 5.7 Radim Blazek.
00008  *
00009  * PURPOSE:      Lower level functions for reading/writing/manipulating vectors.
00010  *
00011  * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00012  *
00013  *               This program is free software under the GNU General Public
00014  *              License (>=v2). Read the file COPYING that comes with GRASS
00015  *              for details.
00016  *
00017  *****************************************************************************/
00018 /*  @(#)prune.c 3.0  2/19/98  */
00019 /* by Michel Wurtz for GRASS 4.2.1 - <mw@engees.u-strasbg.fr>
00020  * This is a complete rewriting of the previous dig_prune subroutine.
00021  * The goal remains : it resamples a dense string of x,y coordinates to
00022  * produce a set of coordinates that approaches hand digitizing.
00023  * That is, the density of points is very low on straight lines, and
00024  * highest on tight curves.
00025  *
00026  * The algorithm used is very different, and based on the suppression
00027  * of intermediate points, when they are closer than thresh from a
00028  * moving straight line.  
00029  *
00030  * The distance between a point M                ->   ->
00031  * and a AD segment is given                  || AM ^ AD ||
00032  * by the following formula :            d = ---------------
00033  *                                                  ->
00034  *                                               || AD ||
00035  *
00036  *                     ->   ->                             ->
00037  * When comparing   || AM ^ AD ||   and    t = thresh * || AD ||
00038  *
00039  *                     ->   ->       ->   ->
00040  * we call  sqdist = | AM ^ AD | = | OA ^ OM + beta | 
00041  *
00042  *                  ->   ->
00043  *  with     beta = OA ^ OD 
00044  *
00045  * The implementation is based on an old integer routine (optimised
00046  * for machine without math coprocessor), itself inspired by a PL/1
00047  * routine written after a fortran programm on some prehistoric
00048  * hardware (IBM 360 probably). Yeah, I'm older than before :-)
00049  *
00050  * The algorithm used doesn't eliminate "duplicate" points (following
00051  * points with same coordinates).  So we should clean the set of points
00052  * before.  As a side effect, dig_prune can be called with a null thresh
00053  * value.  In this case only cleaning is made. The command v.prune is to
00054  * be modified accordingly.
00055  *
00056  * Another important note : Don't try too big threshold, this subroutine
00057  * may produce strange things with some pattern (mainly loops, or crossing
00058  * of level curves): Try the set { 0,0 -5,0 -4,10 -6,20 -5,30 -5,20 10,10}
00059  * with a thershold of 5. This isn't a programmation, but a conceptal bug ;-)
00060  *
00061  * Input parameters :
00062  * points->x, ->y   - double precision sets of coordinates.
00063  * points->n_points - the total number of points in the set.
00064  * thresh           - the distance that a string must wander from a straight
00065  *                    line before another point is selected.
00066  *
00067  * Value returned : - the final number of points in the pruned set.
00068  */
00069 
00070 #include <stdio.h>
00071 #include <grass/Vect.h>
00072 #include <math.h>
00073 
00074 int dig_prune(struct line_pnts *points, double thresh)
00075 {
00076     double *ox, *oy, *nx, *ny;
00077     double cur_x, cur_y;
00078     int o_num;
00079     int n_num;                  /* points left */
00080     int at_num;
00081     int ij = 0,                 /* position of farest point */
00082         ja, jd, i, j, k, n, inu, it;    /* indicateur de parcours du segment */
00083 
00084     double sqdist;              /* square of distance */
00085     double fpdist;              /* square of distance from chord to farest point */
00086     double t, beta;             /* as explained in commented algorithm  */
00087 
00088     double dx, dy;              /* temporary variables */
00089 
00090     double sx[18], sy[18];      /* temporary table for processing points */
00091     int nt[17], nu[17];
00092 
00093     /* nothing to do if less than 3 points ! */
00094     if (points->n_points <= 2)
00095         return (points->n_points);
00096 
00097     ox = points->x;
00098     oy = points->y;
00099     nx = points->x;
00100     ny = points->y;
00101 
00102     o_num = points->n_points;
00103     n_num = 0;
00104 
00105     /* Eliminate duplicate points */
00106 
00107     at_num = 0;
00108     while (at_num < o_num) {
00109         if (nx != ox) {
00110             *nx = *ox++;
00111             *ny = *oy++;
00112         }
00113         else {
00114             ox++;
00115             oy++;
00116         }
00117         cur_x = *nx++;
00118         cur_y = *ny++;
00119         n_num++;
00120         at_num++;
00121 
00122         while (*ox == cur_x && *oy == cur_y) {
00123             if (at_num == o_num)
00124                 break;
00125             at_num++;
00126             ox++;
00127             oy++;
00128         }
00129     }
00130 
00131     /*  Return if less than 3 points left.  When all points are identical,
00132      *  output only one point (is this valid for calling function ?) */
00133 
00134     if (n_num <= 2)
00135         return n_num;
00136 
00137     if (thresh == 0.0)          /* Thresh is null, nothing more to do */
00138         return n_num;
00139 
00140     /* some (re)initialisations */
00141 
00142     o_num = n_num;
00143     ox = points->x;
00144     oy = points->y;
00145 
00146     sx[0] = ox[0];
00147     sy[0] = oy[0];
00148     n_num = 1;
00149     at_num = 2;
00150     k = 1;
00151     sx[1] = ox[1];
00152     sy[1] = oy[1];
00153     nu[0] = 9;
00154     nu[1] = 0;
00155     inu = 2;
00156 
00157     while (at_num < o_num) {    /* Position of last point to be    */
00158         if (o_num - at_num > 14)        /* processed in a segment.         */
00159             n = at_num + 9;     /* There must be at least 6 points */
00160         else                    /* in the current segment.         */
00161             n = o_num;
00162         sx[0] = sx[nu[1]];      /* Last point written becomes      */
00163         sy[0] = sy[nu[1]];      /* first of new segment.           */
00164         if (inu > 1) {          /* One point was keeped in the     *//* previous segment :              */
00165             sx[1] = sx[k];      /* Last point of the old segment   */
00166             sy[1] = sy[k];      /* becomes second of the new one.  */
00167             k = 1;
00168         }
00169         else {                  /* No point keeped : farest point  */
00170             sx[1] = sx[ij];     /* is loaded in second position    */
00171             sy[1] = sy[ij];     /* to avoid cutting lines with     */
00172             sx[2] = sx[k];      /* small cuvature.                 */
00173             sy[2] = sy[k];      /* First point of previous segment */
00174             k = 2;              /* becomes the third one.          */
00175         }
00176         /* Loding remaining points         */
00177         for (j = at_num; j < n; j++) {
00178             k++;
00179             sx[k] = ox[j];
00180             sy[k] = oy[j];
00181         }
00182 
00183         jd = 0;
00184         ja = k;
00185         nt[0] = 0;
00186         nu[0] = k;
00187         inu = 0;
00188         it = 0;
00189         for (;;) {
00190             if (jd + 1 == ja)   /* Exploration of segment terminated */
00191                 goto endseg;
00192 
00193             dx = sx[ja] - sx[jd];
00194             dy = sy[ja] - sy[jd];
00195             t = thresh * hypot(dx, dy);
00196             beta = sx[jd] * sy[ja] - sx[ja] * sy[jd];
00197 
00198             /* Initializing ij, we don't take 0 as initial value
00199              ** for fpdist, in case ja and jd are the same
00200              */
00201             ij = (ja + jd + 1) >> 1;
00202             fpdist = 1.0;
00203 
00204             for (j = jd + 1; j < ja; j++) {
00205                 sqdist = fabs(dx * sy[j] - dy * sx[j] + beta);
00206                 if (sqdist > fpdist) {
00207                     ij = j;
00208                     fpdist = sqdist;
00209                 }
00210             }
00211             if (fpdist > t) {   /* We found a point to be keeped    *//* Restart from farest point        */
00212                 jd = ij;
00213                 nt[++it] = ij;
00214             }
00215             else
00216           endseg:{              /* All points are inside threshold. */
00217                 /* Former start becomes new end     */
00218                 nu[++inu] = jd;
00219                 if (--it < 0)
00220                     break;
00221                 ja = jd;
00222                 jd = nt[it];
00223             }
00224         }
00225         for (j = inu - 1; j > 0; j--) { /* Copy of segment's keeped points  */
00226             i = nu[j];
00227             ox[n_num] = sx[i];
00228             oy[n_num] = sy[i];
00229             n_num++;
00230         }
00231         at_num = n;
00232     }
00233     i = nu[0];
00234     ox[n_num] = sx[i];
00235     oy[n_num] = sy[i];
00236     n_num++;
00237     return n_num;
00238 }
Generated on Tue Apr 6 13:28:10 2010 for GRASS Programmer's Manual by  doxygen 1.6.3