GRASS Programmer's Manual
6.4.1(2011)
|
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 }