← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/vector/Vlib/break_polygons.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:
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
32
 
                                 *   this is stored because points are automaticaly marked as cross, even if not used 
 
58
                                 *   this is stored because points are automatically marked as cross, even if not used 
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++) {
 
159
        G_percent(i, nlines, 1);
94
160
        G_debug(3, "i =  %d", i);
95
161
        if (!Vect_line_alive(Map, i))
96
162
            continue;
130
196
 
131
197
            /* Already in DB? */
132
198
            fpoint = -1;
133
 
            RTreeSearch(RTree, &rect, (void *)srch, 0);
 
199
            RTreeSearch(RTree, &rect, (void *)srch, NULL);
134
200
            G_debug(3, "fpoint =  %d", fpoint);
135
201
 
136
202
            if (Points->n_points <= 2 ||
158
224
            }
159
225
 
160
226
            if (fpoint > 0) {   /* Found */
161
 
                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)
162
231
                    continue;   /* already marked */
163
232
 
164
233
                /* Check angles */
165
234
                if (cross) {
166
 
                    XPnts[fpoint].cross = 1;
 
235
                    XPnt.cross = 1;
167
236
                    nmarks++;
 
237
                    /* write point */
 
238
                    lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
 
239
                    write(xpntfd, &XPnt, sizeof(XPNT2));
168
240
                }
169
241
                else {
170
242
                    G_debug(3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1,
171
 
                            XPnts[fpoint].a1, a2, XPnts[fpoint].a2);
172
 
                    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 */
173
246
 
174
247
                    }
175
248
                    else {
176
 
                        XPnts[fpoint].cross = 1;
 
249
                        XPnt.cross = 1;
177
250
                        nmarks++;
 
251
                        /* write point */
 
252
                        lseek(xpntfd, (off_t) (fpoint - 1) * sizeof(XPNT2), SEEK_SET);
 
253
                        write(xpntfd, &XPnt, sizeof(XPNT2));
178
254
                    }
179
255
                }
180
256
            }
181
257
            else {
182
258
                /* Add to tree and to structure */
183
 
                RTreeInsertRect(&rect, npoints, &RTree, 0);
184
 
                if (npoints >= apoints) {
185
 
                    apoints += 10000;
186
 
                    XPnts =
187
 
                        (XPNT *) G_realloc(XPnts,
188
 
                                           (apoints + 1) * sizeof(XPNT));
189
 
                }
190
 
                XPnts[npoints].x = Points->x[j];
191
 
                XPnts[npoints].y = Points->y[j];
192
 
                XPnts[npoints].used = 0;
 
259
                RTreeInsertRect(&rect, npoints, RTree);
193
260
                if (j == 0 || j == (Points->n_points - 1) ||
194
261
                    Points->n_points < 3) {
195
 
                    XPnts[npoints].a1 = 0;
196
 
                    XPnts[npoints].a2 = 0;
197
 
                    XPnts[npoints].cross = 1;
 
262
                    XPnt.a1 = 0;
 
263
                    XPnt.a2 = 0;
 
264
                    XPnt.cross = 1;
198
265
                    nmarks++;
199
266
                }
200
267
                else {
201
 
                    XPnts[npoints].a1 = a1;
202
 
                    XPnts[npoints].a2 = a2;
203
 
                    XPnts[npoints].cross = 0;
 
268
                    XPnt.a1 = a1;
 
269
                    XPnt.a2 = a2;
 
270
                    XPnt.cross = 0;
204
271
                }
 
272
                /* write point */
 
273
                lseek(xpntfd, (off_t) (npoints - 1) * sizeof(XPNT2), SEEK_SET);
 
274
                write(xpntfd, &XPnt, sizeof(XPNT2));
205
275
 
206
276
                npoints++;
207
277
            }
208
278
        }
209
279
    }
210
280
 
211
 
    /* G_sleep (10); */
212
 
 
213
281
    nbreaks = 0;
214
282
 
215
283
    /* Second loop through lines (existing when loop is started, no need to process lines written again)
216
284
     * and break at points marked for break */
 
285
 
 
286
    G_message(_("Breaking polygons (pass 2: break at selected points)..."));
 
287
 
217
288
    for (i = 1; i <= nlines; i++) {
218
289
        int n_orig_points;
219
290
 
 
291
        G_percent(i, nlines, 1);
220
292
        G_debug(3, "i =  %d", i);
221
293
        if (!Vect_line_alive(Map, i))
222
294
            continue;
255
327
            RTreeSearch(RTree, &rect, (void *)srch, 0);
256
328
            G_debug(3, "fpoint =  %d", fpoint);
257
329
 
258
 
            if (XPnts[fpoint].cross) {  /* realy use to break line */
259
 
                XPnts[fpoint].used = 1;
260
 
            }
261
 
 
262
 
            /* break or write last segment of broken line */
263
 
            if ((j == (Points->n_points - 1) && broken) ||
264
 
                XPnts[fpoint].cross) {
265
 
                Vect_reset_line(BPoints);
266
 
                for (k = last; k <= j; k++) {
267
 
                    Vect_append_point(BPoints, Points->x[k], Points->y[k],
268
 
                                      Points->z[k]);
269
 
                }
270
 
 
271
 
                /* Result may collapse to one point */
272
 
                Vect_line_prune(BPoints);
273
 
                if (BPoints->n_points > 1) {
274
 
                    ret = Vect_write_line(Map, ltype, BPoints, Cats);
275
 
                    G_debug(3,
276
 
                            "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
277
 
                            ret, j, Points->n_points, BPoints->n_points);
278
 
                }
279
 
 
280
 
                if (!broken)
281
 
                    Vect_delete_line(Map, i);   /* not yet deleted */
282
 
 
283
 
                last = j;
284
 
                broken = 1;
285
 
                nbreaks++;
286
 
            }
287
 
        }
288
 
        if (!broken && n_orig_points > Points->n_points) {      /* was pruned before -> rewrite */
289
 
            if (Points->n_points > 1) {
290
 
                Vect_rewrite_line(Map, i, ltype, Points, Cats);
291
 
                G_debug(3, "Line %d pruned, npoints = %d", i,
292
 
                        Points->n_points);
293
 
            }
294
 
            else {
295
 
                Vect_delete_line(Map, i);
296
 
                G_debug(3, "Line %d was deleted", i);
297
 
            }
298
 
        }
299
 
        else {
300
 
            G_debug(3, "Line %d was not changed", i);
301
 
        }
302
 
    }
303
 
 
304
 
    /* Write points on breaks */
305
 
    if (Err) {
306
 
        Vect_reset_cats(Cats);
307
 
        for (i = 1; i < npoints; i++) {
308
 
            if (XPnts[i].used) {
309
 
                Vect_reset_line(Points);
310
 
                Vect_append_point(Points, XPnts[i].x, XPnts[i].y, 0);
311
 
                Vect_write_line(Err, GV_POINT, Points, Cats);
312
 
            }
313
 
        }
314
 
    }
315
 
 
316
 
    G_free(XPnts);
317
 
    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);
318
677
}