1 /*
2 ****************************************************************************
3 *
4 * MODULE: Vector library
5 *
6 * AUTHOR(S): Original author CERL, probably Dave Gerdes.
7 * Update to GRASS 5.7 Radim Blazek.
8 *
9 * PURPOSE: Lower level functions for reading/writing/manipulating vectors.
10 *
11 * COPYRIGHT: (C) 2001 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18 /* @(#)prune.c 3.0 2/19/98 */
19 /* by Michel Wurtz for GRASS 4.2.1 - <mw@engees.u-strasbg.fr>
20 * This is a complete rewriting of the previous dig_prune subroutine.
21 * The goal remains : it resamples a dense string of x,y coordinates to
22 * produce a set of coordinates that approaches hand digitizing.
23 * That is, the density of points is very low on straight lines, and
24 * highest on tight curves.
25 *
26 * The algorithm used is very different, and based on the suppression
27 * of intermediate points, when they are closer than thresh from a
28 * moving straight line.
29 *
30 * The distance between a point M -> ->
31 * and a AD segment is given || AM ^ AD ||
32 * by the following formula : d = ---------------
33 * ->
34 * || AD ||
35 *
36 * -> -> ->
37 * When comparing || AM ^ AD || and t = thresh * || AD ||
38 *
39 * -> -> -> ->
40 * we call sqdist = | AM ^ AD | = | OA ^ OM + beta |
41 *
42 * -> ->
43 * with beta = OA ^ OD
44 *
45 * The implementation is based on an old integer routine (optimised
46 * for machine without math coprocessor), itself inspired by a PL/1
47 * routine written after a fortran program on some prehistoric
48 * hardware (IBM 360 probably). Yeah, I'm older than before :-)
49 *
50 * The algorithm used doesn't eliminate "duplicate" points (following
51 * points with same coordinates). So we should clean the set of points
52 * before. As a side effect, dig_prune can be called with a null thresh
53 * value. In this case only cleaning is made. The command v.prune is to
54 * be modified accordingly.
55 *
56 * Another important note : Don't try too big threshold, this subroutine
57 * may produce strange things with some pattern (mainly loops, or crossing
58 * of level curves): Try the set { 0,0 -5,0 -4,10 -6,20 -5,30 -5,20 10,10}
59 * with a thershold of 5. This isn't a programmation, but a conceptal bug ;-)
60 *
61 * Input parameters :
62 * points->x, ->y - double precision sets of coordinates.
63 * points->n_points - the total number of points in the set.
64 * thresh - the distance that a string must wander from a straight
65 * line before another point is selected.
66 *
67 * Value returned : - the final number of points in the pruned set.
68 */
69
70 #include <stdio.h>
71 #include <grass/vector.h>
72 #include <math.h>
73
dig_prune(struct line_pnts * points,double thresh)74 int dig_prune(struct line_pnts *points, double thresh)
75 {
76 double *ox, *oy, *nx, *ny;
77 double cur_x, cur_y;
78 int o_num;
79 int n_num; /* points left */
80 int at_num;
81 int ij = 0, /* position of farest point */
82 ja, jd, i, j, k, n, inu, it; /* indicateur de parcours du segment */
83
84 double sqdist; /* square of distance */
85 double fpdist; /* square of distance from chord to farest point */
86 double t, beta; /* as explained in commented algorithm */
87
88 double dx, dy; /* temporary variables */
89
90 double sx[18], sy[18]; /* temporary table for processing points */
91 int nt[17], nu[17];
92
93 /* nothing to do if less than 3 points ! */
94 if (points->n_points <= 2)
95 return (points->n_points);
96
97 ox = points->x;
98 oy = points->y;
99 nx = points->x;
100 ny = points->y;
101
102 o_num = points->n_points;
103 n_num = 0;
104
105 /* Eliminate duplicate points */
106
107 at_num = 0;
108 while (at_num < o_num) {
109 if (nx != ox) {
110 *nx = *ox++;
111 *ny = *oy++;
112 }
113 else {
114 ox++;
115 oy++;
116 }
117 cur_x = *nx++;
118 cur_y = *ny++;
119 n_num++;
120 at_num++;
121
122 while (*ox == cur_x && *oy == cur_y) {
123 if (at_num == o_num)
124 break;
125 at_num++;
126 ox++;
127 oy++;
128 }
129 }
130
131 /* Return if less than 3 points left. When all points are identical,
132 * output only one point (is this valid for calling function ?) */
133
134 if (n_num <= 2)
135 return n_num;
136
137 if (thresh == 0.0) /* Thresh is null, nothing more to do */
138 return n_num;
139
140 /* some (re)initialisations */
141
142 o_num = n_num;
143 ox = points->x;
144 oy = points->y;
145
146 sx[0] = ox[0];
147 sy[0] = oy[0];
148 n_num = 1;
149 at_num = 2;
150 k = 1;
151 sx[1] = ox[1];
152 sy[1] = oy[1];
153 nu[0] = 9;
154 nu[1] = 0;
155 inu = 2;
156
157 while (at_num < o_num) { /* Position of last point to be */
158 if (o_num - at_num > 14) /* processed in a segment. */
159 n = at_num + 9; /* There must be at least 6 points */
160 else /* in the current segment. */
161 n = o_num;
162 sx[0] = sx[nu[1]]; /* Last point written becomes */
163 sy[0] = sy[nu[1]]; /* first of new segment. */
164 if (inu > 1) { /* One point was keeped in the *//* previous segment : */
165 sx[1] = sx[k]; /* Last point of the old segment */
166 sy[1] = sy[k]; /* becomes second of the new one. */
167 k = 1;
168 }
169 else { /* No point keeped : farest point */
170 sx[1] = sx[ij]; /* is loaded in second position */
171 sy[1] = sy[ij]; /* to avoid cutting lines with */
172 sx[2] = sx[k]; /* small cuvature. */
173 sy[2] = sy[k]; /* First point of previous segment */
174 k = 2; /* becomes the third one. */
175 }
176 /* Loding remaining points */
177 for (j = at_num; j < n; j++) {
178 k++;
179 sx[k] = ox[j];
180 sy[k] = oy[j];
181 }
182
183 jd = 0;
184 ja = k;
185 nt[0] = 0;
186 nu[0] = k;
187 inu = 0;
188 it = 0;
189 for (;;) {
190 if (jd + 1 == ja) /* Exploration of segment terminated */
191 goto endseg;
192
193 dx = sx[ja] - sx[jd];
194 dy = sy[ja] - sy[jd];
195 t = thresh * hypot(dx, dy);
196 beta = sx[jd] * sy[ja] - sx[ja] * sy[jd];
197
198 /* Initializing ij, we don't take 0 as initial value
199 ** for fpdist, in case ja and jd are the same
200 */
201 ij = (ja + jd + 1) >> 1;
202 fpdist = 1.0;
203
204 for (j = jd + 1; j < ja; j++) {
205 sqdist = fabs(dx * sy[j] - dy * sx[j] + beta);
206 if (sqdist > fpdist) {
207 ij = j;
208 fpdist = sqdist;
209 }
210 }
211 if (fpdist > t) { /* We found a point to be keeped *//* Restart from farest point */
212 jd = ij;
213 nt[++it] = ij;
214 }
215 else
216 endseg:{ /* All points are inside threshold. */
217 /* Former start becomes new end */
218 nu[++inu] = jd;
219 if (--it < 0)
220 break;
221 ja = jd;
222 jd = nt[it];
223 }
224 }
225 for (j = inu - 1; j > 0; j--) { /* Copy of segment's keeped points */
226 i = nu[j];
227 ox[n_num] = sx[i];
228 oy[n_num] = sy[i];
229 n_num++;
230 }
231 at_num = n;
232 }
233 i = nu[0];
234 ox[n_num] = sx[i];
235 oy[n_num] = sy[i];
236 n_num++;
237 return n_num;
238 }
239