← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/vector/Vlib/break_polygons.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/vector/Vlib/break_polygons.c
1
1
/*!
2
 
   \file break_polygons.c
 
2
   \file lib/vector/Vlib/break_polygons.c
3
3
 
4
4
   \brief Vector library - clean geometry (break polygons)
5
5
 
6
6
   Higher level functions for reading/writing/manipulating vectors.
7
7
 
8
 
   (C) 2001-2008 by the GRASS Development Team
 
8
   (C) 2001-2009 by the GRASS Development Team
9
9
 
10
 
   This program is free software under the 
11
 
   GNU General Public License (>=v2). 
12
 
   Read the file COPYING that comes with GRASS
13
 
   for details.
 
10
   This program is free software under the GNU General Public License
 
11
   (>=v2).  Read the file COPYING that comes with GRASS for details.
14
12
 
15
13
   \author Radim Blazek
16
 
 
17
 
   \date 2001-2008
 
14
   \author Update for GRASS 7 Markus Metz
18
15
 */
19
16
 
20
17
#include <stdlib.h>
 
18
#include <sys/stat.h>
 
19
#include <fcntl.h>
 
20
#include <unistd.h>
21
21
#include <math.h>
22
 
#include <grass/gis.h>
23
 
#include <grass/Vect.h>
 
22
#include <grass/vector.h>
24
23
#include <grass/glocale.h>
25
24
 
 
25
/* TODO: 3D support
 
26
 *
 
27
 * atan2() gives angle from x-axis
 
28
 * this is unambiguous only in 2D, not in 3D
 
29
 *
 
30
 * one possibility would be to store unit vectors of length 1
 
31
 * in struct XPNT
 
32
 * double a1[3], a2[3];
 
33
 * 
 
34
 * length = sqrt(dx * dx + dy * dy + dz * dz);
 
35
 * dx /= length; dy /= length; dz /=length;
 
36
 * a1[0] = dx; a1[1] = dy; a1[2] = dz;
 
37
 *
 
38
 * get second dx, dy, dz
 
39
 * length = sqrt(dx * dx + dy * dy + dz * dz);
 
40
 * dx /= length; dy /= length; dz /=length;
 
41
 * a2[0] = dx; a2[1] = dy; a2[2] = dz;
 
42
 *
 
43
 * equal angles
 
44
 * if (a1[0] == a2[0] && a1[1] == a2[1] && a1[2] == a2[2])
 
45
 *
 
46
 * disadvantage: increased memory consumption
 
47
 *
 
48
 * new function Vect_break_faces() ?
 
49
 * 
 
50
 */
 
51
 
26
52
typedef struct
27
53
{
28
 
    double x, y;
 
54
    double x, y;                /* coords */
29
55
    double a1, a2;              /* angles */
30
56
    char cross;                 /* 0 - do not break, 1 - break */
31
57
    char used;                  /* 0 - was not used to break line, 1 - was used to break line
33
59
                                 *   later to break lines */
34
60
} XPNT;
35
61
 
 
62
typedef struct
 
63
{
 
64
    double a1, a2;              /* angles */
 
65
    char cross;                 /* 0 - do not break, 1 - break */
 
66
    char used;                  /* 0 - was not used to break line, 1 - was used to break line
 
67
                                 *   this is stored because points are automatically marked as cross, even if not used 
 
68
                                 *   later to break lines */
 
69
} XPNT2;
 
70
 
36
71
static int fpoint;
37
72
 
38
73
/* Function called from RTreeSearch for point found */
39
 
void srch(int id, int *arg)
 
74
static int srch(int id, const struct RTree_Rect *rect, int *arg)
40
75
{
41
76
    fpoint = id;
42
 
}
43
 
 
44
 
/*!
45
 
   \brief Break polygons in vector map.
46
 
 
47
 
   Breaks lines specified by type in vector map. Points at intersections may be optionally 
48
 
   written to error map. Input map must be opened on level 2 for update at least on GV_BUILD_BASE.
49
 
 
50
 
   Function is optimized for closed polygons rigs (e.g. imported from OGR) but with clean geometry -
51
 
   adjacent polygons mostly have identical boundary. Function creates database of ALL points
52
 
   in the map, and then is looking for those where polygons should be broken. 
53
 
   Lines may be broken only at points existing in input map! 
54
 
 
55
 
   \param Map input map where polygons will be broken
56
 
   \param type type of line to be broken
57
 
   \param Err vector map where points at intersections will be written or NULL
58
 
 
59
 
   \return
60
 
 */
61
 
void
62
 
Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
 
77
    
 
78
    return 0;   /* stop searching */
 
79
}
 
80
 
 
81
/* function used by binary tree to compare items */
 
82
static int compare_xpnts(const void *Xpnta, const void *Xpntb)
 
83
{
 
84
    XPNT *a, *b;
 
85
 
 
86
    a = (XPNT *)Xpnta;
 
87
    b = (XPNT *)Xpntb;
 
88
 
 
89
    if (a->x > b->x)
 
90
        return 1;
 
91
    else if (a->x < b->x)
 
92
        return -1;
 
93
    else {
 
94
        if (a->y > b->y)
 
95
            return 1;
 
96
        else if (a->y < b->y)
 
97
            return -1;
 
98
        else
 
99
            return 0;
 
100
    }
 
101
 
 
102
    G_warning(_("Break polygons: Bug in binary tree!"));
 
103
    return 1;
 
104
}
 
105
 
 
106
/* break polygons using a file-based search index */
 
107
void Vect_break_polygons_file(struct Map_info *Map, int type, struct Map_info *Err)
63
108
{
64
109
    struct line_pnts *BPoints, *Points;
65
 
    struct line_cats *Cats;
 
110
    struct line_cats *Cats, *ErrCats;
66
111
    int i, j, k, ret, ltype, broken, last, nlines;
67
112
    int nbreaks;
68
 
    struct Node *RTree;
69
 
    int apoints, npoints, nallpoints, nmarks;
70
 
    XPNT *XPnts;
71
 
    struct Rect rect;
 
113
    struct RTree *RTree;
 
114
    int npoints, nallpoints, nmarks;
 
115
    XPNT2 XPnt;
72
116
    double dx, dy, a1 = 0, a2 = 0;
73
 
    int closed, last_point, cross;
74
 
 
75
 
    RTree = RTreeNewIndex();
 
117
    int closed, last_point;
 
118
    char cross;
 
119
    int fd, xpntfd;
 
120
    char *filename;
 
121
    static struct RTree_Rect rect;
 
122
    static int rect_init = 0;
 
123
 
 
124
    if (!rect_init) {
 
125
        rect.boundary = G_malloc(6 * sizeof(RectReal));
 
126
        rect_init = 6;
 
127
    }
 
128
    
 
129
    G_debug(1, "File-based version of Vect_break_polygons()");
 
130
 
 
131
    filename = G_tempfile();
 
132
    fd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
 
133
    RTree = RTreeCreateTree(fd, 0, 2);
 
134
    remove(filename);
 
135
 
 
136
    filename = G_tempfile();
 
137
    xpntfd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
 
138
    remove(filename);
76
139
 
77
140
    BPoints = Vect_new_line_struct();
78
141
    Points = Vect_new_line_struct();
79
142
    Cats = Vect_new_cats_struct();
 
143
    ErrCats = Vect_new_cats_struct();
80
144
 
81
145
    nlines = Vect_get_num_lines(Map);
82
146
 
84
148
    /* Go through all lines in vector, and add each point to structure of points,
85
149
     * if such point already exists check angles of segments and if differ mark for break */
86
150
 
87
 
    apoints = 0;
88
151
    nmarks = 0;
89
152
    npoints = 1;                /* index starts from 1 ! */
90
153
    nallpoints = 0;
91
 
    XPnts = NULL;
 
154
    XPnt.used = 0;
 
155
 
 
156
    G_message(_("Breaking polygons (pass 1: select break points)..."));
92
157
 
93
158
    for (i = 1; i <= nlines; i++) {
94
159
        G_percent(i, nlines, 1);
131
196
 
132
197
            /* Already in DB? */
133
198
            fpoint = -1;
134
 
            RTreeSearch(RTree, &rect, (void *)srch, 0);
 
199
            RTreeSearch(RTree, &rect, (void *)srch, NULL);
135
200
            G_debug(3, "fpoint =  %d", fpoint);
136
201
 
137
202
            if (Points->n_points <= 2 ||
159
224
            }
160
225
 
161
226
            if (fpoint > 0) {   /* Found */
162
 
                if (XPnts[fpoint].cross == 1)
 
227
                /* read point */
 
228
                lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
 
229
                read(xpntfd, &XPnt, sizeof(XPNT2));
 
230
                if (XPnt.cross == 1)
163
231
                    continue;   /* already marked */
164
232
 
165
233
                /* Check angles */
166
234
                if (cross) {
167
 
                    XPnts[fpoint].cross = 1;
 
235
                    XPnt.cross = 1;
168
236
                    nmarks++;
 
237
                    /* write point */
 
238
                    lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
 
239
                    write(xpntfd, &XPnt, sizeof(XPNT2));
169
240
                }
170
241
                else {
171
242
                    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 */
 
243
                            XPnt.a1, a2, XPnt.a2);
 
244
                    if ((a1 == XPnt.a1 && a2 == XPnt.a2) ||
 
245
                        (a1 == XPnt.a2 && a2 == XPnt.a1)) {     /* identical */
174
246
 
175
247
                    }
176
248
                    else {
177
 
                        XPnts[fpoint].cross = 1;
 
249
                        XPnt.cross = 1;
178
250
                        nmarks++;
 
251
                        /* write point */
 
252
                        lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
 
253
                        write(xpntfd, &XPnt, sizeof(XPNT2));
179
254
                    }
180
255
                }
181
256
            }
182
257
            else {
183
258
                /* 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;
 
259
                RTreeInsertRect(&rect, npoints, RTree);
194
260
                if (j == 0 || j == (Points->n_points - 1) ||
195
261
                    Points->n_points < 3) {
196
 
                    XPnts[npoints].a1 = 0;
197
 
                    XPnts[npoints].a2 = 0;
198
 
                    XPnts[npoints].cross = 1;
 
262
                    XPnt.a1 = 0;
 
263
                    XPnt.a2 = 0;
 
264
                    XPnt.cross = 1;
199
265
                    nmarks++;
200
266
                }
201
267
                else {
202
 
                    XPnts[npoints].a1 = a1;
203
 
                    XPnts[npoints].a2 = a2;
204
 
                    XPnts[npoints].cross = 0;
 
268
                    XPnt.a1 = a1;
 
269
                    XPnt.a2 = a2;
 
270
                    XPnt.cross = 0;
205
271
                }
 
272
                /* write point */
 
273
                lseek(xpntfd, (off_t) (npoints - 1) * sizeof(XPNT2), SEEK_SET);
 
274
                write(xpntfd, &XPnt, sizeof(XPNT2));
206
275
 
207
276
                npoints++;
208
277
            }
209
278
        }
210
279
    }
211
280
 
212
 
    /* G_sleep (10); */
213
 
 
214
281
    nbreaks = 0;
215
282
 
216
283
    /* Second loop through lines (existing when loop is started, no need to process lines written again)
217
284
     * and break at points marked for break */
 
285
 
 
286
    G_message(_("Breaking polygons (pass 2: break at selected points)..."));
 
287
 
218
288
    for (i = 1; i <= nlines; i++) {
219
289
        int n_orig_points;
220
290
 
257
327
            RTreeSearch(RTree, &rect, (void *)srch, 0);
258
328
            G_debug(3, "fpoint =  %d", fpoint);
259
329
 
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);
 
330
            /* read point */
 
331
            lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
 
332
            read(xpntfd, &XPnt, sizeof(XPNT2));
 
333
 
 
334
            /* break or write last segment of broken line */
 
335
            if ((j == (Points->n_points - 1) && broken) ||
 
336
                XPnt.cross) {
 
337
                Vect_reset_line(BPoints);
 
338
                for (k = last; k <= j; k++) {
 
339
                    Vect_append_point(BPoints, Points->x[k], Points->y[k],
 
340
                                      Points->z[k]);
 
341
                }
 
342
 
 
343
                /* Result may collapse to one point */
 
344
                Vect_line_prune(BPoints);
 
345
                if (BPoints->n_points > 1) {
 
346
                    ret = Vect_write_line(Map, ltype, BPoints, Cats);
 
347
                    G_debug(3,
 
348
                            "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
 
349
                            ret, j, Points->n_points, BPoints->n_points);
 
350
                }
 
351
 
 
352
                if (!broken)
 
353
                    Vect_delete_line(Map, i);   /* not yet deleted */
 
354
 
 
355
                /* Write points on breaks */
 
356
                if (Err) {
 
357
                    if (XPnt.cross && !XPnt.used) {
 
358
                        Vect_reset_line(BPoints);
 
359
                        Vect_append_point(BPoints, Points->x[j], Points->y[j], 0);
 
360
                        Vect_write_line(Err, GV_POINT, BPoints, ErrCats);
 
361
                    }
 
362
                    if (!XPnt.used) {
 
363
                        XPnt.used = 1;
 
364
                        /* write point */
 
365
                        lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
 
366
                        write(xpntfd, &XPnt, sizeof(XPNT2));
 
367
                    }
 
368
                }
 
369
 
 
370
                last = j;
 
371
                broken = 1;
 
372
                nbreaks++;
 
373
            }
 
374
        }
 
375
        if (!broken && n_orig_points > Points->n_points) {      /* was pruned before -> rewrite */
 
376
            if (Points->n_points > 1) {
 
377
                Vect_rewrite_line(Map, i, ltype, Points, Cats);
 
378
                G_debug(3, "Line %d pruned, npoints = %d", i,
 
379
                        Points->n_points);
 
380
            }
 
381
            else {
 
382
                Vect_delete_line(Map, i);
 
383
                G_debug(3, "Line %d was deleted", i);
 
384
            }
 
385
        }
 
386
        else {
 
387
            G_debug(3, "Line %d was not changed", i);
 
388
        }
 
389
    }
 
390
 
 
391
    close(RTree->fd);
 
392
    RTreeDestroyTree(RTree);
 
393
    close(xpntfd);
 
394
    Vect_destroy_line_struct(Points);
 
395
    Vect_destroy_line_struct(BPoints);
 
396
    Vect_destroy_cats_struct(Cats);
 
397
    Vect_destroy_cats_struct(ErrCats);
 
398
    G_verbose_message(_("Breaks: %d"), nbreaks);
 
399
}
 
400
 
 
401
 
 
402
/* break polygons using a memory-based search index */
 
403
void Vect_break_polygons_mem(struct Map_info *Map, int type, struct Map_info *Err)
 
404
{
 
405
    struct line_pnts *BPoints, *Points;
 
406
    struct line_cats *Cats, *ErrCats;
 
407
    int i, j, k, ret, ltype, broken, last, nlines;
 
408
    int nbreaks;
 
409
    struct RB_TREE *RBTree;
 
410
    int npoints, nallpoints, nmarks;
 
411
    XPNT *XPnt_found, XPnt_search;
 
412
    double dx, dy, a1 = 0, a2 = 0;
 
413
    int closed, last_point, cross;
 
414
 
 
415
    G_debug(1, "Memory-based version of Vect_break_polygons()");
 
416
 
 
417
    RBTree = rbtree_create(compare_xpnts, sizeof(XPNT));
 
418
 
 
419
    BPoints = Vect_new_line_struct();
 
420
    Points = Vect_new_line_struct();
 
421
    Cats = Vect_new_cats_struct();
 
422
    ErrCats = Vect_new_cats_struct();
 
423
 
 
424
    nlines = Vect_get_num_lines(Map);
 
425
 
 
426
    G_debug(3, "nlines =  %d", nlines);
 
427
    /* Go through all lines in vector, and add each point to structure of points,
 
428
     * if such point already exists check angles of segments and if differ mark for break */
 
429
 
 
430
    nmarks = 0;
 
431
    npoints = 0;
 
432
    nallpoints = 0;
 
433
    XPnt_search.used = 0;
 
434
 
 
435
    G_message(_("Breaking polygons (pass 1: select break points)..."));
 
436
 
 
437
    for (i = 1; i <= nlines; i++) {
 
438
        G_percent(i, nlines, 1);
 
439
        G_debug(3, "i =  %d", i);
 
440
        if (!Vect_line_alive(Map, i))
 
441
            continue;
 
442
 
 
443
        ltype = Vect_read_line(Map, Points, Cats, i);
 
444
        if (!(ltype & type))
 
445
            continue;
 
446
 
 
447
        /* This would be confused by duplicate coordinates (angle cannot be calculated) ->
 
448
         * prune line first */
 
449
        Vect_line_prune(Points);
 
450
 
 
451
        /* If first and last point are identical it is close polygon, we dont need to register last point
 
452
         * and we can calculate angle for first.
 
453
         * If first and last point are not identical we have to mark for break both */
 
454
        last_point = Points->n_points - 1;
 
455
        if (Points->x[0] == Points->x[last_point] &&
 
456
            Points->y[0] == Points->y[last_point])
 
457
            closed = 1;
 
458
        else
 
459
            closed = 0;
 
460
 
 
461
        for (j = 0; j < Points->n_points; j++) {
 
462
            G_debug(3, "j =  %d", j);
 
463
            nallpoints++;
 
464
 
 
465
            if (j == last_point && closed)
 
466
                continue;       /* do not register last of close polygon */
 
467
 
 
468
            XPnt_search.x = Points->x[j];
 
469
            XPnt_search.y = Points->y[j];
 
470
 
 
471
            /* Already in DB? */
 
472
            XPnt_found = rbtree_find(RBTree, &XPnt_search);
 
473
 
 
474
            if (Points->n_points <= 2 ||
 
475
                (!closed && (j == 0 || j == last_point))) {
 
476
                cross = 1;      /* mark for cross in any case */
 
477
            }
 
478
            else {              /* calculate angles */
 
479
                cross = 0;
 
480
                if (j == 0 && closed) { /* closed polygon */
 
481
                    dx = Points->x[last_point] - Points->x[0];
 
482
                    dy = Points->y[last_point] - Points->y[0];
 
483
                    a1 = atan2(dy, dx);
 
484
                    dx = Points->x[1] - Points->x[0];
 
485
                    dy = Points->y[1] - Points->y[0];
 
486
                    a2 = atan2(dy, dx);
 
487
                }
 
488
                else {
 
489
                    dx = Points->x[j - 1] - Points->x[j];
 
490
                    dy = Points->y[j - 1] - Points->y[j];
 
491
                    a1 = atan2(dy, dx);
 
492
                    dx = Points->x[j + 1] - Points->x[j];
 
493
                    dy = Points->y[j + 1] - Points->y[j];
 
494
                    a2 = atan2(dy, dx);
 
495
                }
 
496
            }
 
497
 
 
498
            if (XPnt_found) {   /* found */
 
499
                if (XPnt_found->cross == 1)
 
500
                    continue;   /* already marked */
 
501
 
 
502
                /* check angles */
 
503
                if (cross) {
 
504
                    XPnt_found->cross = 1;
 
505
                    nmarks++;
 
506
                }
 
507
                else {
 
508
                    G_debug(3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1,
 
509
                            XPnt_found->a1, a2, XPnt_found->a2);
 
510
                    if ((a1 == XPnt_found->a1 && a2 == XPnt_found->a2) ||
 
511
                        (a1 == XPnt_found->a2 && a2 == XPnt_found->a1)) {       /* identical */
 
512
 
 
513
                    }
 
514
                    else {
 
515
                        XPnt_found->cross = 1;
 
516
                        nmarks++;
 
517
                    }
 
518
                }
 
519
            }
 
520
            else {
 
521
                if (j == 0 || j == (Points->n_points - 1) ||
 
522
                    Points->n_points < 3) {
 
523
                    XPnt_search.a1 = 0;
 
524
                    XPnt_search.a2 = 0;
 
525
                    XPnt_search.cross = 1;
 
526
                    nmarks++;
 
527
                }
 
528
                else {
 
529
                    XPnt_search.a1 = a1;
 
530
                    XPnt_search.a2 = a2;
 
531
                    XPnt_search.cross = 0;
 
532
                }
 
533
 
 
534
                /* Add to tree */
 
535
                rbtree_insert(RBTree, &XPnt_search);
 
536
                npoints++;
 
537
            }
 
538
        }
 
539
    }
 
540
 
 
541
    nbreaks = 0;
 
542
    nallpoints = 0;
 
543
    G_debug(2, "Break polygons: unique vertices: %ld", (long int)RBTree->count);
 
544
 
 
545
    /* uncomment to check if search tree is healthy */
 
546
    /* if (rbtree_debug(RBTree, RBTree->root) == 0)
 
547
        G_warning("Break polygons: RBTree not ok"); */
 
548
 
 
549
    /* Second loop through lines (existing when loop is started, no need to process lines written again)
 
550
     * and break at points marked for break */
 
551
 
 
552
    G_message(_("Breaking polygons (pass 2: break at selected points)..."));
 
553
 
 
554
    for (i = 1; i <= nlines; i++) {
 
555
        int n_orig_points;
 
556
 
 
557
        G_percent(i, nlines, 1);
 
558
        G_debug(3, "i =  %d", i);
 
559
        if (!Vect_line_alive(Map, i))
 
560
            continue;
 
561
 
 
562
        ltype = Vect_read_line(Map, Points, Cats, i);
 
563
        if (!(ltype & type))
 
564
            continue;
 
565
        if (!(ltype & GV_LINES))
 
566
            continue;           /* Nonsense to break points */
 
567
 
 
568
        /* Duplicates would result in zero length lines -> prune line first */
 
569
        n_orig_points = Points->n_points;
 
570
        Vect_line_prune(Points);
 
571
 
 
572
        broken = 0;
 
573
        last = 0;
 
574
        G_debug(3, "n_points =  %d", Points->n_points);
 
575
        for (j = 1; j < Points->n_points; j++) {
 
576
            G_debug(3, "j =  %d", j);
 
577
            nallpoints++;
 
578
 
 
579
            if (Points->n_points <= 1 ||
 
580
                (j == (Points->n_points - 1) && !broken))
 
581
                break;
 
582
            /* One point only or 
 
583
             * last point and line is not broken, do nothing */
 
584
 
 
585
            XPnt_search.x = Points->x[j];
 
586
            XPnt_search.y = Points->y[j];
 
587
 
 
588
            XPnt_found = rbtree_find(RBTree, &XPnt_search);
 
589
 
 
590
            /* all points must be in the search tree, without duplicates */
 
591
            if (XPnt_found == NULL)
 
592
                G_fatal_error(_("Point not in search tree!"));
 
593
 
 
594
            /* break or write last segment of broken line */
 
595
            if ((j == (Points->n_points - 1) && broken) ||
 
596
                XPnt_found->cross) {
 
597
                Vect_reset_line(BPoints);
 
598
                for (k = last; k <= j; k++) {
 
599
                    Vect_append_point(BPoints, Points->x[k], Points->y[k],
 
600
                                      Points->z[k]);
 
601
                }
 
602
 
 
603
                /* Result may collapse to one point */
 
604
                Vect_line_prune(BPoints);
 
605
                if (BPoints->n_points > 1) {
 
606
                    ret = Vect_write_line(Map, ltype, BPoints, Cats);
 
607
                    G_debug(3,
 
608
                            "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
 
609
                            ret, j, Points->n_points, BPoints->n_points);
 
610
                }
 
611
 
 
612
                if (!broken)
 
613
                    Vect_delete_line(Map, i);   /* not yet deleted */
 
614
 
 
615
                /* Write points on breaks */
 
616
                if (Err) {
 
617
                    if (XPnt_found->cross && !XPnt_found->used) {
 
618
                        Vect_reset_line(BPoints);
 
619
                        Vect_append_point(BPoints, Points->x[j], Points->y[j], 0);
 
620
                        Vect_write_line(Err, GV_POINT, BPoints, ErrCats);
 
621
                    }
 
622
                    XPnt_found->used = 1;
 
623
                }
 
624
 
 
625
                last = j;
 
626
                broken = 1;
 
627
                nbreaks++;
 
628
            }
 
629
        }
 
630
        if (!broken && n_orig_points > Points->n_points) {      /* was pruned before -> rewrite */
 
631
            if (Points->n_points > 1) {
 
632
                Vect_rewrite_line(Map, i, ltype, Points, Cats);
 
633
                G_debug(3, "Line %d pruned, npoints = %d", i,
 
634
                        Points->n_points);
 
635
            }
 
636
            else {
 
637
                Vect_delete_line(Map, i);
 
638
                G_debug(3, "Line %d was deleted", i);
 
639
            }
 
640
        }
 
641
        else {
 
642
            G_debug(3, "Line %d was not changed", i);
 
643
        }
 
644
    }
 
645
 
 
646
    rbtree_destroy(RBTree);
 
647
    Vect_destroy_line_struct(Points);
 
648
    Vect_destroy_line_struct(BPoints);
 
649
    Vect_destroy_cats_struct(Cats);
 
650
    G_verbose_message(_("Breaks: %d"), nbreaks);
 
651
}
 
652
 
 
653
/*!
 
654
   \brief Break polygons in vector map
 
655
 
 
656
   Breaks lines specified by type in vector map. Points at
 
657
   intersections may be optionally written to error map. Input vector
 
658
   map must be opened on level 2 for update at least on GV_BUILD_BASE.
 
659
 
 
660
   Function is optimized for closed polygons rings (e.g. imported from
 
661
   OGR) but with clean geometry - adjacent polygons mostly have
 
662
   identical boundary. Function creates database of ALL points in the
 
663
   vector map, and then is looking for those where polygons should be
 
664
   broken.  Lines may be broken only at points existing in input
 
665
   vector map!
 
666
 
 
667
   \param Map input map where polygons will be broken
 
668
   \param type type of line to be broken (GV_LINE or GV_BOUNDARY)
 
669
   \param Err vector map where points at intersections will be written or NULL
 
670
 */
 
671
void Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
 
672
{
 
673
    if (getenv("GRASS_VECTOR_LOWMEM"))
 
674
        return Vect_break_polygons_file(Map, type, Err);
 
675
    else
 
676
        return Vect_break_polygons_mem(Map, type, Err);
320
677
}