GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-f4d8c62acd
prune.c
Go to the documentation of this file.
1 /*****************************************************************************
2  *
3  * MODULE: Vector library
4  *
5  * AUTHOR(S): Original author CERL, probably Dave Gerdes.
6  * Update to GRASS 5.7 Radim Blazek.
7  *
8  * PURPOSE: Lower level functions for reading/writing/manipulating vectors.
9  *
10  * COPYRIGHT: (C) 2001 by the GRASS Development Team
11  *
12  * This program is free software under the GNU General Public
13  * License (>=v2). Read the file COPYING that comes with GRASS
14  * for details.
15  *
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 
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 farthest 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 farthest 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 >
165  1) { /* One point was keeped in the */ /* previous segment : */
166  sx[1] = sx[k]; /* Last point of the old segment */
167  sy[1] = sy[k]; /* becomes second of the new one. */
168  k = 1;
169  }
170  else { /* No point keeped : farthest point */
171  sx[1] = sx[ij]; /* is loaded in second position */
172  sy[1] = sy[ij]; /* to avoid cutting lines with */
173  sx[2] = sx[k]; /* small cuvature. */
174  sy[2] = sy[k]; /* First point of previous segment */
175  k = 2; /* becomes the third one. */
176  }
177  /* Loading remaining points */
178  for (j = at_num; j < n; j++) {
179  k++;
180  sx[k] = ox[j];
181  sy[k] = oy[j];
182  }
183 
184  jd = 0;
185  ja = k;
186  nt[0] = 0;
187  nu[0] = k;
188  inu = 0;
189  it = 0;
190  for (;;) {
191  if (jd + 1 == ja) /* Exploration of segment terminated */
192  goto endseg;
193 
194  dx = sx[ja] - sx[jd];
195  dy = sy[ja] - sy[jd];
196  t = thresh * hypot(dx, dy);
197  beta = sx[jd] * sy[ja] - sx[ja] * sy[jd];
198 
199  /* Initializing ij, we don't take 0 as initial value
200  ** for fpdist, in case ja and jd are the same
201  */
202  ij = (ja + jd + 1) >> 1;
203  fpdist = 1.0;
204 
205  for (j = jd + 1; j < ja; j++) {
206  sqdist = fabs(dx * sy[j] - dy * sx[j] + beta);
207  if (sqdist > fpdist) {
208  ij = j;
209  fpdist = sqdist;
210  }
211  }
212  if (fpdist >
213  t) { /* We found a point to be keeped */ /* Restart from
214  farthest point */
215  jd = ij;
216  nt[++it] = ij;
217  }
218  else
219  endseg : { /* All points are inside threshold. */
220  /* Former start becomes new end */
221  nu[++inu] = jd;
222  if (--it < 0)
223  break;
224  ja = jd;
225  jd = nt[it];
226  }
227  }
228  for (j = inu - 1; j > 0; j--) { /* Copy of segment's keeped points */
229  i = nu[j];
230  ox[n_num] = sx[i];
231  oy[n_num] = sy[i];
232  n_num++;
233  }
234  at_num = n;
235  }
236  i = nu[0];
237  ox[n_num] = sx[i];
238  oy[n_num] = sy[i];
239  n_num++;
240  return n_num;
241 }
double cur_x
Definition: driver/init.c:32
double cur_y
Definition: driver/init.c:33
int dig_prune(struct line_pnts *points, double thresh)
Definition: prune.c:74
double t
Definition: r_raster.c:39
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
double * x
Array of X coordinates.
Definition: dig_structs.h:1655
int n_points
Number of points.
Definition: dig_structs.h:1667