1 /******************************************************************************
2 * $Id$
3 *
4 * Project: MapServer
5 * Purpose: Various geospatial search operations.
6 * Author: Steve Lime and the MapServer team.
7 *
8 * Notes: For information on point in polygon function please see:
9 *
10 * http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
11 *
12 * The appropriate copyright notice accompanies the funtion definition.
13 *
14 ******************************************************************************
15 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
16 *
17 * Permission is hereby granted, free of charge, to any person obtaining a
18 * copy of this software and associated documentation files (the "Software"),
19 * to deal in the Software without restriction, including without limitation
20 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
21 * and/or sell copies of the Software, and to permit persons to whom the
22 * Software is furnished to do so, subject to the following conditions:
23 *
24 * The above copyright notice and this permission notice shall be included in
25 * all copies of this Software or works derived from this Software.
26 *
27 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
28 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
30 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
31 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
32 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
33 * DEALINGS IN THE SOFTWARE.
34 ****************************************************************************/
35
36 #include "mapserver.h"
37
38
39
40 #define LASTVERT(v,n) ((v) == 0 ? n-2 : v-1)
41 #define NEXTVERT(v,n) ((v) == n-2 ? 0 : v+1)
42
43 /*
44 ** Returns MS_TRUE if rectangles a and b overlap
45 */
msRectOverlap(const rectObj * a,const rectObj * b)46 int msRectOverlap(const rectObj *a, const rectObj *b)
47 {
48 if(a->minx > b->maxx) return(MS_FALSE);
49 if(a->maxx < b->minx) return(MS_FALSE);
50 if(a->miny > b->maxy) return(MS_FALSE);
51 if(a->maxy < b->miny) return(MS_FALSE);
52 return(MS_TRUE);
53 }
54
55 /*
56 ** Computes the intersection of two rectangles, updating the first
57 ** to be only the intersection of the two. Returns MS_FALSE if
58 ** the intersection is empty.
59 */
msRectIntersect(rectObj * a,const rectObj * b)60 int msRectIntersect( rectObj *a, const rectObj *b )
61 {
62 if( a->maxx > b->maxx )
63 a->maxx = b->maxx;
64 if( a->minx < b->minx )
65 a->minx = b->minx;
66 if( a->maxy > b->maxy )
67 a->maxy = b->maxy;
68 if( a->miny < b->miny )
69 a->miny = b->miny;
70
71 if( a->maxx < a->minx || b->maxx < b->minx )
72 return MS_FALSE;
73 else
74 return MS_TRUE;
75 }
76
77 /*
78 ** Returns MS_TRUE if rectangle a is contained in rectangle b
79 */
msRectContained(const rectObj * a,const rectObj * b)80 int msRectContained(const rectObj *a, const rectObj *b)
81 {
82 if(a->minx >= b->minx && a->maxx <= b->maxx)
83 if(a->miny >= b->miny && a->maxy <= b->maxy)
84 return(MS_TRUE);
85 return(MS_FALSE);
86 }
87
88 /*
89 ** Merges rect b into rect a. Rect a changes, b does not.
90 */
msMergeRect(rectObj * a,rectObj * b)91 void msMergeRect(rectObj *a, rectObj *b)
92 {
93 a->minx = MS_MIN(a->minx, b->minx);
94 a->maxx = MS_MAX(a->maxx, b->maxx);
95 a->miny = MS_MIN(a->miny, b->miny);
96 a->maxy = MS_MAX(a->maxy, b->maxy);
97 }
98
msPointInRect(const pointObj * p,const rectObj * rect)99 int msPointInRect(const pointObj *p, const rectObj *rect)
100 {
101 if(p->x < rect->minx) return(MS_FALSE);
102 if(p->x > rect->maxx) return(MS_FALSE);
103 if(p->y < rect->miny) return(MS_FALSE);
104 if(p->y > rect->maxy) return(MS_FALSE);
105 return(MS_TRUE);
106 }
107
msPolygonDirection(lineObj * c)108 int msPolygonDirection(lineObj *c)
109 {
110 double mx, my, area;
111 int i, v=0, lv, nv;
112
113 /* first find lowest, rightmost point of polygon */
114 mx = c->point[0].x;
115 my = c->point[0].y;
116
117 for(i=0; i<c->numpoints-1; i++) {
118 if((c->point[i].y < my) || ((c->point[i].y == my) && (c->point[i].x > mx))) {
119 v = i;
120 mx = c->point[i].x;
121 my = c->point[i].y;
122 }
123 }
124
125 lv = LASTVERT(v,c->numpoints);
126 nv = NEXTVERT(v,c->numpoints);
127
128 area = c->point[lv].x*c->point[v].y - c->point[lv].y*c->point[v].x + c->point[lv].y*c->point[nv].x - c->point[lv].x*c->point[nv].y + c->point[v].x*c->point[nv].y - c->point[nv].x*c->point[v].y;
129 if(area > 0)
130 return(1); /* counter clockwise orientation */
131 else if(area < 0) /* clockwise orientation */
132 return(-1);
133 else
134 return(0); /* shouldn't happen unless the polygon is self intersecting */
135 }
136
137 /*
138 ** Copyright (c) 1970-2003, Wm. Randolph Franklin
139 **
140 ** Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
141 ** associated documentation files (the "Software"), to deal in the Software without restriction, including
142 ** without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
143 ** copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the
144 ** following conditions:
145 **
146 ** 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the
147 ** following disclaimers.
148 ** 2. Redistributions in binary form must reproduce the above copyright notice in the documentation and/or
149 ** other materials provided with the distribution.
150 ** 3. The name of W. Randolph Franklin may not be used to endorse or promote products derived from this
151 ** Software without specific prior written permission.
152 **
153 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT
154 ** LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
155 ** NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
156 ** WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
157 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
158 */
msPointInPolygon(pointObj * p,lineObj * c)159 int msPointInPolygon(pointObj *p, lineObj *c)
160 {
161 int i, j, status = MS_FALSE;
162
163 for (i = 0, j = c->numpoints-1; i < c->numpoints; j = i++) {
164 if ((((c->point[i].y<=p->y) && (p->y<c->point[j].y)) || ((c->point[j].y<=p->y) && (p->y<c->point[i].y))) && (p->x < (c->point[j].x - c->point[i].x) * (p->y - c->point[i].y) / (c->point[j].y - c->point[i].y) + c->point[i].x))
165 status = !status;
166 }
167 return status;
168 }
169
170 /*
171 ** Note: the following functions are brute force implementations. Some fancy
172 ** computational geometry algorithms would speed things up nicely in many
173 ** cases. In due time... -SDL-
174 */
175
msIntersectSegments(const pointObj * a,const pointObj * b,const pointObj * c,const pointObj * d)176 int msIntersectSegments(const pointObj *a, const pointObj *b, const pointObj *c, const pointObj *d) /* from comp.graphics.alogorithms FAQ */
177 {
178
179 double r, s;
180 double denominator, numerator;
181
182 numerator = ((a->y-c->y)*(d->x-c->x) - (a->x-c->x)*(d->y-c->y));
183 denominator = ((b->x-a->x)*(d->y-c->y) - (b->y-a->y)*(d->x-c->x));
184
185 if((denominator == 0) && (numerator == 0)) { /* lines are coincident, intersection is a line segement if it exists */
186 if(a->y == c->y) { /* coincident horizontally, check x's */
187 if(((a->x >= MS_MIN(c->x,d->x)) && (a->x <= MS_MAX(c->x,d->x))) || ((b->x >= MS_MIN(c->x,d->x)) && (b->x <= MS_MAX(c->x,d->x))))
188 return(MS_TRUE);
189 else
190 return(MS_FALSE);
191 } else { /* test for y's will work fine for remaining cases */
192 if(((a->y >= MS_MIN(c->y,d->y)) && (a->y <= MS_MAX(c->y,d->y))) || ((b->y >= MS_MIN(c->y,d->y)) && (b->y <= MS_MAX(c->y,d->y))))
193 return(MS_TRUE);
194 else
195 return(MS_FALSE);
196 }
197 }
198
199 if(denominator == 0) /* lines are parallel, can't intersect */
200 return(MS_FALSE);
201
202 r = numerator/denominator;
203
204 if((r<0) || (r>1))
205 return(MS_FALSE); /* no intersection */
206
207 numerator = ((a->y-c->y)*(b->x-a->x) - (a->x-c->x)*(b->y-a->y));
208 s = numerator/denominator;
209
210 if((s<0) || (s>1))
211 return(MS_FALSE); /* no intersection */
212
213 return(MS_TRUE);
214 }
215
216 /*
217 ** Instead of using ring orientation we count the number of parts the
218 ** point falls in. If odd the point is in the polygon, if 0 or even
219 ** then the point is in a hole or completely outside.
220 */
msIntersectPointPolygon(pointObj * point,shapeObj * poly)221 int msIntersectPointPolygon(pointObj *point, shapeObj *poly)
222 {
223 int i;
224 int status=MS_FALSE;
225
226 for(i=0; i<poly->numlines; i++) {
227 if(msPointInPolygon(point, &poly->line[i]) == MS_TRUE) /* ok, the point is in a line */
228 status = !status;
229 }
230
231 return(status);
232 }
233
msIntersectMultipointPolygon(shapeObj * multipoint,shapeObj * poly)234 int msIntersectMultipointPolygon(shapeObj *multipoint, shapeObj *poly)
235 {
236 int i,j;
237
238 /* The change to loop through all the lines has been made for ticket
239 * #2443 but is no more needed since ticket #2762. PostGIS now put all
240 * points into a single line. */
241 for(i=0; i<multipoint->numlines; i++ ) {
242 lineObj points = multipoint->line[i];
243 for(j=0; j<points.numpoints; j++) {
244 if(msIntersectPointPolygon(&(points.point[j]), poly) == MS_TRUE)
245 return(MS_TRUE);
246 }
247 }
248
249 return(MS_FALSE);
250 }
251
msIntersectPolylines(shapeObj * line1,shapeObj * line2)252 int msIntersectPolylines(shapeObj *line1, shapeObj *line2)
253 {
254 int c1,v1,c2,v2;
255
256 for(c1=0; c1<line1->numlines; c1++)
257 for(v1=1; v1<line1->line[c1].numpoints; v1++)
258 for(c2=0; c2<line2->numlines; c2++)
259 for(v2=1; v2<line2->line[c2].numpoints; v2++)
260 if(msIntersectSegments(&(line1->line[c1].point[v1-1]), &(line1->line[c1].point[v1]),
261 &(line2->line[c2].point[v2-1]), &(line2->line[c2].point[v2])) == MS_TRUE)
262 return(MS_TRUE);
263
264 return(MS_FALSE);
265 }
266
msIntersectPolylinePolygon(shapeObj * line,shapeObj * poly)267 int msIntersectPolylinePolygon(shapeObj *line, shapeObj *poly)
268 {
269 int i;
270
271 /* STEP 1: polygon might competely contain the polyline or one of it's parts (only need to check one point from each part) */
272 for(i=0; i<line->numlines; i++) {
273 if(msIntersectPointPolygon(&(line->line[i].point[0]), poly) == MS_TRUE) /* this considers holes and multiple parts */
274 return(MS_TRUE);
275 }
276
277 /* STEP 2: look for intersecting line segments */
278 if (msIntersectPolylines(line, poly) == MS_TRUE)
279 return (MS_TRUE);
280
281 return(MS_FALSE);
282 }
283
msIntersectPolygons(shapeObj * p1,shapeObj * p2)284 int msIntersectPolygons(shapeObj *p1, shapeObj *p2)
285 {
286 int i;
287
288 /* STEP 1: polygon 1 completely contains 2 (only need to check one point from each part) */
289 for(i=0; i<p2->numlines; i++) {
290 if(msIntersectPointPolygon(&(p2->line[i].point[0]), p1) == MS_TRUE) /* this considers holes and multiple parts */
291 return(MS_TRUE);
292 }
293
294 /* STEP 2: polygon 2 completely contains 1 (only need to check one point from each part) */
295 for(i=0; i<p1->numlines; i++) {
296 if(msIntersectPointPolygon(&(p1->line[i].point[0]), p2) == MS_TRUE) /* this considers holes and multiple parts */
297 return(MS_TRUE);
298 }
299
300 /* STEP 3: look for intersecting line segments */
301 if (msIntersectPolylines(p1, p2) == MS_TRUE)
302 return(MS_TRUE);
303
304 /*
305 ** At this point we know there are are no intersections between edges. There may be other tests necessary
306 ** but I haven't run into any cases that require them.
307 */
308
309 return(MS_FALSE);
310 }
311
312
313 /*
314 ** Distance computations
315 */
316
msDistancePointToPoint(pointObj * a,pointObj * b)317 double msDistancePointToPoint(pointObj *a, pointObj *b)
318 {
319 double d;
320
321 d = sqrt(msSquareDistancePointToPoint(a, b));
322
323 return(d);
324 }
325
326 /*
327 ** Quickly compute the square of the distance; avoids expensive sqrt() call on each invocation
328 */
msSquareDistancePointToPoint(pointObj * a,pointObj * b)329 double msSquareDistancePointToPoint(pointObj *a, pointObj *b)
330 {
331 double dx, dy;
332
333 dx = a->x - b->x;
334 dy = a->y - b->y;
335
336 return(dx*dx + dy*dy);
337 }
338
msDistancePointToSegment(pointObj * p,pointObj * a,pointObj * b)339 double msDistancePointToSegment(pointObj *p, pointObj *a, pointObj *b)
340 {
341 return (sqrt(msSquareDistancePointToSegment(p, a, b)));
342 }
343
msSquareDistancePointToSegment(pointObj * p,pointObj * a,pointObj * b)344 double msSquareDistancePointToSegment(pointObj *p, pointObj *a, pointObj *b)
345 {
346 double l_squared; /* squared length of line ab */
347 double r,s;
348
349 l_squared = msSquareDistancePointToPoint(a,b);
350
351 if(l_squared == 0.0) /* a = b */
352 return(msSquareDistancePointToPoint(a,p));
353
354 r = ((a->y - p->y)*(a->y - b->y) - (a->x - p->x)*(b->x - a->x))/(l_squared);
355
356 if(r > 1) /* perpendicular projection of P is on the forward extention of AB */
357 return(MS_MIN(msSquareDistancePointToPoint(p, b),msSquareDistancePointToPoint(p, a)));
358 if(r < 0) /* perpendicular projection of P is on the backward extention of AB */
359 return(MS_MIN(msSquareDistancePointToPoint(p, b),msSquareDistancePointToPoint(p, a)));
360
361 s = ((a->y - p->y)*(b->x - a->x) - (a->x - p->x)*(b->y - a->y))/l_squared;
362
363 return(fabs(s*s*l_squared));
364 }
365
366 #define SMALL_NUMBER 0.00000001
367 #define dot(u,v) ((u).x *(v).x + (u).y *(v).y) /* vector dot product */
368 #define norm(v) sqrt(dot(v,v))
369
370 #define slope(a,b) (((a)->y - (b)->y)/((a)->x - (b)->x))
371
372 /* Segment to segment distance code is a modified version of that found at: */
373 /* */
374 /* http://www.geometryalgorithms.com/Archive/algorithm_0106/algorithm_0106.htm */
375 /* */
376 /* Copyright 2001, softSurfer (www.softsurfer.com) */
377 /* This code may be freely used and modified for any purpose */
378 /* providing that this copyright notice is included with it. */
379 /* SoftSurfer makes no warranty for this code, and cannot be held */
380 /* liable for any real or imagined damage resulting from its use. */
381 /* Users of this code must verify correctness for their application. */
382
msDistanceSegmentToSegment(pointObj * pa,pointObj * pb,pointObj * pc,pointObj * pd)383 double msDistanceSegmentToSegment(pointObj *pa, pointObj *pb, pointObj *pc, pointObj *pd)
384 {
385 vectorObj dP;
386 vectorObj u, v, w;
387 double a, b, c, d, e;
388 double D;
389 double sc, sN, sD; /* N=numerator, D=demoninator */
390 double tc, tN, tD;
391
392 /* check for strictly parallel segments first */
393 /* if(((pa->x == pb->x) && (pc->x == pd->x)) || (slope(pa,pb) == slope(pc,pd))) { // vertical (infinite slope) || otherwise parallel */
394 /* D = msDistancePointToSegment(pa, pc, pd); */
395 /* D = MS_MIN(D, msDistancePointToSegment(pb, pc, pd)); */
396 /* D = MS_MIN(D, msDistancePointToSegment(pc, pa, pb)); */
397 /* return(MS_MIN(D, msDistancePointToSegment(pd, pa, pb))); */
398 /* } */
399
400 u.x = pb->x - pa->x; /* u = pb - pa */
401 u.y = pb->y - pa->y;
402 v.x = pd->x - pc->x; /* v = pd - pc */
403 v.y = pd->y - pc->y;
404 w.x = pa->x - pc->x; /* w = pa - pc */
405 w.y = pa->y - pc->y;
406
407 a = dot(u,u);
408 b = dot(u,v);
409 c = dot(v,v);
410 d = dot(u,w);
411 e = dot(v,w);
412
413 D = a*c - b*b;
414 sc = sN = sD = D;
415 tc = tN = tD = D;
416
417 /* compute the line parameters of the two closest points */
418 if(D < SMALL_NUMBER) { /* lines are parallel or almost parallel */
419 sN = 0.0;
420 sD = 1.0;
421 tN = e;
422 tD = c;
423 } else { /* get the closest points on the infinite lines */
424 sN = b*e - c*d;
425 tN = a*e - b*d;
426 if(sN < 0) {
427 sN = 0.0;
428 tN = e;
429 tD = c;
430 } else if(sN > sD) {
431 sN = sD;
432 tN = e + b;
433 tD = c;
434 }
435 }
436
437 if(tN < 0) {
438 tN = 0.0;
439 if(-d < 0)
440 sN = 0.0;
441 else if(-d > a)
442 sN = sD;
443 else {
444 sN = -d;
445 sD = a;
446 }
447 } else if(tN > tD) {
448 tN = tD;
449 if((-d + b) < 0)
450 sN = 0.0;
451 else if((-d + b) > a)
452 sN = sD;
453 else {
454 sN = (-d + b);
455 sD = a;
456 }
457 }
458
459 /* finally do the division to get sc and tc */
460 sc = sN/sD;
461 tc = tN/tD;
462
463 dP.x = w.x + (sc*u.x) - (tc*v.x);
464 dP.y = w.y + (sc*u.y) - (tc*v.y);
465
466 return(norm(dP));
467 }
468
msDistancePointToShape(pointObj * point,shapeObj * shape)469 double msDistancePointToShape(pointObj *point, shapeObj *shape)
470 {
471 double d;
472
473 d = msSquareDistancePointToShape(point, shape);
474
475 return(sqrt(d));
476 }
477
478 /*
479 ** As msDistancePointToShape; avoid expensive sqrt calls
480 */
msSquareDistancePointToShape(pointObj * point,shapeObj * shape)481 double msSquareDistancePointToShape(pointObj *point, shapeObj *shape)
482 {
483 int i, j;
484 double dist, minDist=-1;
485
486 switch(shape->type) {
487 case(MS_SHAPE_POINT):
488 for(j=0; j<shape->numlines; j++) {
489 for(i=0; i<shape->line[j].numpoints; i++) {
490 dist = msSquareDistancePointToPoint(point, &(shape->line[j].point[i]));
491 if((dist < minDist) || (minDist < 0)) minDist = dist;
492 }
493 }
494 break;
495 case(MS_SHAPE_LINE):
496 for(j=0; j<shape->numlines; j++) {
497 for(i=1; i<shape->line[j].numpoints; i++) {
498 dist = msSquareDistancePointToSegment(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
499 if((dist < minDist) || (minDist < 0)) minDist = dist;
500 }
501 }
502 break;
503 case(MS_SHAPE_POLYGON):
504 if(msIntersectPointPolygon(point, shape))
505 minDist = 0; /* point is IN the shape */
506 else { /* treat shape just like a line */
507 for(j=0; j<shape->numlines; j++) {
508 for(i=1; i<shape->line[j].numpoints; i++) {
509 dist = msSquareDistancePointToSegment(point, &(shape->line[j].point[i-1]), &(shape->line[j].point[i]));
510 if((dist < minDist) || (minDist < 0)) minDist = dist;
511 }
512 }
513 }
514 break;
515 default:
516 break;
517 }
518
519 return(minDist);
520 }
521
msDistanceShapeToShape(shapeObj * shape1,shapeObj * shape2)522 double msDistanceShapeToShape(shapeObj *shape1, shapeObj *shape2)
523 {
524 int i,j,k,l;
525 double dist, minDist=-1;
526
527 switch(shape1->type) {
528 case(MS_SHAPE_POINT): /* shape1 */
529 for(i=0; i<shape1->numlines; i++) {
530 for(j=0; j<shape1->line[i].numpoints; j++) {
531 dist = msSquareDistancePointToShape(&(shape1->line[i].point[j]), shape2);
532 if((dist < minDist) || (minDist < 0))
533 minDist = dist;
534 }
535 }
536 minDist = sqrt(minDist);
537 break;
538 case(MS_SHAPE_LINE): /* shape1 */
539 switch(shape2->type) {
540 case(MS_SHAPE_POINT):
541 for(i=0; i<shape2->numlines; i++) {
542 for(j=0; j<shape2->line[i].numpoints; j++) {
543 dist = msSquareDistancePointToShape(&(shape2->line[i].point[j]), shape1);
544 if((dist < minDist) || (minDist < 0))
545 minDist = dist;
546 }
547 }
548 minDist = sqrt(minDist);
549 break;
550 case(MS_SHAPE_LINE):
551 for(i=0; i<shape1->numlines; i++) {
552 for(j=1; j<shape1->line[i].numpoints; j++) {
553 for(k=0; k<shape2->numlines; k++) {
554 for(l=1; l<shape2->line[k].numpoints; l++) {
555 /* check intersection (i.e. dist=0) */
556 if(msIntersectSegments(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l])) == MS_TRUE)
557 return(0);
558
559 /* no intersection, compute distance */
560 dist = msDistanceSegmentToSegment(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l]));
561 if((dist < minDist) || (minDist < 0))
562 minDist = dist;
563 }
564 }
565 }
566 }
567 break;
568 case(MS_SHAPE_POLYGON):
569 /* shape2 (the polygon) could contain shape1 or one of it's parts */
570 for(i=0; i<shape1->numlines; i++) {
571 if(msIntersectPointPolygon(&(shape1->line[0].point[0]), shape2) == MS_TRUE) /* this considers holes and multiple parts */
572 return(0);
573 }
574
575 /* check segment intersection and, if necessary, distance between segments */
576 for(i=0; i<shape1->numlines; i++) {
577 for(j=1; j<shape1->line[i].numpoints; j++) {
578 for(k=0; k<shape2->numlines; k++) {
579 for(l=1; l<shape2->line[k].numpoints; l++) {
580 /* check intersection (i.e. dist=0) */
581 if(msIntersectSegments(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l])) == MS_TRUE)
582 return(0);
583
584 /* no intersection, compute distance */
585 dist = msDistanceSegmentToSegment(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l]));
586 if((dist < minDist) || (minDist < 0))
587 minDist = dist;
588 }
589 }
590 }
591 }
592 break;
593 }
594 break;
595 case(MS_SHAPE_POLYGON): /* shape1 */
596 switch(shape2->type) {
597 case(MS_SHAPE_POINT):
598 for(i=0; i<shape2->numlines; i++) {
599 for(j=0; j<shape2->line[i].numpoints; j++) {
600 dist = msSquareDistancePointToShape(&(shape2->line[i].point[j]), shape1);
601 if((dist < minDist) || (minDist < 0))
602 minDist = dist;
603 }
604 }
605 minDist = sqrt(minDist);
606 break;
607 case(MS_SHAPE_LINE):
608 /* shape1 (the polygon) could contain shape2 or one of it's parts */
609 for(i=0; i<shape2->numlines; i++) {
610 if(msIntersectPointPolygon(&(shape2->line[i].point[0]), shape1) == MS_TRUE) /* this considers holes and multiple parts */
611 return(0);
612 }
613
614 /* check segment intersection and, if necessary, distance between segments */
615 for(i=0; i<shape1->numlines; i++) {
616 for(j=1; j<shape1->line[i].numpoints; j++) {
617 for(k=0; k<shape2->numlines; k++) {
618 for(l=1; l<shape2->line[k].numpoints; l++) {
619 /* check intersection (i.e. dist=0) */
620 if(msIntersectSegments(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l])) == MS_TRUE)
621 return(0);
622
623 /* no intersection, compute distance */
624 dist = msDistanceSegmentToSegment(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l]));
625 if((dist < minDist) || (minDist < 0))
626 minDist = dist;
627 }
628 }
629 }
630 }
631 break;
632 case(MS_SHAPE_POLYGON):
633 /* shape1 completely contains shape2 (only need to check one point from each part) */
634 for(i=0; i<shape2->numlines; i++) {
635 if(msIntersectPointPolygon(&(shape2->line[i].point[0]), shape1) == MS_TRUE) /* this considers holes and multiple parts */
636 return(0);
637 }
638
639 /* shape2 completely contains shape1 (only need to check one point from each part) */
640 for(i=0; i<shape1->numlines; i++) {
641 if(msIntersectPointPolygon(&(shape1->line[i].point[0]), shape2) == MS_TRUE) /* this considers holes and multiple parts */
642 return(0);
643 }
644
645 /* check segment intersection and, if necessary, distance between segments */
646 for(i=0; i<shape1->numlines; i++) {
647 for(j=1; j<shape1->line[i].numpoints; j++) {
648 for(k=0; k<shape2->numlines; k++) {
649 for(l=1; l<shape2->line[k].numpoints; l++) {
650 /* check intersection (i.e. dist=0) */
651 if(msIntersectSegments(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l])) == MS_TRUE)
652 return(0);
653
654 /* no intersection, compute distance */
655 dist = msDistanceSegmentToSegment(&(shape1->line[i].point[j-1]), &(shape1->line[i].point[j]), &(shape2->line[k].point[l-1]), &(shape2->line[k].point[l]));
656 if((dist < minDist) || (minDist < 0))
657 minDist = dist;
658 }
659 }
660 }
661 }
662 break;
663 }
664 break;
665 }
666
667 return(minDist);
668 }
669