1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <grass/gis.h>
4 #include <grass/vector.h>
5 #include <grass/glocale.h>
6 #include "defs.h"
7
next_dist(int line,int side,double mf)8 static int next_dist(int line, int side, double mf)
9 {
10 double d, dist, nextdist, totaldist;
11 int nextline, nlines;
12 int node;
13 static struct line_pnts *Points = NULL;
14
15 G_debug(3, "next_dist()");
16
17 if (!Points)
18 Points = Vect_new_line_struct();
19
20 Vect_read_line(&Out, Points, NULL, abs(line));
21 dist = Vect_line_length(Points);
22
23 if (line < 0)
24 Vect_get_line_nodes(&Out, -line, &node, NULL);
25 else
26 Vect_get_line_nodes(&Out, line, NULL, &node);
27
28 nlines = Vect_get_node_n_lines(&Out, node);
29
30 if (nlines == 1)
31 return 1;
32
33 nextdist = totaldist = 0.;
34 while (nlines > 1) {
35 nextline = dig_angle_next_line(&(Out.plus), -line, side, GV_LINE, NULL);
36 Vect_read_line(&Out, Points, NULL, abs(nextline));
37 d = Vect_line_length(Points);
38 nextdist += d;
39 totaldist += d;
40
41 if (nextline < 0)
42 Vect_get_line_nodes(&Out, -nextline, &node, NULL);
43 else
44 Vect_get_line_nodes(&Out, nextline, NULL, &node);
45
46 nlines = Vect_get_node_n_lines(&Out, node);
47 if (nlines > 2)
48 nextdist = 0.;
49 line = nextline;
50 }
51
52 if (totaldist > nextdist && dist > nextdist)
53 return (totaldist < mf * dist);
54
55 return (dist > nextdist);
56 }
57
loop_test(int line,int node,struct line_pnts * Points,double l,double mf)58 static int loop_test(int line, int node, struct line_pnts *Points,
59 double l, double mf)
60 {
61 int line1, line2, nextline;
62 int n1, n2, nspikes, nout;
63 double l1, minl, maxl;
64
65 if (Vect_get_node_n_lines(&Out, node) != 3)
66 return 0;
67
68 line1 = line2 = 0;
69 line1 = Vect_get_node_line(&Out, node, 0);
70 if (abs(line1) == abs(line))
71 line1 = 0;
72 line2 = Vect_get_node_line(&Out, node, 1);
73 if (abs(line2) == abs(line))
74 line2 = 0;
75 if (line1 == 0) {
76 line1 = line2;
77 line2 = 0;
78 }
79 if (line2 == 0)
80 line2 = Vect_get_node_line(&Out, node, 2);
81
82 if (abs(line1) == abs(line2))
83 return 1;
84
85 nextline = dig_angle_next_line(&(Out.plus), -line, GV_LEFT, GV_LINE, NULL);
86 line1 = nextline;
87
88 nspikes = 1;
89 nout = 0;
90 minl = 0;
91 maxl = 0;
92 do {
93 if (nextline < 0)
94 Vect_get_line_nodes(&Out, -nextline, &n1, NULL);
95 else
96 Vect_get_line_nodes(&Out, nextline, NULL, &n1);
97
98 if (Vect_get_node_n_lines(&Out, n1) == 1)
99 return 0;
100
101 if (n1 != node && Vect_get_node_n_lines(&Out, n1) > 2) {
102 nspikes++;
103 line2 = dig_angle_next_line(&(Out.plus), -nextline, GV_LEFT, GV_LINE, NULL);
104
105 if (line2 < 0)
106 Vect_get_line_nodes(&Out, -line2, &n2, NULL);
107 else
108 Vect_get_line_nodes(&Out, line2, NULL, &n2);
109
110 if (Vect_get_node_n_lines(&Out, n2) == 1) {
111 Vect_read_line(&Out, Points, NULL, abs(line2));
112 l1 = Vect_line_length(Points);
113 if (minl == 0 || minl > l1) {
114 minl = l1;
115 }
116 if (maxl < l1) {
117 maxl = l1;
118 }
119 }
120 else
121 nout++;
122 }
123
124 nextline = dig_angle_next_line(&(Out.plus), -nextline, GV_RIGHT, GV_LINE, NULL);
125
126 } while (abs(nextline) != abs(line1));
127
128 if (minl == 0)
129 minl = l;
130 if (maxl == 0)
131 maxl = mf * l;
132
133 nspikes -= nout;
134
135 if (mf > 1) {
136 return (nspikes < 3 || (nspikes == 3 && (l > minl || mf * l > maxl)));
137
138 if (l > minl)
139 return 1;
140 if (nspikes == 3 && l < minl)
141 return (mf * l > maxl);
142 else
143 return (nspikes < 3);
144 }
145
146 return (nspikes < 3 || (nspikes > 2 && l > minl));
147 }
148
break_loop(int line,int node,struct line_pnts * Points)149 static int break_loop(int line, int node, struct line_pnts *Points)
150 {
151 int line1, line2, firstline, nextline;
152 int n1;
153 double l1, l2;
154
155 if (Vect_get_node_n_lines(&Out, node) != 3)
156 return 0;
157
158 line1 = line2 = 0;
159 line1 = Vect_get_node_line(&Out, node, 0);
160 if (abs(line1) == abs(line))
161 line1 = 0;
162 line2 = Vect_get_node_line(&Out, node, 1);
163 if (abs(line2) == abs(line))
164 line2 = 0;
165 if (line1 == 0) {
166 line1 = line2;
167 line2 = 0;
168 }
169 if (line2 == 0)
170 line2 = Vect_get_node_line(&Out, node, 2);
171
172 if (abs(line1) == abs(line2))
173 return 1;
174
175 nextline = dig_angle_next_line(&(Out.plus), -line, GV_LEFT, GV_LINE, NULL);
176 firstline = nextline;
177
178 do {
179 if (nextline < 0)
180 Vect_get_line_nodes(&Out, -nextline, &n1, NULL);
181 else
182 Vect_get_line_nodes(&Out, nextline, NULL, &n1);
183
184 if (Vect_get_node_n_lines(&Out, n1) == 1)
185 return 0;
186
187 nextline = dig_angle_next_line(&(Out.plus), -nextline, GV_RIGHT, GV_LINE, NULL);
188
189 } while (abs(nextline) != abs(firstline));
190
191 if (abs(nextline) != abs(firstline)) {
192 G_warning("no loop ???");
193 return 0;
194 }
195
196 Vect_read_line(&Out, Points, NULL, abs(line1));
197 l1 = Vect_line_length(Points);
198
199 Vect_read_line(&Out, Points, NULL, abs(line2));
200 l2 = Vect_line_length(Points);
201
202 if (l1 > l2)
203 Vect_delete_line(&Out, abs(line1));
204 else
205 Vect_delete_line(&Out, abs(line2));
206
207 Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
208
209 return 1;
210 }
211
length_test(int line,int node,struct line_pnts * Points,double l,double mf)212 static int length_test(int line, int node, struct line_pnts *Points,
213 double l, double mf)
214 {
215 int line1, line2, nlines;
216 int n1;
217 double l1, l2;
218
219 if (Vect_get_node_n_lines(&Out, node) != 3)
220 return 0;
221
222 line1 = Vect_get_node_line(&Out, node, 0);
223 line2 = 0;
224 if (abs(line1) == abs(line))
225 line1 = 0;
226 line2 = Vect_get_node_line(&Out, node, 1);
227 if (abs(line2) == abs(line))
228 line2 = 0;
229 if (line1 == 0) {
230 line1 = line2;
231 line2 = 0;
232 }
233 if (line2 == 0)
234 line2 = Vect_get_node_line(&Out, node, 2);
235
236 if (line1 < 0)
237 Vect_get_line_nodes(&Out, -line1, &n1, NULL);
238 else
239 Vect_get_line_nodes(&Out, line1, NULL, &n1);
240
241 nlines = Vect_get_node_n_lines(&Out, n1);
242 if (nlines > 1)
243 return 0;
244
245 if (line2 < 0)
246 Vect_get_line_nodes(&Out, -line2, &n1, NULL);
247 else
248 Vect_get_line_nodes(&Out, line2, NULL, &n1);
249
250 nlines = Vect_get_node_n_lines(&Out, n1);
251 if (nlines > 1)
252 return 0;
253
254 Vect_read_line(&Out, Points, NULL, abs(line1));
255 l1 = Vect_line_length(Points);
256
257 Vect_read_line(&Out, Points, NULL, abs(line2));
258 l2 = Vect_line_length(Points);
259
260 if (l1 > mf * l2) {
261 if (mf * l < l1 && l < l2)
262 return 0;
263 }
264 if (l2 > mf * l1) {
265 if (mf * l < l2 && l < l1)
266 return 0;
267 }
268
269 return (mf * l > l1 || mf * l > l2);
270 }
271
272 /* thin the skeletons */
thin_skeleton(double thresh)273 int thin_skeleton(double thresh)
274 {
275 int i;
276 int node, n1, n2;
277 struct line_pnts *Points;
278 struct ilist *list;
279 double l, minl, maxl;
280 int nlines, line, minline, line2;
281 int counter = 1;
282 double morphof = 1.6180339887;
283
284 Points = Vect_new_line_struct();
285 list = Vect_new_list();
286
287 if (thresh < 0)
288 morphof = 1;
289
290 Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
291
292 while (TRUE) {
293 G_verbose_message(_("Pass %d"), counter++);
294 Vect_reset_list(list);
295
296 for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
297 if (!Vect_node_alive(&Out, node))
298 continue;
299 nlines = Vect_get_node_n_lines(&Out, node);
300
301 if (nlines > 1)
302 continue;
303
304 line = Vect_get_node_line(&Out, node, 0);
305 if (line < 0)
306 Vect_get_line_nodes(&Out, -line, &n1, NULL);
307 else
308 Vect_get_line_nodes(&Out, line, NULL, &n1);
309
310 nlines = Vect_get_node_n_lines(&Out, n1);
311
312 if (nlines < 3)
313 continue;
314
315 Vect_read_line(&Out, Points, NULL, abs(line));
316 minline = line;
317 minl = Vect_line_length(Points);
318
319 if (nlines == 3) {
320 if (loop_test(line, n1, Points, minl, morphof))
321 continue;
322 if (length_test(line, n1, Points, minl, morphof))
323 continue;
324 }
325
326 maxl = 0;
327 for (i = 0; i < nlines; i++) {
328 line2 = Vect_get_node_line(&Out, n1, i);
329 if (abs(line2) == abs(minline) || abs(line2) == abs(line))
330 continue;
331 if (line2 < 0)
332 Vect_get_line_nodes(&Out, -line2, &n2, NULL);
333 else
334 Vect_get_line_nodes(&Out, line2, NULL, &n2);
335 if (Vect_get_node_n_lines(&Out, n2) > 1)
336 continue;
337 Vect_read_line(&Out, Points, NULL, abs(line2));
338 l = Vect_line_length(Points);
339 if (minl > l) {
340 minl = l;
341 minline = line2;
342 }
343 if (maxl == 0 || maxl < l) {
344 maxl = l;
345 }
346 }
347 if (thresh < 0) {
348 G_ilist_add(list, minline);
349 }
350 else {
351 if (minl < thresh)
352 G_ilist_add(list, minline);
353 }
354 }
355 if (list->n_values == 0)
356 break;
357 nlines = 0;
358 for (i = 0; i < list->n_values; i++) {
359 line = list->value[i];
360 if (Vect_line_alive(&Out, abs(line))) {
361 if (next_dist(line, GV_RIGHT, morphof))
362 continue;
363 if (next_dist(line, GV_LEFT, morphof))
364 continue;
365 Vect_delete_line(&Out, abs(line));
366 nlines++;
367 }
368 }
369 if (nlines == 0)
370 break;
371 else
372 Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
373 }
374
375 if (thresh >= 0)
376 return 0;
377
378 for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
379 if (!Vect_node_alive(&Out, node))
380 continue;
381 nlines = Vect_get_node_n_lines(&Out, node);
382
383 if (nlines > 1)
384 continue;
385
386 line = Vect_get_node_line(&Out, node, 0);
387 if (line < 0)
388 Vect_get_line_nodes(&Out, -line, &n1, NULL);
389 else
390 Vect_get_line_nodes(&Out, line, NULL, &n1);
391
392 nlines = Vect_get_node_n_lines(&Out, n1);
393
394 if (nlines == 3) {
395 break_loop(line, n1, Points);
396 }
397 }
398
399 while (TRUE) {
400 G_verbose_message(_("Pass %d"), counter++);
401 Vect_reset_list(list);
402
403 for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
404 if (!Vect_node_alive(&Out, node))
405 continue;
406 nlines = Vect_get_node_n_lines(&Out, node);
407
408 if (nlines > 1)
409 continue;
410
411 line = Vect_get_node_line(&Out, node, 0);
412 if (line < 0)
413 Vect_get_line_nodes(&Out, -line, &n1, NULL);
414 else
415 Vect_get_line_nodes(&Out, line, NULL, &n1);
416
417 nlines = Vect_get_node_n_lines(&Out, n1);
418
419 if (nlines < 3)
420 continue;
421
422 Vect_read_line(&Out, Points, NULL, abs(line));
423 minline = line;
424 minl = Vect_line_length(Points);
425
426 if (nlines == 3) {
427 if (break_loop(line, n1, Points))
428 continue;
429 }
430
431 for (i = 0; i < nlines; i++) {
432 line2 = Vect_get_node_line(&Out, n1, i);
433 if (abs(line2) == abs(minline) || abs(line2) == abs(line))
434 continue;
435 if (line2 < 0)
436 Vect_get_line_nodes(&Out, -line2, &n2, NULL);
437 else
438 Vect_get_line_nodes(&Out, line2, NULL, &n2);
439 if (Vect_get_node_n_lines(&Out, n2) > 1)
440 continue;
441 Vect_read_line(&Out, Points, NULL, abs(line2));
442 l = Vect_line_length(Points);
443 if (minl > l) {
444 minl = l;
445 minline = line2;
446 }
447 }
448 if (thresh < 0 || minl < thresh)
449 G_ilist_add(list, minline);
450 }
451 if (list->n_values == 0)
452 break;
453 nlines = 0;
454 for (i = 0; i < list->n_values; i++) {
455 line = list->value[i];
456 if (Vect_line_alive(&Out, abs(line))) {
457 if (next_dist(line, GV_RIGHT, morphof))
458 continue;
459 if (next_dist(line, GV_LEFT, morphof))
460 continue;
461 Vect_delete_line(&Out, abs(line));
462 nlines++;
463 }
464 }
465 if (nlines == 0)
466 break;
467 else
468 Vect_merge_lines(&Out, GV_LINE, NULL, NULL);
469 }
470
471 return 0;
472 }
473
tie_up(void)474 int tie_up(void)
475 {
476 int i;
477 int node;
478 int nlines;
479 double xmin, ymin, x, y;
480 double dx, dy, dist, distmin;
481 struct line_pnts *Points;
482 struct line_pnts **IPoints;
483 struct line_cats *Cats;
484 int isl_allocated;
485 int area, isle, n_isles;
486 int ntied = 0;
487
488 Points = Vect_new_line_struct();
489 isl_allocated = 10;
490 IPoints = G_malloc(isl_allocated * sizeof(struct line_pnts *));
491 for (i = 0; i < isl_allocated; i++)
492 IPoints[i] = Vect_new_line_struct();
493 Cats = Vect_new_cats_struct();
494
495 for (node = 1; node <= Vect_get_num_nodes(&Out); node++) {
496 if (!Vect_node_alive(&Out, node))
497 continue;
498 nlines = Vect_get_node_n_lines(&Out, node);
499
500 if (nlines > 1)
501 continue;
502
503 Vect_get_node_coor(&Out, node, &x, &y, NULL);
504
505 /* find area for this node */
506 area = Vect_find_area(&In, x, y);
507
508 if (area == 0)
509 G_fatal_error(_("Node is outside any input area"));
510
511 /* get area outer ring */
512 Vect_get_area_points(&In, area, Points);
513
514 /* get area inner rings */
515 n_isles = Vect_get_area_num_isles(&In, area);
516 if (n_isles > isl_allocated) {
517 IPoints = (struct line_pnts **)
518 G_realloc(IPoints, (1 + n_isles) * sizeof(struct line_pnts *));
519 for (i = isl_allocated; i < n_isles; i++)
520 IPoints[i] = Vect_new_line_struct();
521 isl_allocated = n_isles;
522 }
523 for (i = 0; i < n_isles; i++) {
524 Vect_get_isle_points(&In, Vect_get_area_isle(&In, area, i),
525 IPoints[i]);
526 }
527
528 distmin = 1. / 0.; /* +inf */
529 xmin = x;
530 ymin = y;
531
532 /* find closest point to outer ring */
533 /* must be an existing vertex */
534 for (i = 0; i < Points->n_points - 1; i++) {
535 dx = x - Points->x[i];
536 dy = y - Points->y[i];
537
538 dist = dx * dx + dy * dy;
539
540 if (distmin > dist) {
541 distmin = dist;
542 xmin = Points->x[i];
543 ymin = Points->y[i];
544 }
545 }
546
547 /* find closest point to inner rings */
548 /* must be an existing vertex */
549 for (isle = 0; isle < n_isles; isle++) {
550 for (i = 0; i < IPoints[isle]->n_points - 1; i++) {
551 dx = x - IPoints[isle]->x[i];
552 dy = y - IPoints[isle]->y[i];
553
554 dist = dx * dx + dy * dy;
555
556 if (distmin > dist) {
557 distmin = dist;
558 xmin = IPoints[isle]->x[i];
559 ymin = IPoints[isle]->y[i];
560 }
561 }
562 }
563
564 if (xmin != x && ymin != y) {
565 Vect_get_area_cats(&In, area, Cats);
566 Vect_reset_line(Points);
567 Vect_append_point(Points, xmin, ymin, 0);
568 Vect_append_point(Points, x, y, 0);
569 Vect_write_line(&Out, GV_LINE, Points, Cats);
570 ntied++;
571 }
572 }
573
574 return ntied;
575 }
576