← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to vector/v.voronoi/skeleton.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/v.voronoi/skeleton.c
 
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
 
 
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
 
 
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
 
 
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
 
 
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 */
 
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
 
 
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
}