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 }