00001
00018 #include <math.h>
00019 #include <grass/gis.h>
00020 #include <grass/Vect.h>
00021
00022 #ifndef HUGE_VAL
00023 #define HUGE_VAL 9999999999999.0
00024 #endif
00025
00037 int
00038 Vect_find_node(struct Map_info *Map,
00039 double ux, double uy, double uz, double maxdist, int with_z)
00040 {
00041 int i, nnodes, node;
00042 BOUND_BOX box;
00043 struct ilist *NList;
00044 double x, y, z;
00045 double cur_dist, dist;
00046
00047 G_debug(3, "Vect_find_node() for %f %f %f maxdist = %f", ux, uy, uz,
00048 maxdist);
00049 NList = Vect_new_list();
00050
00051
00052 box.N = uy + maxdist;
00053 box.S = uy - maxdist;
00054 box.E = ux + maxdist;
00055 box.W = ux - maxdist;
00056 if (with_z) {
00057 box.T = uz + maxdist;
00058 box.B = uz - maxdist;
00059 }
00060 else {
00061 box.T = HUGE_VAL;
00062 box.B = -HUGE_VAL;
00063 }
00064
00065 nnodes = Vect_select_nodes_by_box(Map, &box, NList);
00066 G_debug(3, " %d nodes in box", nnodes);
00067
00068 if (nnodes == 0)
00069 return 0;
00070
00071
00072 cur_dist = PORT_DOUBLE_MAX;
00073 node = 0;
00074 for (i = 0; i < nnodes; i++) {
00075 Vect_get_node_coor(Map, NList->value[i], &x, &y, &z);
00076 dist = Vect_points_distance(ux, uy, uz, x, y, z, with_z);
00077 if (dist < cur_dist) {
00078 cur_dist = dist;
00079 node = i;
00080 }
00081 }
00082 G_debug(3, " nearest node %d in distance %f", NList->value[node],
00083 cur_dist);
00084
00085
00086 if (cur_dist <= maxdist)
00087 return (NList->value[node]);
00088 else
00089 return 0;
00090 }
00091
00108
00109 int
00110 Vect_find_line(struct Map_info *map,
00111 double ux, double uy, double uz,
00112 int type, double maxdist, int with_z, int exclude)
00113 {
00114 int line;
00115 struct ilist *exclude_list;
00116
00117 exclude_list = Vect_new_list();
00118
00119 Vect_list_append(exclude_list, exclude);
00120
00121 line = Vect_find_line_list(map, ux, uy, uz,
00122 type, maxdist, with_z, exclude_list, NULL);
00123
00124 Vect_destroy_list(exclude_list);
00125
00126 return line;
00127 }
00128
00144 int
00145 Vect_find_line_list(struct Map_info *map,
00146 double ux, double uy, double uz,
00147 int type, double maxdist, int with_z,
00148 struct ilist *exclude, struct ilist *found)
00149 {
00150 int choice;
00151 double new_dist;
00152 double cur_dist;
00153 int gotone;
00154 int i, line;
00155 static struct line_pnts *Points;
00156 static int first_time = 1;
00157 struct Plus_head *Plus;
00158 BOUND_BOX box;
00159 struct ilist *List;
00160
00161 G_debug(3, "Vect_find_line_list() for %f %f %f type = %d maxdist = %f",
00162 ux, uy, uz, type, maxdist);
00163
00164 if (first_time) {
00165 Points = Vect_new_line_struct();
00166 first_time = 0;
00167 }
00168
00169 Plus = &(map->plus);
00170 gotone = 0;
00171 choice = 0;
00172 cur_dist = HUGE_VAL;
00173
00174 box.N = uy + maxdist;
00175 box.S = uy - maxdist;
00176 box.E = ux + maxdist;
00177 box.W = ux - maxdist;
00178 if (with_z) {
00179 box.T = uz + maxdist;
00180 box.B = uz - maxdist;
00181 }
00182 else {
00183 box.T = PORT_DOUBLE_MAX;
00184 box.B = -PORT_DOUBLE_MAX;
00185 }
00186
00187 List = Vect_new_list();
00188
00189 if (found)
00190 Vect_reset_list(found);
00191
00192 Vect_select_lines_by_box(map, &box, type, List);
00193 for (i = 0; i < List->n_values; i++) {
00194 line = List->value[i];
00195 if (Vect_val_in_list(exclude, line)) {
00196 G_debug(3, " line = %d exclude", line);
00197 continue;
00198 }
00199
00200
00201
00202
00203
00204
00205
00206
00207 Vect_read_line(map, Points, NULL, line);
00208
00209 Vect_line_distance(Points, ux, uy, uz, with_z, NULL, NULL, NULL,
00210 &new_dist, NULL, NULL);
00211 G_debug(3, " line = %d distance = %f", line, new_dist);
00212
00213 if (found && new_dist <= maxdist) {
00214 Vect_list_append(found, line);
00215 }
00216
00217 if ((++gotone == 1) || (new_dist <= cur_dist)) {
00218 if (new_dist == cur_dist) {
00219
00220
00221 continue;
00222 }
00223
00224 choice = line;
00225 cur_dist = new_dist;
00226 }
00227 }
00228
00229 G_debug(3, "min distance found = %f", cur_dist);
00230 if (cur_dist > maxdist)
00231 choice = 0;
00232
00233 Vect_destroy_list(List);
00234
00235 return (choice);
00236 }
00237
00247
00248 int Vect_find_area(struct Map_info *Map, double x, double y)
00249 {
00250 int i, ret, area;
00251 static int first = 1;
00252 BOUND_BOX box;
00253 static struct ilist *List;
00254
00255 G_debug(3, "Vect_find_area() x = %f y = %f", x, y);
00256
00257 if (first) {
00258 List = Vect_new_list();
00259 first = 0;
00260 }
00261
00262
00263 box.E = x;
00264 box.W = x;
00265 box.N = y;
00266 box.S = y;
00267 box.T = PORT_DOUBLE_MAX;
00268 box.B = -PORT_DOUBLE_MAX;
00269 Vect_select_areas_by_box(Map, &box, List);
00270 G_debug(3, " %d areas selected by box", List->n_values);
00271
00272 for (i = 0; i < List->n_values; i++) {
00273 area = List->value[i];
00274 ret = Vect_point_in_area(Map, area, x, y);
00275
00276 G_debug(3, " area = %d Vect_point_in_area() = %d", area, ret);
00277
00278 if (ret >= 1)
00279 return (area);
00280 }
00281
00282 return 0;
00283 }
00284
00294
00295 int Vect_find_island(struct Map_info *Map, double x, double y)
00296 {
00297 int i, ret, island, current, current_size, size;
00298 static int first = 1;
00299 BOUND_BOX box;
00300 static struct ilist *List;
00301 static struct line_pnts *Points;
00302
00303 G_debug(3, "Vect_find_island() x = %f y = %f", x, y);
00304
00305 if (first) {
00306 List = Vect_new_list();
00307 Points = Vect_new_line_struct();
00308 first = 0;
00309 }
00310
00311
00312 box.E = x;
00313 box.W = x;
00314 box.N = y;
00315 box.S = y;
00316 box.T = PORT_DOUBLE_MAX;
00317 box.B = -PORT_DOUBLE_MAX;
00318 Vect_select_isles_by_box(Map, &box, List);
00319 G_debug(3, " %d islands selected by box", List->n_values);
00320
00321 current_size = -1;
00322 current = 0;
00323 for (i = 0; i < List->n_values; i++) {
00324 island = List->value[i];
00325 ret = Vect_point_in_island(x, y, Map, island);
00326
00327 if (ret >= 1) {
00328 if (current > 0) {
00329 if (current_size == -1) {
00330 G_begin_polygon_area_calculations();
00331 Vect_get_isle_points(Map, current, Points);
00332 current_size =
00333 G_area_of_polygon(Points->x, Points->y,
00334 Points->n_points);
00335 }
00336
00337 Vect_get_isle_points(Map, island, Points);
00338 size =
00339 G_area_of_polygon(Points->x, Points->y, Points->n_points);
00340
00341 if (size < current_size) {
00342 current = island;
00343 current_size = size;
00344 }
00345 }
00346 else {
00347 current = island;
00348 }
00349 }
00350 }
00351
00352 return current;
00353 }