GRASS Programmer's Manual
6.4.1(2011)
|
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 /* Select all nodes in box */ 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 /* find nearest */ 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 /* Check if in max distance */ 00086 if (cur_dist <= maxdist) 00087 return (NList->value[node]); 00088 else 00089 return 0; 00090 } 00091 00108 /* original dig_point_to_line() in grass50 */ 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 /* No more needed */ 00201 /* 00202 Line = Plus->Line[line]; 00203 if ( Line == NULL ) continue; 00204 if ( !(type & Line->type) ) continue; 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 /* TODO */ 00220 /* choice = dig_center_check (map->Line, choice, a, ux, uy); */ 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 /* original dig_point_to_area() in grass50 */ 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 /* select areas by box */ 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 /* original dig_point_to_area() in grass50 */ 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 /* select islands by box */ 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) { /* inside */ 00328 if (current > 0) { /* not first */ 00329 if (current_size == -1) { /* second */ 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 { /* first */ 00347 current = island; 00348 } 00349 } 00350 } 00351 00352 return current; 00353 }