GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-56b21d1bb
prune.c
Go to the documentation of this file.
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  * ->
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 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 }
double cur_y
Definition: driver/init.c:33
double cur_x
Definition: driver/init.c:32
int n_points
Number of points.
Definition: dig_structs.h:1692
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
double t
Definition: r_raster.c:39
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
int dig_prune(struct line_pnts *points, double thresh)
Definition: prune.c:74