GRASS Programmer's Manual  6.5.svn(2014)-r66266
 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 /* TODO: 3D support
27  *
28  * atan2() gives angle from x-axis
29  * this is unambiguous only in 2D, not in 3D
30  *
31  * one possibility would be to store unit vectors of length 1
32  * in struct XPNT
33  * double a1[3], a2[3];
34  *
35  * length = sqrt(dx * dx + dy * dy + dz * dz);
36  * dx /= length; dy /= length; dz /=length;
37  * a1[0] = dx; a1[1] = dy; a1[2] = dz;
38  *
39  * get second dx, dy, dz
40  * length = sqrt(dx * dx + dy * dy + dz * dz);
41  * dx /= length; dy /= length; dz /=length;
42  * a2[0] = dx; a2[1] = dy; a2[2] = dz;
43  *
44  * equal angles
45  * if (a1[0] == a2[0] && a1[1] == a2[1] && a1[2] == a2[2])
46  *
47  * disadvantage: increased memory consumption
48  *
49  * new function Vect_break_faces() ?
50  *
51  */
52 
53 typedef struct
54 {
55  double x, y;
56  double a1, a2; /* angles */
57  char cross; /* 0 - do not break, 1 - break */
58  char used; /* 0 - was not used to break line, 1 - was used to break line
59  * this is stored because points are automatically marked as cross, even if not used
60  * later to break lines */
61 } XPNT;
62 
63 static int fpoint;
64 
65 /* Function called from RTreeSearch for point found */
66 int srch(int id, int *arg)
67 {
68  fpoint = id;
69 
70  return 0;
71 }
72 
90 void
91 Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
92 {
93  struct line_pnts *BPoints, *Points;
94  struct line_cats *Cats;
95  int i, j, k, ret, ltype, broken, last, nlines;
96  int nbreaks;
97  struct Node *RTree;
98  int apoints, npoints, nallpoints, nmarks;
99  XPNT *XPnts;
100  struct Rect rect;
101  double dx, dy, a1 = 0, a2 = 0;
102  int closed, last_point, cross;
103 
104  RTree = RTreeNewIndex();
105 
106  BPoints = Vect_new_line_struct();
107  Points = Vect_new_line_struct();
108  Cats = Vect_new_cats_struct();
109 
110  nlines = Vect_get_num_lines(Map);
111 
112  G_debug(3, "nlines = %d", nlines);
113  /* Go through all lines in vector, and add each point to structure of points,
114  * if such point already exists check angles of segments and if differ mark for break */
115 
116  apoints = 0;
117  nmarks = 0;
118  npoints = 1; /* index starts from 1 ! */
119  nallpoints = 0;
120  XPnts = NULL;
121 
122  G_verbose_message(_("Break polygons Pass 1: select break points"));
123 
124  for (i = 1; i <= nlines; i++) {
125  G_percent(i, nlines, 1);
126  G_debug(3, "i = %d", i);
127  if (!Vect_line_alive(Map, i))
128  continue;
129 
130  ltype = Vect_read_line(Map, Points, Cats, i);
131  if (!(ltype & type))
132  continue;
133 
134  /* This would be confused by duplicate coordinates (angle cannot be calculated) ->
135  * prune line first */
136  Vect_line_prune(Points);
137 
138  /* If first and last point are identical it is close polygon, we dont need to register last point
139  * and we can calculate angle for first.
140  * If first and last point are not identical we have to mark for break both */
141  last_point = Points->n_points - 1;
142  if (Points->x[0] == Points->x[last_point] &&
143  Points->y[0] == Points->y[last_point])
144  closed = 1;
145  else
146  closed = 0;
147 
148  for (j = 0; j < Points->n_points; j++) {
149  G_debug(3, "j = %d", j);
150  nallpoints++;
151 
152  if (j == last_point && closed)
153  continue; /* do not register last of close polygon */
154 
155  /* Box */
156  rect.boundary[0] = Points->x[j];
157  rect.boundary[3] = Points->x[j];
158  rect.boundary[1] = Points->y[j];
159  rect.boundary[4] = Points->y[j];
160  rect.boundary[2] = 0;
161  rect.boundary[5] = 0;
162 
163  /* Already in DB? */
164  fpoint = -1;
165  RTreeSearch(RTree, &rect, (void *)srch, 0);
166  G_debug(3, "fpoint = %d", fpoint);
167 
168  if (Points->n_points <= 2 ||
169  (!closed && (j == 0 || j == last_point))) {
170  cross = 1; /* mark for cross in any case */
171  }
172  else { /* calculate angles */
173  cross = 0;
174  if (j == 0 && closed) { /* closed polygon */
175  dx = Points->x[last_point] - Points->x[0];
176  dy = Points->y[last_point] - Points->y[0];
177  a1 = atan2(dy, dx);
178  dx = Points->x[1] - Points->x[0];
179  dy = Points->y[1] - Points->y[0];
180  a2 = atan2(dy, dx);
181  }
182  else {
183  dx = Points->x[j - 1] - Points->x[j];
184  dy = Points->y[j - 1] - Points->y[j];
185  a1 = atan2(dy, dx);
186  dx = Points->x[j + 1] - Points->x[j];
187  dy = Points->y[j + 1] - Points->y[j];
188  a2 = atan2(dy, dx);
189  }
190  }
191 
192  if (fpoint > 0) { /* Found */
193  if (XPnts[fpoint].cross == 1)
194  continue; /* already marked */
195 
196  /* Check angles */
197  if (cross) {
198  XPnts[fpoint].cross = 1;
199  nmarks++;
200  }
201  else {
202  G_debug(3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1,
203  XPnts[fpoint].a1, a2, XPnts[fpoint].a2);
204  if ((a1 == XPnts[fpoint].a1 && a2 == XPnts[fpoint].a2) || (a1 == XPnts[fpoint].a2 && a2 == XPnts[fpoint].a1)) { /* identical */
205 
206  }
207  else {
208  XPnts[fpoint].cross = 1;
209  nmarks++;
210  }
211  }
212  }
213  else {
214  /* Add to tree and to structure */
215  RTreeInsertRect(&rect, npoints, &RTree, 0);
216  if (npoints >= apoints) {
217  apoints += 10000;
218  XPnts =
219  (XPNT *) G_realloc(XPnts,
220  (apoints + 1) * sizeof(XPNT));
221  }
222  XPnts[npoints].x = Points->x[j];
223  XPnts[npoints].y = Points->y[j];
224  XPnts[npoints].used = 0;
225  if (j == 0 || j == (Points->n_points - 1) ||
226  Points->n_points < 3) {
227  XPnts[npoints].a1 = 0;
228  XPnts[npoints].a2 = 0;
229  XPnts[npoints].cross = 1;
230  nmarks++;
231  }
232  else {
233  XPnts[npoints].a1 = a1;
234  XPnts[npoints].a2 = a2;
235  XPnts[npoints].cross = 0;
236  }
237 
238  npoints++;
239  }
240  }
241  }
242 
243  /* G_sleep (10); */
244 
245  nbreaks = 0;
246 
247  /* Second loop through lines (existing when loop is started, no need to process lines written again)
248  * and break at points marked for break */
249 
250  G_verbose_message(_("Break polygons Pass 2: break at selected points"));
251 
252  for (i = 1; i <= nlines; i++) {
253  int n_orig_points;
254 
255  G_percent(i, nlines, 1);
256  G_debug(3, "i = %d", i);
257  if (!Vect_line_alive(Map, i))
258  continue;
259 
260  ltype = Vect_read_line(Map, Points, Cats, i);
261  if (!(ltype & type))
262  continue;
263  if (!(ltype & GV_LINES))
264  continue; /* Nonsense to break points */
265 
266  /* Duplicates would result in zero length lines -> prune line first */
267  n_orig_points = Points->n_points;
268  Vect_line_prune(Points);
269 
270  broken = 0;
271  last = 0;
272  G_debug(3, "n_points = %d", Points->n_points);
273  for (j = 1; j < Points->n_points; j++) {
274  G_debug(3, "j = %d", j);
275  nallpoints++;
276 
277  /* Box */
278  rect.boundary[0] = Points->x[j];
279  rect.boundary[3] = Points->x[j];
280  rect.boundary[1] = Points->y[j];
281  rect.boundary[4] = Points->y[j];
282  rect.boundary[2] = 0;
283  rect.boundary[5] = 0;
284 
285  if (Points->n_points <= 1 ||
286  (j == (Points->n_points - 1) && !broken))
287  break;
288  /* One point only or
289  * last point and line is not broken, do nothing */
290 
291  RTreeSearch(RTree, &rect, (void *)srch, 0);
292  G_debug(3, "fpoint = %d", fpoint);
293 
294  if (XPnts[fpoint].cross) { /* realy use to break line */
295  XPnts[fpoint].used = 1;
296  }
297 
298  /* break or write last segment of broken line */
299  if ((j == (Points->n_points - 1) && broken) ||
300  XPnts[fpoint].cross) {
301  Vect_reset_line(BPoints);
302  for (k = last; k <= j; k++) {
303  Vect_append_point(BPoints, Points->x[k], Points->y[k],
304  Points->z[k]);
305  }
306 
307  /* Result may collapse to one point */
308  Vect_line_prune(BPoints);
309  if (BPoints->n_points > 1) {
310  ret = Vect_write_line(Map, ltype, BPoints, Cats);
311  G_debug(3,
312  "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
313  ret, j, Points->n_points, BPoints->n_points);
314  }
315 
316  if (!broken)
317  Vect_delete_line(Map, i); /* not yet deleted */
318 
319  last = j;
320  broken = 1;
321  nbreaks++;
322  }
323  }
324  if (!broken && n_orig_points > Points->n_points) { /* was pruned before -> rewrite */
325  if (Points->n_points > 1) {
326  Vect_rewrite_line(Map, i, ltype, Points, Cats);
327  G_debug(3, "Line %d pruned, npoints = %d", i,
328  Points->n_points);
329  }
330  else {
331  Vect_delete_line(Map, i);
332  G_debug(3, "Line %d was deleted", i);
333  }
334  }
335  else {
336  G_debug(3, "Line %d was not changed", i);
337  }
338  }
339 
340  /* Write points on breaks */
341  if (Err) {
342  Vect_reset_cats(Cats);
343  for (i = 1; i < npoints; i++) {
344  if (XPnts[i].used) {
345  Vect_reset_line(Points);
346  Vect_append_point(Points, XPnts[i].x, XPnts[i].y, 0);
347  Vect_write_line(Err, GV_POINT, Points, Cats);
348  }
349  }
350  }
351 
352  G_free(XPnts);
353  RTreeDestroyNode(RTree);
354 }
RectReal boundary[NUMSIDES]
Definition: index.h:42
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
double y
double x
double a2
void Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
Break polygons in vector map.
struct line_pnts * Vect_new_line_struct()
Creates and initializes a struct line_pnts.
Definition: line.c:57
char used
int RTreeSearch(struct Node *N, struct Rect *R, SearchHitCallback shcb, void *cbarg)
Definition: index.h:56
void RTreeDestroyNode(struct Node *)
Definition: node.c:217
int Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:148
int Vect_reset_cats(struct line_cats *Cats)
Reset category structure to make sure cats structure is clean to be re-used.
int RTreeInsertRect(struct Rect *R, int Tid, struct Node **Root, int Level)
int y
Definition: plot.c:34
double a1
char cross
int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
Appends one point to the end of a line.
Definition: line.c:168
int G_percent(long n, long d, int s)
Print percent complete messages.
Definition: percent.c:63
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
Definition: lib/gis/error.c:95
tuple id
self.OnVectorSurface(event)
Definition: tools.py:3426
int Vect_rewrite_line(struct Map_info *Map, int line, int type, struct line_pnts *points, struct line_cats *cats)
Rewrites feature info at the given offset.
int GV_LINES
Definition: vdigit/main.py:24
int Vect_line_alive(struct Map_info *Map, int line)
Check if feature is alive or dead.
struct line_cats * Vect_new_cats_struct()
Creates and initializes line_cats structure.
return NULL
Definition: dbfopen.c:1394
int Vect_line_prune(struct line_pnts *Points)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:255
int Vect_get_num_lines(struct Map_info *map)
Fetch number of features (points, lines, boundaries, centroids) in vector map.
Definition: level_two.c:69
tuple Map
Definition: render.py:1310
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
int Vect_delete_line(struct Map_info *Map, int line)
Delete feature.
struct Node * RTreeNewIndex(void)
long Vect_write_line(struct Map_info *Map, int type, struct line_pnts *points, struct line_cats *cats)
Writes new feature to the end of file (table)
Definition: index.h:40
int srch(int id, int *arg)
int Vect_read_line(struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c, int line)
Read vector feature.