source: grass/trunk/vector/v.distance/distance.c

Last change on this file was 66460, checked in by wenzeslaus, 9 years ago

v.distance: use bbox in 2D and 3D accordingly (see #2734)

Check if we are considering z and use appropriate point in box function. This should
fix #2734 where point inside an area is classified as outsie of that area when point
3D but the area is 2D. The expected behavior is reduction to 2D.

The test works with this commit but its 3D part fails with code prior to this commit.

This commit depens on r66459 which adds Vect_point_in_box_2d() function.
  • Property svn:eol-style set to native
  • Property svn:mime-type set to text/x-csrc
File size: 11.3 KB
Line 
1#include <math.h>
2#include <grass/gis.h>
3#include <grass/glocale.h>
4#include <grass/vector.h>
5#include "local_proto.h"
6
7dist_func *line_distance;
8
9int get_line_box(const struct line_pnts *Points,
10 struct bound_box *box)
11{
12 int i;
13
14 if (Points->n_points == 0) {
15 box->E = box->W = box->N = box->S = box->T = box->B = 0.0 / 0.0;
16 return 0;
17 }
18
19 box->E = box->W = Points->x[0];
20 box->N = box->S = Points->y[0];
21 box->T = box->B = Points->z[0];
22
23 for (i = 1; i < Points->n_points; i++) {
24 if (box->E < Points->x[i])
25 box->E = Points->x[i];
26 if (box->W > Points->x[i])
27 box->W = Points->x[i];
28 if (box->N < Points->y[i])
29 box->N = Points->y[i];
30 if (box->S > Points->y[i])
31 box->S = Points->y[i];
32 if (box->T < Points->z[i])
33 box->T = Points->z[i];
34 if (box->B > Points->z[i])
35 box->B = Points->z[i];
36 }
37
38 return 1;
39}
40
41/* segment angle */
42double sangle(struct line_pnts *Points, int segment)
43{
44 double dx, dy;
45
46 if (Points->n_points < 2 || segment < 1)
47 return -9;
48 if (segment >= Points->n_points)
49 G_fatal_error(_("Invalid segment number %d for %d points"),
50 segment, Points->n_points);
51
52 dx = Points->x[segment] - Points->x[segment - 1];
53 dy = Points->y[segment] - Points->y[segment - 1];
54
55 return atan2(dy, dx);
56}
57
58/* calculate distance parameters between two primitives
59 * return 1 point to point
60 * return 2 point to line
61 * return 1 line to line
62 */
63int line2line(struct line_pnts *FPoints, int ftype,
64 struct line_pnts *TPoints, int ttype,
65 double *fx, double *fy, double *fz,
66 double *falong, double *fangle,
67 double *tx, double *ty, double *tz,
68 double *talong, double *tangle,
69 double *dist,
70 int with_z)
71{
72 int i, fseg, tseg, tmp_seg;
73 double tmp_dist, tmp_x, tmp_y, tmp_z, tmp_along;
74 int ret = 1;
75 static struct line_pnts *iPoints = NULL;
76
77 if (!iPoints) {
78 iPoints = Vect_new_line_struct();
79 }
80
81 *dist = PORT_DOUBLE_MAX;
82
83 /* fangle and tangle are angles in radians, counter clockwise from x axis
84 * initialize to invalid angle */
85 *fangle = *tangle = -9.;
86 *falong = *talong = 0.;
87
88 *fx = FPoints->x[0];
89 *fy = FPoints->y[0];
90 *fz = FPoints->z[0];
91
92 *tx = TPoints->x[0];
93 *ty = TPoints->y[0];
94 *tz = TPoints->z[0];
95
96 tmp_z = 0;
97
98 /* point -> point */
99 if ((ftype & GV_POINTS) && (ttype & GV_POINTS)) {
100 line_distance(TPoints, FPoints->x[0], FPoints->y[0],
101 FPoints->z[0], with_z, tx, ty, tz, dist,
102 NULL, talong);
103 }
104
105 /* point -> line and line -> line */
106 if ((ttype & GV_LINES)) {
107
108 fseg = tseg = 0;
109 /* calculate the min distance between each point in fline with tline */
110 for (i = 0; i < FPoints->n_points; i++) {
111
112 tmp_seg = line_distance(TPoints, FPoints->x[i],
113 FPoints->y[i], FPoints->z[i],
114 with_z, &tmp_x, &tmp_y, &tmp_z,
115 &tmp_dist, NULL, &tmp_along);
116 if (*dist > tmp_dist) {
117 *dist = tmp_dist;
118 *fx = FPoints->x[i];
119 *fy = FPoints->y[i];
120 *fz = FPoints->z[i];
121 *tx = tmp_x;
122 *ty = tmp_y;
123 *tz = tmp_z;
124 *talong = tmp_along;
125 tseg = tmp_seg;
126 fseg = i + 1;
127 }
128 }
129 *tangle = sangle(TPoints, tseg);
130
131 if (FPoints->n_points > 1 && fseg > 0) {
132 int np = FPoints->n_points;
133
134 fseg--;
135
136 if (fseg > 0) {
137 FPoints->n_points = fseg + 1;
138 *falong = Vect_line_length(FPoints);
139 FPoints->n_points = np;
140 }
141 np = fseg;
142 if (fseg == 0)
143 fseg = 1;
144 *fangle = sangle(FPoints, fseg);
145 }
146
147 ret++;
148 }
149
150 /* line -> point and line -> line */
151 if (ftype & GV_LINES) {
152
153 fseg = tseg = 0;
154
155 /* calculate the min distance between each point in tline with fline */
156 for (i = 0; i < TPoints->n_points; i++) {
157
158 tmp_seg = line_distance(FPoints, TPoints->x[i],
159 TPoints->y[i], TPoints->z[i],
160 with_z, &tmp_x, &tmp_y, &tmp_z,
161 &tmp_dist, NULL, &tmp_along);
162 if (*dist > tmp_dist) {
163 *dist = tmp_dist;
164 *fx = tmp_x;
165 *fy = tmp_y;
166 *fz = tmp_z;
167 *falong = tmp_along;
168 *tx = TPoints->x[i];
169 *ty = TPoints->y[i];
170 *tz = TPoints->z[i];
171 fseg = tmp_seg;
172 tseg = i + 1;
173 }
174 }
175 *fangle = sangle(FPoints, fseg);
176
177 if (TPoints->n_points > 1 && tseg > 0) {
178 int np = TPoints->n_points;
179
180 tseg--;
181
182 if (tseg > 0) {
183 TPoints->n_points = tseg + 1;
184 *talong = Vect_line_length(TPoints);
185 TPoints->n_points = np;
186 }
187 np = tseg;
188 if (tseg == 0)
189 tseg = 1;
190 *tangle = sangle(TPoints, tseg);
191 }
192
193 ret++;
194
195 if ((ttype & GV_LINES) && *dist > 0) {
196 /* check for line intersection */
197 struct bound_box fbox, tbox;
198
199 get_line_box(FPoints, &fbox);
200 get_line_box(TPoints, &tbox);
201
202 if (Vect_box_overlap(&fbox, &tbox)) {
203 Vect_reset_line(iPoints);
204 Vect_line_get_intersections(FPoints, TPoints, iPoints, with_z);
205 if (iPoints->n_points) {
206 *dist = 0;
207 *fx = *tx = iPoints->x[0];
208 *fy = *ty = iPoints->y[0];
209 *fz = *tz = iPoints->z[0];
210
211 /* falong, talong */
212 fseg = line_distance(FPoints, iPoints->x[0],
213 iPoints->y[0], iPoints->z[0],
214 with_z, NULL, NULL, NULL,
215 NULL, NULL, falong);
216 tseg = line_distance(TPoints, iPoints->x[0],
217 iPoints->y[0], iPoints->z[0],
218 with_z, NULL, NULL, NULL,
219 NULL, NULL, talong);
220 /* fangle, tangle */
221 *fangle = sangle(FPoints, fseg);
222 *tangle = sangle(TPoints, tseg);
223 }
224 }
225 }
226 }
227
228 return ret;
229}
230
231/* shortest distance between line and area
232 * return 1 inside area
233 * return 2 inside isle of area
234 * return 3 outside area */
235int line2area(const struct Map_info *To,
236 struct line_pnts *Points, int type,
237 int area, const struct bound_box *abox,
238 double *fx, double *fy, double *fz,
239 double *falong, double *fangle,
240 double *tx, double *ty, double *tz,
241 double *talong, double *tangle,
242 double *dist,
243 int with_z)
244{
245 int i, j;
246 double tmp_dist;
247 int isle, nisles;
248 int all_inside_outer, all_outside_outer, all_outside_inner;
249 static struct line_pnts *aPoints = NULL;
250 static struct line_pnts **iPoints = NULL;
251 static struct bound_box *ibox = NULL;
252 static int isle_alloc = 0;
253
254 if (!aPoints)
255 aPoints = Vect_new_line_struct();
256
257 *dist = PORT_DOUBLE_MAX;
258
259 /* fangle and tangle are angles in radians, counter clockwise from x axis
260 * initialize to invalid angle */
261 *fangle = *tangle = -9.;
262 *falong = *talong = 0.;
263
264 *fx = Points->x[0];
265 *fy = Points->y[0];
266 *fz = Points->z[0];
267
268 *tx = Points->x[0];
269 *ty = Points->y[0];
270 *tz = Points->z[0];
271
272 Vect_get_area_points(To, area, aPoints);
273 nisles = Vect_get_area_num_isles(To, area);
274
275 if (nisles > isle_alloc) {
276 iPoints = G_realloc(iPoints, nisles * sizeof(struct line_pnts *));
277 ibox = G_realloc(ibox, nisles * sizeof(struct bound_box));
278 for (i = isle_alloc; i < nisles; i++)
279 iPoints[i] = Vect_new_line_struct();
280 isle_alloc = nisles;
281 }
282 for (i = 0; i < nisles; i++) {
283 isle = Vect_get_area_isle(To, area, i);
284 Vect_get_isle_points(To, isle, iPoints[i]);
285 Vect_get_isle_box(To, isle, &ibox[i]);
286 }
287
288 /* inside area ? */
289 all_inside_outer = all_outside_outer = 1;
290 all_outside_inner = 1;
291
292 int in_box;
293 for (i = 0; i < Points->n_points; i++) {
294 if (with_z)
295 in_box = Vect_point_in_box(Points->x[i], Points->y[i],
296 Points->z[i], abox);
297 else
298 in_box = Vect_point_in_box_2d(Points->x[i], Points->y[i], abox);
299 if (in_box) {
300
301 int poly;
302
303 poly = Vect_point_in_poly(Points->x[i], Points->y[i], aPoints);
304
305 if (poly > 0) {
306 /* inside outer ring */
307 all_outside_outer = 0;
308 }
309 else {
310 /* outside outer ring */
311 all_inside_outer = 0;
312 }
313
314 /* exactly on boundary */
315 if (poly == 2) {
316 line2line(Points, type, aPoints, GV_BOUNDARY,
317 fx, fy, fz, falong, fangle,
318 tx, ty, tz, talong, tangle,
319 dist, with_z);
320
321 *talong = 0;
322 *tangle = -9;
323
324 return 1;
325 }
326 /* inside outer ring */
327 else if (poly == 1) {
328 int inside_isle = 0;
329
330 for (j = 0; j < nisles; j++) {
331 if (with_z)
332 in_box = Vect_point_in_box(Points->x[i], Points->y[i],
333 Points->z[i], &ibox[j]);
334 else
335 in_box = Vect_point_in_box_2d(Points->x[i],
336 Points->y[i], &ibox[j]);
337 if (in_box) {
338
339 poly = Vect_point_in_poly(Points->x[i], Points->y[i], iPoints[j]);
340
341 /* inside or exactly on boundary */
342 if (poly > 0) {
343 double tmp_fx, tmp_fy, tmp_fz, tmp_fangle, tmp_falong;
344 double tmp_tx, tmp_ty, tmp_tz, tmp_tangle, tmp_talong;
345
346 /* pass all points of the line,
347 * this will catch an intersection */
348 line2line(Points, type, iPoints[j], GV_BOUNDARY,
349 &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
350 &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
351 &tmp_dist, with_z);
352
353 if (*dist > tmp_dist) {
354 *dist = tmp_dist;
355
356 *fx = tmp_fx;
357 *fy = tmp_fy;
358 *fz = tmp_fz;
359 *falong = tmp_falong;
360 *fangle = tmp_fangle;
361
362 *tx = tmp_tx;
363 *ty = tmp_ty;
364 *tz = tmp_tz;
365 *talong = 0;
366 *tangle = tmp_tangle;
367 }
368
369 if (poly == 1) /* excludes isle boundary */
370 inside_isle = 1;
371
372 }
373 }
374 if (*dist == 0)
375 break;
376 }
377 /* inside area (inside outer ring, outside inner rings
378 * or exactly on one of the inner rings) */
379 if (!inside_isle) {
380 *fx = Points->x[i];
381 *fy = Points->y[i];
382 *fz = Points->z[i];
383
384 *tx = Points->x[i];
385 *ty = Points->y[i];
386 *tz = Points->z[i];
387
388 *fangle = *tangle = -9.;
389 *falong = *talong = 0.;
390
391 *dist = 0;
392
393 return 1;
394 }
395 else {
396 /* inside one of the islands */
397 all_outside_inner = 0;
398 if (*dist == 0) {
399 /* the line intersected with the isle boundary
400 * -> line is partially inside the area */
401
402 *fangle = *tangle = -9.;
403 *falong = *talong = 0.;
404
405 return 1;
406 }
407 /* else continue with next point */
408 }
409 } /* end inside outer ring */
410 }
411 else {
412 /* point not in box of outer ring */
413 all_inside_outer = 0;
414 }
415 /* exactly on boundary */
416 if (*dist == 0)
417 return 1;
418 }
419
420 /* if all points are inside the outer ring and inside inner rings,
421 * there could still be an intersection with one of the inner rings */
422 if (all_inside_outer) {
423 if (all_outside_inner) {
424 /* at least one point is really inside the area!
425 * that should have been detected above */
426 G_fatal_error(_("At least one point is really inside the area!"));
427 }
428 /* else all points are inside one of the area isles
429 * and we already have the minimum distance */
430 return 2;
431 }
432
433 /* if at least one point was found to be inside the outer ring,
434 * but no point really inside the area,
435 * and at least one point outside,
436 * then there must be an intersection of the line with both
437 * the outer ring and one of the isle boundaries */
438
439 /* if all line points are outside of the area,
440 * intersection is still possible */
441
442 line2line(Points, type, aPoints, GV_BOUNDARY,
443 fx, fy, fz, falong, fangle,
444 tx, ty, tz, talong, tangle,
445 dist, with_z);
446
447 *talong = 0;
448
449 if (*dist == 0)
450 return 1;
451
452 return 3;
453}
Note: See TracBrowser for help on using the repository browser.