← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/vector/Vlib/select.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/select.c
1
1
/*!
2
 
   \file select.c
3
 
 
4
 
   \brief Vector library - select vector features
5
 
 
6
 
   Higher level functions for reading/writing/manipulating vectors.
7
 
 
8
 
   (C) 2001-2008 by the GRASS Development Team
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.
 
2
   \file lib/vector/Vlib/select.c
 
3
 
 
4
   \brief Vector library - spatial index
 
5
 
 
6
   Higher level functions for a custom spatial index.
 
7
 
 
8
   (C) 2001-2009 by the GRASS Development Team
 
9
 
 
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
18
14
 */
19
15
 
20
16
#include <stdlib.h>
21
 
#include <grass/gis.h>
22
 
#include <grass/Vect.h>
23
 
 
24
 
/*!
25
 
   \brief Select lines by box.
26
 
 
27
 
   Select lines whose boxes overlap specified box!!!
28
 
   It means that selected line may or may not overlap the box.
29
 
 
30
 
   \param Map vector map
31
 
   \param Box bounding box
32
 
   \param type line type
33
 
   \param[out] list output list, must be initialized
34
 
 
35
 
   \return number of lines
36
 
 */
37
 
int
38
 
Vect_select_lines_by_box(struct Map_info *Map, BOUND_BOX * Box,
39
 
                         int type, struct ilist *list)
40
 
{
41
 
    int i, line, nlines;
42
 
    struct Plus_head *plus;
43
 
    P_LINE *Line;
44
 
    static struct ilist *LocList = NULL;
45
 
 
46
 
    G_debug(3, "Vect_select_lines_by_box()");
47
 
    G_debug(3, "  Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
48
 
            Box->E, Box->W, Box->T, Box->B);
49
 
    plus = &(Map->plus);
50
 
 
51
 
    if (!(plus->Spidx_built)) {
52
 
        G_debug(3, "Building spatial index.");
53
 
        Vect_build_sidx_from_topo(Map);
54
 
    }
55
 
 
56
 
    list->n_values = 0;
57
 
    if (!LocList)
58
 
        LocList = Vect_new_list();
59
 
 
60
 
    nlines = dig_select_lines(plus, Box, LocList);
61
 
    G_debug(3, "  %d lines selected (all types)", nlines);
62
 
 
63
 
    /* Remove lines of not requested types */
64
 
    for (i = 0; i < nlines; i++) {
65
 
        line = LocList->value[i];
66
 
        if (plus->Line[line] == NULL)
67
 
            continue;           /* Should not happen */
68
 
        Line = plus->Line[line];
69
 
        if (!(Line->type & type))
70
 
            continue;
71
 
        dig_list_add(list, line);
72
 
    }
73
 
 
74
 
    G_debug(3, "  %d lines of requested type", list->n_values);
75
 
 
76
 
    return list->n_values;
77
 
}
78
 
 
79
 
/*!
80
 
   \brief Select areas by box.
81
 
 
82
 
   Select areas whose boxes overlap specified box!!!
83
 
   It means that selected area may or may not overlap the box.
84
 
 
85
 
   \param Map vector map
86
 
   \param Box bounding box
87
 
   \param[out] output list, must be initialized
88
 
 
89
 
   \return number of areas
90
 
 */
91
 
int
92
 
Vect_select_areas_by_box(struct Map_info *Map, BOUND_BOX * Box,
93
 
                         struct ilist *list)
94
 
{
95
 
    int i;
96
 
 
97
 
    G_debug(3, "Vect_select_areas_by_box()");
98
 
    G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
99
 
            Box->E, Box->W, Box->T, Box->B);
100
 
 
101
 
    if (!(Map->plus.Spidx_built)) {
102
 
        G_debug(3, "Building spatial index.");
103
 
        Vect_build_sidx_from_topo(Map);
104
 
    }
105
 
 
106
 
    dig_select_areas(&(Map->plus), Box, list);
107
 
    G_debug(3, "  %d areas selected", list->n_values);
108
 
    for (i = 0; i < list->n_values; i++) {
109
 
        G_debug(3, "  area = %d pointer to area structure = %lx",
110
 
                list->value[i],
111
 
                (unsigned long)Map->plus.Area[list->value[i]]);
112
 
 
113
 
    }
114
 
    return list->n_values;
115
 
}
116
 
 
117
 
 
118
 
/*!
119
 
   \brief Select isles by box.
120
 
 
121
 
   Select isles whose boxes overlap specified box!!!
122
 
   It means that selected isle may or may not overlap the box.
123
 
 
124
 
   \param Map vector map
125
 
   \param Box bounding box
126
 
   \param[out] list output list, must be initialized
127
 
 
128
 
   \return number of isles
129
 
 */
130
 
int
131
 
Vect_select_isles_by_box(struct Map_info *Map, BOUND_BOX * Box,
132
 
                         struct ilist *list)
133
 
{
134
 
    G_debug(3, "Vect_select_isles_by_box()");
135
 
    G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
136
 
            Box->E, Box->W, Box->T, Box->B);
137
 
 
138
 
    if (!(Map->plus.Spidx_built)) {
139
 
        G_debug(3, "Building spatial index.");
140
 
        Vect_build_sidx_from_topo(Map);
141
 
    }
142
 
 
143
 
    dig_select_isles(&(Map->plus), Box, list);
144
 
    G_debug(3, "  %d isles selected", list->n_values);
145
 
 
146
 
    return list->n_values;
147
 
}
148
 
 
149
 
/*!
150
 
   \brief Select nodes by box.
151
 
 
152
 
   \param Map vector map
153
 
   \param Box bounding box
154
 
   \param[out] list output list, must be initialized
155
 
 
156
 
   \return number of nodes
157
 
 */
158
 
int
159
 
Vect_select_nodes_by_box(struct Map_info *Map, BOUND_BOX * Box,
160
 
                         struct ilist *list)
161
 
{
162
 
    struct Plus_head *plus;
163
 
 
164
 
    G_debug(3, "Vect_select_nodes_by_box()");
165
 
    G_debug(3, "Box(N,S,E,W,T,B): %e, %e, %e, %e, %e, %e", Box->N, Box->S,
166
 
            Box->E, Box->W, Box->T, Box->B);
167
 
 
168
 
    plus = &(Map->plus);
169
 
 
170
 
    if (!(plus->Spidx_built)) {
171
 
        G_debug(3, "Building spatial index.");
172
 
        Vect_build_sidx_from_topo(Map);
173
 
    }
174
 
 
175
 
    list->n_values = 0;
176
 
 
177
 
    dig_select_nodes(plus, Box, list);
178
 
    G_debug(3, "  %d nodes selected", list->n_values);
179
 
 
180
 
    return list->n_values;
181
 
}
182
 
 
183
 
/*!
184
 
   \brief Select lines by Polygon with optional isles. 
185
 
 
186
 
   Polygons should be closed, i.e. first and last points must be identical.
187
 
 
188
 
   \param Map vector map
189
 
   \param Polygon outer ring
190
 
   \param nisles number of islands or 0
191
 
   \param Isles array of islands or NULL
192
 
   \param type line type
193
 
   \param[out] list output list, must be initialised
194
 
 
195
 
   \return number of lines
196
 
 */
197
 
int
198
 
Vect_select_lines_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
199
 
                             int nisles, struct line_pnts **Isles, int type,
200
 
                             struct ilist *List)
201
 
{
202
 
    int i;
203
 
    BOUND_BOX box;
204
 
    static struct line_pnts *LPoints = NULL;
205
 
    static struct ilist *LocList = NULL;
206
 
 
207
 
    /* TODO: this function was not tested with isles */
208
 
    G_debug(3, "Vect_select_lines_by_polygon() nisles = %d", nisles);
209
 
 
210
 
    List->n_values = 0;
211
 
    if (!LPoints)
212
 
        LPoints = Vect_new_line_struct();
213
 
    if (!LocList)
214
 
        LocList = Vect_new_list();
215
 
 
216
 
    /* Select first all lines by box */
217
 
    dig_line_box(Polygon, &box);
218
 
    Vect_select_lines_by_box(Map, &box, type, LocList);
219
 
    G_debug(3, "  %d lines selected by box", LocList->n_values);
220
 
 
221
 
    /* Check all lines if intersect the polygon */
222
 
    for (i = 0; i < LocList->n_values; i++) {
223
 
        int j, line, intersect = 0;
224
 
 
225
 
        line = LocList->value[i];
226
 
        /* Read line points */
227
 
        Vect_read_line(Map, LPoints, NULL, line);
228
 
 
229
 
        /* Check if any of line vertices is within polygon */
230
 
        for (j = 0; j < LPoints->n_points; j++) {
231
 
            if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Polygon) >= 1) {       /* inside polygon */
232
 
                int k, inisle = 0;
233
 
 
234
 
                for (k = 0; k < nisles; k++) {
235
 
                    if (Vect_point_in_poly(LPoints->x[j], LPoints->y[j], Isles[k]) >= 1) {      /* in isle */
236
 
                        inisle = 1;
237
 
                        break;
238
 
                    }
239
 
                }
240
 
 
241
 
                if (!inisle) {  /* inside polygon, outside isles -> select */
242
 
                    intersect = 1;
243
 
                    break;
244
 
                }
245
 
            }
246
 
        }
247
 
        if (intersect) {
248
 
            dig_list_add(List, line);
249
 
            continue;
250
 
        }
251
 
 
252
 
        /* Check intersections of the line with area/isles boundary */
253
 
        /* Outer boundary */
254
 
        if (Vect_line_check_intersection(LPoints, Polygon, 0)) {
255
 
            dig_list_add(List, line);
256
 
            continue;
257
 
        }
258
 
 
259
 
        /* Islands */
260
 
        for (j = 0; j < nisles; j++) {
261
 
            if (Vect_line_check_intersection(LPoints, Isles[j], 0)) {
262
 
                intersect = 1;
263
 
                break;
264
 
            }
265
 
        }
266
 
        if (intersect) {
267
 
            dig_list_add(List, line);
268
 
        }
269
 
    }
270
 
 
271
 
    G_debug(4, "  %d lines selected by polygon", List->n_values);
272
 
 
273
 
    return List->n_values;
274
 
}
275
 
 
276
 
 
277
 
/*!
278
 
   \brief Select areas by Polygon with optional isles. 
279
 
 
280
 
   Polygons should be closed, i.e. first and last points must be identical.
281
 
 
282
 
   Warning : values in list may be duplicate!
283
 
 
284
 
 
285
 
   \param Map vector map
286
 
   \param Polygon outer ring
287
 
   \param nisles number of islands or 0
288
 
   \param Isles array of islands or NULL
289
 
   \param[out] list output list, must be initialised
290
 
 
291
 
   \return number of areas
292
 
 */
293
 
int
294
 
Vect_select_areas_by_polygon(struct Map_info *Map, struct line_pnts *Polygon,
295
 
                             int nisles, struct line_pnts **Isles,
296
 
                             struct ilist *List)
297
 
{
298
 
    int i, area;
299
 
    static struct ilist *BoundList = NULL;
300
 
 
301
 
    /* TODO: this function was not tested with isles */
302
 
    G_debug(3, "Vect_select_areas_by_polygon() nisles = %d", nisles);
303
 
 
304
 
    List->n_values = 0;
305
 
    if (!BoundList)
306
 
        BoundList = Vect_new_list();
307
 
 
308
 
    /* Select boundaries by polygon */
309
 
    Vect_select_lines_by_polygon(Map, Polygon, nisles, Isles, GV_BOUNDARY,
310
 
                                 BoundList);
311
 
 
312
 
    /* Add areas on left/right side of selected boundaries */
313
 
    for (i = 0; i < BoundList->n_values; i++) {
314
 
        int line, left, right;
315
 
 
316
 
        line = BoundList->value[i];
317
 
 
318
 
        Vect_get_line_areas(Map, line, &left, &right);
319
 
        G_debug(4, "boundary = %d left = %d right = %d", line, left, right);
320
 
 
321
 
        if (left > 0) {
322
 
            dig_list_add(List, left);
323
 
        }
324
 
        else if (left < 0) {    /* island */
325
 
            area = Vect_get_isle_area(Map, abs(left));
326
 
            G_debug(4, "  left island -> area = %d", area);
327
 
            if (area > 0)
328
 
                dig_list_add(List, area);
329
 
        }
330
 
 
331
 
        if (right > 0) {
332
 
            dig_list_add(List, right);
333
 
        }
334
 
        else if (right < 0) {   /* island */
335
 
            area = Vect_get_isle_area(Map, abs(right));
336
 
            G_debug(4, "  right island -> area = %d", area);
337
 
            if (area > 0)
338
 
                dig_list_add(List, area);
339
 
        }
340
 
    }
341
 
 
342
 
    /* But the Polygon may be completely inside the area (only one), in that case 
343
 
     * we find the area by one polygon point and add it to the list */
344
 
    area = Vect_find_area(Map, Polygon->x[0], Polygon->y[0]);
345
 
    if (area > 0)
346
 
        dig_list_add(List, area);
347
 
 
348
 
    G_debug(3, "  %d areas selected by polygon", List->n_values);
349
 
 
350
 
    return List->n_values;
 
17
#include <unistd.h>
 
18
#include <sys/stat.h>
 
19
#include <string.h>
 
20
#include <grass/vector.h>
 
21
#include <grass/glocale.h>
 
22
 
 
23
/*!
 
24
   \brief Initialize spatial index structure
 
25
 
 
26
   \param si pointer to spatial index structure
 
27
 
 
28
   \return void
 
29
 */
 
30
void Vect_spatial_index_init(struct spatial_index * si, int with_z)
 
31
{
 
32
    G_debug(1, "Vect_spatial_index_init()");
 
33
 
 
34
    si->si_tree = RTreeCreateTree(-1, 0, 2 + (with_z != 0));
 
35
}
 
36
 
 
37
/*!
 
38
   \brief Destroy existing spatial index
 
39
 
 
40
   Vect_spatial_index_init() must be call before new use.
 
41
 
 
42
   \param si pointer to spatial index structure
 
43
 
 
44
   \return void
 
45
 */
 
46
void Vect_spatial_index_destroy(struct spatial_index * si)
 
47
{
 
48
    G_debug(1, "Vect_spatial_index_destroy()");
 
49
 
 
50
    RTreeDestroyTree(si->si_tree);
 
51
}
 
52
 
 
53
/*!
 
54
   \brief Add a new item to spatial index structure
 
55
 
 
56
   \param[in,out] si pointer to spatial index structure
 
57
   \param id item identifier
 
58
   \param box pointer to item bounding box
 
59
 
 
60
   \return void
 
61
 */
 
62
void Vect_spatial_index_add_item(struct spatial_index * si, int id,
 
63
                                 const struct bound_box * box)
 
64
{
 
65
    static struct RTree_Rect rect;
 
66
    static int rect_init = 0;
 
67
 
 
68
    if (!rect_init) {
 
69
        rect.boundary = G_malloc(si->si_tree->nsides_alloc * sizeof(RectReal));
 
70
        rect_init = si->si_tree->nsides_alloc;
 
71
    }
 
72
 
 
73
    G_debug(3, "Vect_spatial_index_add_item(): id = %d", id);
 
74
 
 
75
    rect.boundary[0] = box->W;
 
76
    rect.boundary[1] = box->S;
 
77
    rect.boundary[2] = box->B;
 
78
    rect.boundary[3] = box->E;
 
79
    rect.boundary[4] = box->N;
 
80
    rect.boundary[5] = box->T;
 
81
    RTreeInsertRect(&rect, id, si->si_tree);
 
82
}
 
83
 
 
84
/*!
 
85
   \brief Delete item from spatial index structure
 
86
 
 
87
   \param[in,out] si pointer to spatial index structure
 
88
   \param id item identifier
 
89
 
 
90
   \return void
 
91
 */
 
92
void Vect_spatial_index_del_item(struct spatial_index * si, int id,
 
93
                                 const struct bound_box * box)
 
94
{
 
95
    int ret;
 
96
    static struct RTree_Rect rect;
 
97
    static int rect_init = 0;
 
98
 
 
99
    if (!rect_init) {
 
100
        rect.boundary = G_malloc(si->si_tree->nsides_alloc * sizeof(RectReal));
 
101
        rect_init = si->si_tree->nsides_alloc;
 
102
    }
 
103
 
 
104
    G_debug(3, "Vect_spatial_index_del_item(): id = %d", id);
 
105
 
 
106
    rect.boundary[0] = box->W;
 
107
    rect.boundary[1] = box->S;
 
108
    rect.boundary[2] = box->B;
 
109
    rect.boundary[3] = box->E;
 
110
    rect.boundary[4] = box->N;
 
111
    rect.boundary[5] = box->T;
 
112
 
 
113
    ret = RTreeDeleteRect(&rect, id, si->si_tree);
 
114
 
 
115
    if (ret)
 
116
        G_fatal_error(_("Unable to delete item %d from spatial index"), id);
 
117
}
 
118
 
 
119
/************************* SELECT BY BOX *********************************/
 
120
/* This function is called by  RTreeSearch() to add selected item to the list */
 
121
static int _add_item(int id, const struct RTree_Rect *rect, struct ilist *list)
 
122
{
 
123
    G_ilist_add(list, id);
 
124
    return 1;
 
125
}
 
126
 
 
127
/*!
 
128
   \brief Select items by bounding box to list
 
129
 
 
130
   \param si pointer to spatial index structure
 
131
   \param box bounding box
 
132
   \param[out] list pointer to list where selected items are stored
 
133
 
 
134
   \return number of selected items
 
135
 */
 
136
int Vect_spatial_index_select(const struct spatial_index * si, const struct bound_box * box,
 
137
                              struct ilist *list)
 
138
{
 
139
    static struct RTree_Rect rect;
 
140
    static int rect_init = 0;
 
141
 
 
142
    if (!rect_init) {
 
143
        rect.boundary = G_malloc(si->si_tree->nsides_alloc * sizeof(RectReal));
 
144
        rect_init = si->si_tree->nsides_alloc;
 
145
    }
 
146
 
 
147
    Vect_reset_list(list);
 
148
 
 
149
    rect.boundary[0] = box->W;
 
150
    rect.boundary[1] = box->S;
 
151
    rect.boundary[2] = box->B;
 
152
    rect.boundary[3] = box->E;
 
153
    rect.boundary[4] = box->N;
 
154
    rect.boundary[5] = box->T;
 
155
    RTreeSearch(si->si_tree, &rect, (void *)_add_item, list);
 
156
 
 
157
    G_debug(3, "Vect_spatial_index_select(): %d items selected", list->n_values);
 
158
 
 
159
    return list->n_values;
351
160
}