1 /*!
2 * \file lib/vector/Vlib/find.c
3 *
4 * \brief Vector library - Find nearest vector feature.
5 *
6 * Higher level functions for reading/writing/manipulating vectors.
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.
12
13 * \author Original author CERL, probably Dave Gerdes or Mike
14 * Higgins.
15 * \author Update to GRASS 5.7 Radim Blazek and David D. Gray.
16 */
17 #include <stdlib.h>
18 #include <math.h>
19 #include <grass/vector.h>
20
21 #ifndef HUGE_VAL
22 #define HUGE_VAL 9999999999999.0
23 #endif
24
25
26 /* for qsort */
27
28 typedef struct {
29 int i;
30 double size;
31 struct bound_box box;
32 } BOX_SIZE;
33
sort_by_size(const void * a,const void * b)34 static int sort_by_size(const void *a, const void *b)
35 {
36 BOX_SIZE *as = (BOX_SIZE *)a;
37 BOX_SIZE *bs = (BOX_SIZE *)b;
38
39 if (as->size < bs->size)
40 return -1;
41
42 return (as->size > bs->size);
43 }
44
45
46 /*!
47 * \brief Find the nearest node.
48 *
49 * \param Map vector map
50 * \param ux,uy,uz point coordinates
51 * \param maxdist max distance from the line
52 * \param with_z 3D (WITH_Z, WITHOUT_Z)
53 *
54 * \return number of nearest node
55 * \return 0 if not found
56 */
57 int
Vect_find_node(struct Map_info * Map,double ux,double uy,double uz,double maxdist,int with_z)58 Vect_find_node(struct Map_info *Map,
59 double ux, double uy, double uz, double maxdist, int with_z)
60 {
61 int i, nnodes, node;
62 struct bound_box box;
63 struct ilist *NList;
64 double x, y, z;
65 double cur_dist, dist;
66
67 G_debug(3, "Vect_find_node() for %f %f %f maxdist = %f", ux, uy, uz,
68 maxdist);
69 NList = Vect_new_list();
70
71 /* Select all nodes in box */
72 box.N = uy + maxdist;
73 box.S = uy - maxdist;
74 box.E = ux + maxdist;
75 box.W = ux - maxdist;
76 if (with_z) {
77 box.T = uz + maxdist;
78 box.B = uz - maxdist;
79 }
80 else {
81 box.T = HUGE_VAL;
82 box.B = -HUGE_VAL;
83 }
84
85 nnodes = Vect_select_nodes_by_box(Map, &box, NList);
86 G_debug(3, " %d nodes in box", nnodes);
87
88 if (nnodes == 0)
89 return 0;
90
91 /* find nearest */
92 cur_dist = PORT_DOUBLE_MAX;
93 node = 0;
94 for (i = 0; i < nnodes; i++) {
95 Vect_get_node_coor(Map, NList->value[i], &x, &y, &z);
96 dist = Vect_points_distance(ux, uy, uz, x, y, z, with_z);
97 if (dist < cur_dist) {
98 cur_dist = dist;
99 node = i;
100 }
101 }
102 G_debug(3, " nearest node %d in distance %f", NList->value[node],
103 cur_dist);
104
105 /* Check if in max distance */
106 if (cur_dist <= maxdist)
107 return (NList->value[node]);
108 else
109 return 0;
110 }
111
112 /*!
113 * \brief Find the nearest line.
114 *
115 * \param map vector map
116 * \param ux,uy,uz points coordinates
117 * \param type feature type (GV_LINE, GV_POINT, GV_BOUNDARY or GV_CENTROID)
118 * if only want to search certain types of lines or -1 if search all lines
119 * \param maxdist max distance from the line
120 * \param with_z 3D (WITH_Z, WITHOUT_Z)
121 * \param exclude if > 0 number of line which should be excluded from selection.
122 * May be useful if we need line nearest to other one.
123 *
124 * \return number of nearest line
125 * \return 0 if not found
126 *
127 */
128 int
Vect_find_line(struct Map_info * map,double ux,double uy,double uz,int type,double maxdist,int with_z,int exclude)129 Vect_find_line(struct Map_info *map,
130 double ux, double uy, double uz,
131 int type, double maxdist, int with_z, int exclude)
132 {
133 int line;
134 struct ilist *exclude_list;
135
136 exclude_list = Vect_new_list();
137
138 Vect_list_append(exclude_list, exclude);
139
140 line = Vect_find_line_list(map, ux, uy, uz,
141 type, maxdist, with_z, exclude_list, NULL);
142
143 Vect_destroy_list(exclude_list);
144
145 return line;
146 }
147
148 /*!
149 * \brief Find the nearest line(s).
150 *
151 * \param map vector map
152 * \param ux,uy,uz points coordinates
153 * \param type feature type (GV_LINE, GV_POINT, GV_BOUNDARY or GV_CENTROID)
154 * if only want to search certain types of lines or -1 if search all lines
155 * \param maxdist max distance from the line
156 * \param with_z 3D (WITH_Z, WITHOUT_Z)
157 * \param exclude list of lines which should be excluded from selection
158 * \param found list of found lines (or NULL)
159 *
160 * \return number of nearest line
161 * \return 0 if not found
162 */
163 int
Vect_find_line_list(struct Map_info * map,double ux,double uy,double uz,int type,double maxdist,int with_z,const struct ilist * exclude,struct ilist * found)164 Vect_find_line_list(struct Map_info *map,
165 double ux, double uy, double uz,
166 int type, double maxdist, int with_z,
167 const struct ilist *exclude, struct ilist *found)
168 {
169 int choice;
170 double new_dist;
171 double cur_dist;
172 int gotone;
173 int i, line;
174 static struct line_pnts *Points;
175 static int first_time = 1;
176 struct bound_box box;
177 struct boxlist *List;
178
179 G_debug(3, "Vect_find_line_list() for %f %f %f type = %d maxdist = %f",
180 ux, uy, uz, type, maxdist);
181
182 if (first_time) {
183 Points = Vect_new_line_struct();
184 first_time = 0;
185 }
186
187 gotone = 0;
188 choice = 0;
189 cur_dist = HUGE_VAL;
190
191 box.N = uy + maxdist;
192 box.S = uy - maxdist;
193 box.E = ux + maxdist;
194 box.W = ux - maxdist;
195 if (with_z) {
196 box.T = uz + maxdist;
197 box.B = uz - maxdist;
198 }
199 else {
200 box.T = PORT_DOUBLE_MAX;
201 box.B = -PORT_DOUBLE_MAX;
202 }
203
204 List = Vect_new_boxlist(0);
205
206 if (found)
207 Vect_reset_list(found);
208
209 Vect_select_lines_by_box(map, &box, type, List);
210 for (i = 0; i < List->n_values; i++) {
211 line = List->id[i];
212 if (Vect_val_in_list(exclude, line)) {
213 G_debug(3, " line = %d exclude", line);
214 continue;
215 }
216
217 /* No more needed */
218 /*
219 Line = Plus->Line[line];
220 if ( Line == NULL ) continue;
221 if ( !(type & Line->type) ) continue;
222 */
223
224 Vect_read_line(map, Points, NULL, line);
225
226 Vect_line_distance(Points, ux, uy, uz, with_z, NULL, NULL, NULL,
227 &new_dist, NULL, NULL);
228 G_debug(3, " line = %d distance = %f", line, new_dist);
229
230 if (found && new_dist <= maxdist) {
231 Vect_list_append(found, line);
232 }
233
234 if ((++gotone == 1) || (new_dist <= cur_dist)) {
235 if (new_dist == cur_dist) {
236 /* TODO */
237 /* choice = dig_center_check (map->Line, choice, a, ux, uy); */
238 continue;
239 }
240
241 choice = line;
242 cur_dist = new_dist;
243 }
244 }
245
246 G_debug(3, "min distance found = %f", cur_dist);
247 if (cur_dist > maxdist)
248 choice = 0;
249
250 Vect_destroy_boxlist(List);
251
252 return (choice);
253 }
254
255 /*!
256 * \brief Find the nearest area
257 *
258 * \param Map vector map
259 * \param x,y point coordinates
260 *
261 * \return area number
262 * \return 0 if not found
263 */
Vect_find_area(struct Map_info * Map,double x,double y)264 int Vect_find_area(struct Map_info *Map, double x, double y)
265 {
266 int i, j, ret, area, isle;
267 struct bound_box box;
268 static struct boxlist *List = NULL;
269 static BOX_SIZE *size_list;
270 static int alloc_size_list = 0;
271 const struct Plus_head *Plus;
272 struct P_area *Area;
273
274 G_debug(3, "Vect_find_area() x = %f y = %f", x, y);
275
276 if (!List) {
277 List = Vect_new_boxlist(1);
278 alloc_size_list = 10;
279 size_list = G_malloc(alloc_size_list * sizeof(BOX_SIZE));
280 }
281
282 Plus = &(Map->plus);
283
284 /* select areas by box */
285 box.E = x;
286 box.W = x;
287 box.N = y;
288 box.S = y;
289 box.T = PORT_DOUBLE_MAX;
290 box.B = -PORT_DOUBLE_MAX;
291 Vect_select_areas_by_box(Map, &box, List);
292 G_debug(3, " %d areas selected by box", List->n_values);
293
294 /* sort areas by bbox size
295 * get the smallest area that contains the point
296 * using the bbox size is working because if 2 areas both contain
297 * the point, one of these areas must be inside the other area
298 * which means that the bbox of the outer area must be larger than
299 * the bbox of the inner area, and equal bbox sizes are not possible */
300
301 if (alloc_size_list < List->n_values) {
302 alloc_size_list = List->n_values;
303 size_list = G_realloc(size_list, alloc_size_list * sizeof(BOX_SIZE));
304 }
305
306 for (i = 0; i < List->n_values; i++) {
307 size_list[i].i = List->id[i];
308 box = List->box[i];
309 size_list[i].box = List->box[i];
310 size_list[i].size = (box.N - box.S) * (box.E - box.W);
311 }
312
313 if (List->n_values == 2) {
314 /* simple swap */
315 if (size_list[1].size < size_list[0].size) {
316 size_list[0].i = List->id[1];
317 size_list[1].i = List->id[0];
318 size_list[0].box = List->box[1];
319 size_list[1].box = List->box[0];
320 }
321 }
322 else if (List->n_values > 2)
323 qsort(size_list, List->n_values, sizeof(BOX_SIZE), sort_by_size);
324
325 for (i = 0; i < List->n_values; i++) {
326 area = size_list[i].i;
327 /* outer ring */
328 ret = Vect_point_in_area_outer_ring(x, y, Map, area, &size_list[i].box);
329
330 G_debug(3, " area = %d Vect_point_in_area_outer_ring() = %d", area, ret);
331
332 if (ret >= 1) {
333 /* check if in islands */
334 Area = Plus->Area[area];
335 for (j = 0; j < Area->n_isles; j++) {
336 isle = Area->isles[j];
337 Vect_get_isle_box(Map, isle, &box);
338 ret = Vect_point_in_island(x, y, Map, isle, &box);
339
340 G_debug(3, " area = %d Vect_point_in_island() = %d", area, ret);
341
342 if (ret >= 1) {
343 /* point is not in area
344 * point is also not in any inner area, those have
345 * been tested before (sorted list)
346 * -> area inside island could not be built */
347 return 0;
348 }
349 }
350 return (area);
351 }
352 }
353
354 return 0;
355 }
356
357 /*!
358 * \brief Find the nearest island
359 *
360 * \param Map vector map
361 * \param x,y points coordinates
362 *
363 * \return island number,
364 * \return 0 if not found
365 */
Vect_find_island(struct Map_info * Map,double x,double y)366 int Vect_find_island(struct Map_info *Map, double x, double y)
367 {
368 int i, ret, island, current, current_size, size;
369 static int first = 1;
370 struct bound_box box;
371 static struct boxlist *List;
372 static struct line_pnts *Points;
373
374 /* TODO: sync to Vect_find_area() */
375
376 G_debug(3, "Vect_find_island() x = %f y = %f", x, y);
377
378 if (first) {
379 List = Vect_new_boxlist(1);
380 Points = Vect_new_line_struct();
381 first = 0;
382 }
383
384 /* select islands by box */
385 box.E = x;
386 box.W = x;
387 box.N = y;
388 box.S = y;
389 box.T = PORT_DOUBLE_MAX;
390 box.B = -PORT_DOUBLE_MAX;
391 Vect_select_isles_by_box(Map, &box, List);
392 G_debug(3, " %d islands selected by box", List->n_values);
393
394 current_size = -1;
395 current = 0;
396 for (i = 0; i < List->n_values; i++) {
397 island = List->id[i];
398 ret = Vect_point_in_island(x, y, Map, island, &List->box[i]);
399
400 if (ret >= 1) { /* inside */
401 if (current > 0) { /* not first */
402 if (current_size == -1) { /* second */
403 G_begin_polygon_area_calculations();
404 Vect_get_isle_points(Map, current, Points);
405 current_size =
406 G_area_of_polygon(Points->x, Points->y,
407 Points->n_points);
408 }
409
410 Vect_get_isle_points(Map, island, Points);
411 size =
412 G_area_of_polygon(Points->x, Points->y, Points->n_points);
413
414 if (size < current_size) {
415 current = island;
416 current_size = size;
417 }
418 }
419 else { /* first */
420 current = island;
421 }
422 }
423 }
424
425 return current;
426 }
427