prune.c

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

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5