GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
break_polygons.c
Go to the documentation of this file.
1 
20 #include <stdlib.h>
21 #include <math.h>
22 #include <grass/gis.h>
23 #include <grass/Vect.h>
24 #include <grass/glocale.h>
25 
26 typedef struct
27 {
28  double x, y;
29  double a1, a2; /* angles */
30  char cross; /* 0 - do not break, 1 - break */
31  char used; /* 0 - was not used to break line, 1 - was used to break line
32  * this is stored because points are automaticaly marked as cross, even if not used
33  * later to break lines */
34 } XPNT;
35 
36 static int fpoint;
37 
38 /* Function called from RTreeSearch for point found */
39 void srch(int id, int *arg)
40 {
41  fpoint = id;
42 }
43 
61 void
62 Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
63 {
64  struct line_pnts *BPoints, *Points;
65  struct line_cats *Cats;
66  int i, j, k, ret, ltype, broken, last, nlines;
67  int nbreaks;
68  struct Node *RTree;
69  int apoints, npoints, nallpoints, nmarks;
70  XPNT *XPnts;
71  struct Rect rect;
72  double dx, dy, a1 = 0, a2 = 0;
73  int closed, last_point, cross;
74 
75  RTree = RTreeNewIndex();
76 
77  BPoints = Vect_new_line_struct();
78  Points = Vect_new_line_struct();
79  Cats = Vect_new_cats_struct();
80 
81  nlines = Vect_get_num_lines(Map);
82 
83  G_debug(3, "nlines = %d", nlines);
84  /* Go through all lines in vector, and add each point to structure of points,
85  * if such point already exists check angles of segments and if differ mark for break */
86 
87  apoints = 0;
88  nmarks = 0;
89  npoints = 1; /* index starts from 1 ! */
90  nallpoints = 0;
91  XPnts = NULL;
92 
93  for (i = 1; i <= nlines; i++) {
94  G_percent(i, nlines, 1);
95  G_debug(3, "i = %d", i);
96  if (!Vect_line_alive(Map, i))
97  continue;
98 
99  ltype = Vect_read_line(Map, Points, Cats, i);
100  if (!(ltype & type))
101  continue;
102 
103  /* This would be confused by duplicate coordinates (angle cannot be calculated) ->
104  * prune line first */
105  Vect_line_prune(Points);
106 
107  /* If first and last point are identical it is close polygon, we dont need to register last point
108  * and we can calculate angle for first.
109  * If first and last point are not identical we have to mark for break both */
110  last_point = Points->n_points - 1;
111  if (Points->x[0] == Points->x[last_point] &&
112  Points->y[0] == Points->y[last_point])
113  closed = 1;
114  else
115  closed = 0;
116 
117  for (j = 0; j < Points->n_points; j++) {
118  G_debug(3, "j = %d", j);
119  nallpoints++;
120 
121  if (j == last_point && closed)
122  continue; /* do not register last of close polygon */
123 
124  /* Box */
125  rect.boundary[0] = Points->x[j];
126  rect.boundary[3] = Points->x[j];
127  rect.boundary[1] = Points->y[j];
128  rect.boundary[4] = Points->y[j];
129  rect.boundary[2] = 0;
130  rect.boundary[5] = 0;
131 
132  /* Already in DB? */
133  fpoint = -1;
134  RTreeSearch(RTree, &rect, (void *)srch, 0);
135  G_debug(3, "fpoint = %d", fpoint);
136 
137  if (Points->n_points <= 2 ||
138  (!closed && (j == 0 || j == last_point))) {
139  cross = 1; /* mark for cross in any case */
140  }
141  else { /* calculate angles */
142  cross = 0;
143  if (j == 0 && closed) { /* closed polygon */
144  dx = Points->x[last_point] - Points->x[0];
145  dy = Points->y[last_point] - Points->y[0];
146  a1 = atan2(dy, dx);
147  dx = Points->x[1] - Points->x[0];
148  dy = Points->y[1] - Points->y[0];
149  a2 = atan2(dy, dx);
150  }
151  else {
152  dx = Points->x[j - 1] - Points->x[j];
153  dy = Points->y[j - 1] - Points->y[j];
154  a1 = atan2(dy, dx);
155  dx = Points->x[j + 1] - Points->x[j];
156  dy = Points->y[j + 1] - Points->y[j];
157  a2 = atan2(dy, dx);
158  }
159  }
160 
161  if (fpoint > 0) { /* Found */
162  if (XPnts[fpoint].cross == 1)
163  continue; /* already marked */
164 
165  /* Check angles */
166  if (cross) {
167  XPnts[fpoint].cross = 1;
168  nmarks++;
169  }
170  else {
171  G_debug(3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1,
172  XPnts[fpoint].a1, a2, XPnts[fpoint].a2);
173  if ((a1 == XPnts[fpoint].a1 && a2 == XPnts[fpoint].a2) || (a1 == XPnts[fpoint].a2 && a2 == XPnts[fpoint].a1)) { /* identical */
174 
175  }
176  else {
177  XPnts[fpoint].cross = 1;
178  nmarks++;
179  }
180  }
181  }
182  else {
183  /* Add to tree and to structure */
184  RTreeInsertRect(&rect, npoints, &RTree, 0);
185  if (npoints >= apoints) {
186  apoints += 10000;
187  XPnts =
188  (XPNT *) G_realloc(XPnts,
189  (apoints + 1) * sizeof(XPNT));
190  }
191  XPnts[npoints].x = Points->x[j];
192  XPnts[npoints].y = Points->y[j];
193  XPnts[npoints].used = 0;
194  if (j == 0 || j == (Points->n_points - 1) ||
195  Points->n_points < 3) {
196  XPnts[npoints].a1 = 0;
197  XPnts[npoints].a2 = 0;
198  XPnts[npoints].cross = 1;
199  nmarks++;
200  }
201  else {
202  XPnts[npoints].a1 = a1;
203  XPnts[npoints].a2 = a2;
204  XPnts[npoints].cross = 0;
205  }
206 
207  npoints++;
208  }
209  }
210  }
211 
212  /* G_sleep (10); */
213 
214  nbreaks = 0;
215 
216  /* Second loop through lines (existing when loop is started, no need to process lines written again)
217  * and break at points marked for break */
218  for (i = 1; i <= nlines; i++) {
219  int n_orig_points;
220 
221  G_percent(i, nlines, 1);
222  G_debug(3, "i = %d", i);
223  if (!Vect_line_alive(Map, i))
224  continue;
225 
226  ltype = Vect_read_line(Map, Points, Cats, i);
227  if (!(ltype & type))
228  continue;
229  if (!(ltype & GV_LINES))
230  continue; /* Nonsense to break points */
231 
232  /* Duplicates would result in zero length lines -> prune line first */
233  n_orig_points = Points->n_points;
234  Vect_line_prune(Points);
235 
236  broken = 0;
237  last = 0;
238  G_debug(3, "n_points = %d", Points->n_points);
239  for (j = 1; j < Points->n_points; j++) {
240  G_debug(3, "j = %d", j);
241  nallpoints++;
242 
243  /* Box */
244  rect.boundary[0] = Points->x[j];
245  rect.boundary[3] = Points->x[j];
246  rect.boundary[1] = Points->y[j];
247  rect.boundary[4] = Points->y[j];
248  rect.boundary[2] = 0;
249  rect.boundary[5] = 0;
250 
251  if (Points->n_points <= 1 ||
252  (j == (Points->n_points - 1) && !broken))
253  break;
254  /* One point only or
255  * last point and line is not broken, do nothing */
256 
257  RTreeSearch(RTree, &rect, (void *)srch, 0);
258  G_debug(3, "fpoint = %d", fpoint);
259 
260  if (XPnts[fpoint].cross) { /* realy use to break line */
261  XPnts[fpoint].used = 1;
262  }
263 
264  /* break or write last segment of broken line */
265  if ((j == (Points->n_points - 1) && broken) ||
266  XPnts[fpoint].cross) {
267  Vect_reset_line(BPoints);
268  for (k = last; k <= j; k++) {
269  Vect_append_point(BPoints, Points->x[k], Points->y[k],
270  Points->z[k]);
271  }
272 
273  /* Result may collapse to one point */
274  Vect_line_prune(BPoints);
275  if (BPoints->n_points > 1) {
276  ret = Vect_write_line(Map, ltype, BPoints, Cats);
277  G_debug(3,
278  "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
279  ret, j, Points->n_points, BPoints->n_points);
280  }
281 
282  if (!broken)
283  Vect_delete_line(Map, i); /* not yet deleted */
284 
285  last = j;
286  broken = 1;
287  nbreaks++;
288  }
289  }
290  if (!broken && n_orig_points > Points->n_points) { /* was pruned before -> rewrite */
291  if (Points->n_points > 1) {
292  Vect_rewrite_line(Map, i, ltype, Points, Cats);
293  G_debug(3, "Line %d pruned, npoints = %d", i,
294  Points->n_points);
295  }
296  else {
297  Vect_delete_line(Map, i);
298  G_debug(3, "Line %d was deleted", i);
299  }
300  }
301  else {
302  G_debug(3, "Line %d was not changed", i);
303  }
304  }
305 
306  /* Write points on breaks */
307  if (Err) {
308  Vect_reset_cats(Cats);
309  for (i = 1; i < npoints; i++) {
310  if (XPnts[i].used) {
311  Vect_reset_line(Points);
312  Vect_append_point(Points, XPnts[i].x, XPnts[i].y, 0);
313  Vect_write_line(Err, GV_POINT, Points, Cats);
314  }
315  }
316  }
317 
318  G_free(XPnts);
319  RTreeDestroyNode(RTree);
320 }