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