GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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 #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 
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
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
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
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  */
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  */
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 */
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 }
void Vect_destroy_boxlist(struct boxlist *)
Frees all memory associated with a struct boxlist, including the struct itself.
#define G_malloc(n)
Definition: defs/gis.h:112
Bounding box.
Definition: dig_structs.h:65
int * id
Array of ids.
Definition: dig_structs.h:1755
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:125
int Vect_point_in_island(double, double, const struct Map_info *, int, struct bound_box *)
Determines if a point (X,Y) is inside an island.
Definition: Vlib/poly.c:925
double W
West.
Definition: dig_structs.h:82
int Vect_reset_list(struct ilist *)
Reset ilist structure.
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:171
int Vect_select_nodes_by_box(struct Map_info *, const struct bound_box *, struct ilist *)
Select nodes by box.
Definition: sindex.c:194
int Vect_point_in_area_outer_ring(double, double, const 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:856
struct P_area ** Area
Array of areas.
Definition: dig_structs.h:891
int n_points
Number of points.
Definition: dig_structs.h:1692
int Vect_val_in_list(const struct ilist *, int)
Find a given item in the list.
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:651
double E
East.
Definition: dig_structs.h:78
plus_t * isles
1st generation interior islands
Definition: dig_structs.h:1640
double Vect_points_distance(double, double, double, double, double, double, int)
Calculate distance of 2 points.
Definition: line.c:898
#define NULL
Definition: ccmath.h:32
#define x
int n_values
Number of items in the list.
Definition: dig_structs.h:1767
struct bound_box * box
Array of bounding boxes.
Definition: dig_structs.h:1759
int G_begin_polygon_area_calculations(void)
Begin polygon area calculations.
Definition: gis/area.c:120
int Vect_get_node_coor(const struct Map_info *, int, double *, double *, double *)
Get node coordinates.
Definition: level_two.c:278
double N
North.
Definition: dig_structs.h:70
int Vect_find_area(struct Map_info *Map, double x, double y)
Find the nearest area.
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
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:34
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
Basic topology-related info.
Definition: dig_structs.h:784
double b
Definition: r_raster.c:39
void Vect_destroy_list(struct ilist *)
Frees all memory associated with a struct ilist, including the struct itself.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
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_get_isle_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points for given isle.
double B
Bottom.
Definition: dig_structs.h:90
struct Plus_head plus
Plus info (topology, version, ...)
Definition: dig_structs.h:1286
#define HUGE_VAL
double T
Top.
Definition: dig_structs.h:86
int Vect_get_isle_box(const struct Map_info *, int, struct bound_box *)
Get bounding box of isle.
plus_t n_isles
Number of islands inside.
Definition: dig_structs.h:1632
Vector map info.
Definition: dig_structs.h:1259
#define PORT_DOUBLE_MAX
Limits for portable types.
Definition: dig_defines.h:66
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
Area (topology) info.
Definition: dig_structs.h:1605
int Vect_find_node(struct Map_info *Map, double ux, double uy, double uz, double maxdist, int with_z)
Find the nearest node.
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
int Vect_find_island(struct Map_info *Map, double x, double y)
Find the nearest island.
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
List of bounding boxes with id.
Definition: dig_structs.h:1750
double S
South.
Definition: dig_structs.h:74
#define G_realloc(p, n)
Definition: defs/gis.h:114
List of integers.
Definition: gis.h:700
int * value
Array of values.
Definition: gis.h:705
struct boxlist * Vect_new_boxlist(int)
Creates and initializes a struct boxlist.
int Vect_read_line(const struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int G_debug(int, const char *,...) __attribute__((format(printf
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.
double G_area_of_polygon(const double *, const double *, int)
Area in square meters of polygon.
Definition: gis/area.c:159