Ticket #63: max_distance_new_poly_poly.patch
| File max_distance_new_poly_poly.patch, 39.6 KB (added by nicklas, 3 years ago) |
|---|
-
doc/reference.xml
9833 9833 </refsection> 9834 9834 </refentry> 9835 9835 9836 <refentry id="ST_Max_Distance"> 9837 <refnamediv> 9838 <refname>ST_Max_Distance</refname> 9839 9840 <refpurpose>Returns the 2-dimensional largest distance between two geometries in 9841 projected units.</refpurpose> 9842 </refnamediv> 9843 9844 <refsynopsisdiv> 9845 <funcsynopsis> 9846 <funcprototype> 9847 <funcdef>float <function>ST_Max_Distance</function></funcdef> 9848 9849 <paramdef><type>geometry </type> 9850 <parameter>g1</parameter></paramdef> 9851 9852 <paramdef><type>geometry </type> 9853 <parameter>g2</parameter></paramdef> 9854 </funcprototype> 9855 </funcsynopsis> 9856 </refsynopsisdiv> 9857 9858 <refsection> 9859 <title>Description</title> 9860 9861 <para>Returns the 2-dimensional maximum cartesian distance between two linestrings in 9862 projected units. If g1 and g2 is the same geometry the function will return the distance between 9863 the two vertices most far from eachother in that geometry.</para> 9864 9865 </refsection> 9866 9867 <refsection> 9868 <title>Examples</title> 9869 9870 <programlisting>postgis=# SELECT ST_Max_Distance('POINT(0 0)'::geometry, 'LINESTRING ( 2 0, 0 2 )'::geometry); 9871 st_max_distance 9872 ----------------- 9873 2 9874 (1 row)</programlisting> 9875 </refsection> 9876 9877 <refsection> 9878 <title>See Also</title> 9879 9880 <para><xref linkend="ST_Distance"/></para> 9881 </refsection> 9882 </refentry> 9883 9836 9884 <refentry id="ST_Distance_Sphere"> 9837 9885 <refnamediv> 9838 9886 <refname>ST_Distance_Sphere</refname> -
liblwgeom/liblwgeom.h
1 1 2 /********************************************************************** 2 3 * $Id$ 3 4 * … … 19 20 #include <stdarg.h> 20 21 #include <stdio.h> 21 22 23 22 24 /** 23 25 * @file liblwgeom.h 24 26 * … … 290 292 } 291 293 POINT4D; 292 294 295 typedef struct 296 { 297 double d; /*the distance between p1 and p2*/ 298 POINT2D p1; 299 POINT2D p2; 300 int thedir; /*the direction of looking, if thedir = -1 then we look for maxdistance and if it is 1 then we look for mindistance*/ 301 int twisted; /*To preserve the order of incoming points to match the first and secon point in shortest and longest line*/ 302 double tolerance; /*the tolerance for dwithin and dfullywithin*/ 303 } DISTPTS; 293 304 /******************************************************************/ 294 305 295 306 /* … … 1184 1195 extern void lwgeom_force3dz_recursive(uchar *serialized, uchar *optr, size_t *retsize); 1185 1196 extern void lwgeom_force3dm_recursive(uchar *serialized, uchar *optr, size_t *retsize); 1186 1197 extern void lwgeom_force4d_recursive(uchar *serialized, uchar *optr, size_t *retsize); 1198 1187 1199 extern double distance2d_pt_pt(POINT2D *p1, POINT2D *p2); 1200 extern void lw_dist2d_comp_pt_pt(POINT2D *p1, POINT2D *p2, DISTPTS *dl); 1188 1201 extern double distance2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B); 1189 extern double distance2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D); 1190 extern double distance2d_pt_ptarray(POINT2D *p, POINTARRAY *pa); 1191 extern double distance2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2); 1202 extern void lw_dist2d_comp_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl); 1203 extern void lw_dist2d_comp_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl); 1204 extern void lw_dist2d_comp_pt_ptarray(POINT2D *p, POINTARRAY *pa, DISTPTS *dl); 1205 extern void lw_dist2d_comp_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl); 1192 1206 extern int pt_in_ring_2d(POINT2D *p, POINTARRAY *ring); 1193 1207 extern int pt_in_poly_2d(POINT2D *p, LWPOLY *poly); 1194 extern double distance2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly);1195 extern double distance2d_point_point(LWPOINT *point1, LWPOINT *point2);1196 extern double distance2d_point_line(LWPOINT *point, LWLINE *line);1197 extern double distance2d_line_line(LWLINE *line1, LWLINE *line2);1198 extern double distance2d_point_poly(LWPOINT *point, LWPOLY *poly);1199 extern double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2);1200 extern double distance2d_line_poly(LWLINE *line, LWPOLY *poly);1208 extern void lw_dist2d_comp_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl); 1209 extern void lw_dist2d_comp_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl); 1210 extern void lw_dist2d_comp_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl); 1211 extern void lw_dist2d_comp_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl); 1212 extern void lw_dist2d_comp_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl); 1213 extern void lw_dist2d_comp_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl); 1214 extern void lw_dist2d_comp_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl); 1201 1215 extern int azimuth_pt_pt(POINT2D *p1, POINT2D *p2, double *ret); 1202 1216 extern double lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2); 1203 1217 extern double lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance); 1218 extern double lwgeom_maxdistance2d_recursive(uchar *lw1, uchar *lw2); 1219 extern double lwgeom_maxdistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance); 1220 extern void lw_dist2d_comp(uchar *lw1, uchar *lw2, DISTPTS *dl); 1204 1221 extern int lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad); 1205 1222 extern int32 lwgeom_npoints(uchar *serialized); 1206 1223 extern char ptarray_isccw(const POINTARRAY *pa); -
liblwgeom/measures.c
81 81 return (cn&1); /* 0 if even (out), and 1 if odd (in) */ 82 82 } 83 83 84 /*Compares incomming points and stores the points closest to each other or most far away from each other depending on dl->thedir (thedirection) */ 85 void lw_dist2d_comp_pt_pt(POINT2D *thep1, POINT2D *thep2,DISTPTS *dl) 86 { 87 double hside = thep2->x - thep1->x; 88 double vside = thep2->y - thep1->y; 89 double dist = sqrt ( hside*hside + vside*vside ); 90 if (((dl->d - dist)*(dl->thedir))>0) /*multiplication with thedir to handle mindistance (thedir=1) and maxdistance (thedir = (-1)*/ 91 { 92 dl->d = dist; 93 94 if (dl->twisted>0) /*To get the points in right order. twisted is updated between 1 and (-1) every time the order is changed earlier in the chain*/ 95 { 96 dl->p1 = *thep1; 97 dl->p2 = *thep2; 98 } 99 else 100 { 101 dl->p1 = *thep2; 102 dl->p2 = *thep1; 103 } 104 } 105 return; 106 } 107 108 /*The old function nessecary for ptarray_segmentize2d in ptarray.c*/ 84 109 double distance2d_pt_pt(POINT2D *p1, POINT2D *p2) 85 110 { 86 111 double hside = p2->x - p1->x; … … 94 119 ); */ 95 120 } 96 121 97 /*distance2d from p to line A->B */ 122 /*lw_dist2d_comp from p to line A->B 123 This one is now sending every occation to lw_dist2d_comp_pt_pt 124 Berfore it was handling occations where r was between 0 and 1 internally and just returning the distance without identifying the points. 125 To get this points it was nessecary to change and it also showed to be about 10%faster. */ 126 void lw_dist2d_comp_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl) 127 { 128 /*lwnotice("lw_dist2d_comp_pt_seg %e",dl->d);*/ 129 double r; 130 131 /*if start==end, then use pt distance */ 132 if ( ( A->x == B->x) && (A->y == B->y) ) 133 { 134 lw_dist2d_comp_pt_pt(p,A,dl); 135 return; 136 } 137 /* 138 * otherwise, we use comp.graphics.algorithms 139 * Frequently Asked Questions method 140 * 141 * (1) AC dot AB 142 * r = --------- 143 * ||AB||^2 144 * r has the following meaning: 145 * r=0 P = A 146 * r=1 P = B 147 * r<0 P is on the backward extension of AB 148 * r>1 P is on the forward extension of AB 149 * 0<r<1 P is interior to AB 150 */ 151 152 r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) )/( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) ); 153 154 /*This is for finding the maxdistance. 155 the maxdistance have to be between two vertexes as I understand it, 156 compared to mindistance which can be between some point between the vertices and a vertex.*/ 157 if (dl->thedir < 1) 158 { 159 if (r>=0.5) 160 { 161 lw_dist2d_comp_pt_pt(p,A,dl); 162 return; 163 } 164 if (r<0.5) 165 { 166 lw_dist2d_comp_pt_pt(p,B,dl); 167 return; 168 } 169 } 170 171 if (r<0) /*If the first vertex A is closest to the point p*/ 172 { 173 lw_dist2d_comp_pt_pt(p,A,dl); 174 return; 175 } 176 if (r>1) /*If the second vertex B is closest to the point p*/ 177 { 178 179 lw_dist2d_comp_pt_pt(p,B,dl); 180 return; 181 } 182 183 184 /*else if the point p is closer to some point between a and b then we find that point and send it to lw_dist2d_comp_pt_pt*/ 185 POINT2D c; 186 c.x=A->x + r * (B->x-A->x); 187 c.y=A->y + r * (B->y-A->y); 188 189 lw_dist2d_comp_pt_pt(p,&c,dl); 190 return; 191 } 192 193 194 195 196 /*The old function nessecary for ptarray_segmentize2d in ptarray.c*/ 98 197 double distance2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B) 99 198 { 100 199 double r,s; … … 142 241 ); 143 242 } 144 243 145 /* find the minimum 2d distance from AB to CD */ 146 double distance2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D) 244 245 246 247 248 /*This function is changed so it is not doing any comparasion of distance 249 but just sending every possible combination further to lw_dist2d_comp_pt_seg*/ 250 void lw_dist2d_comp_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl) 147 251 { 148 252 149 253 double s_top, s_bot,s; 150 254 double r_top, r_bot,r; 151 255 152 LWDEBUGF(2, "distance2d_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]", 256 257 LWDEBUGF(2, "lw_dist2d_comp_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]", 153 258 A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y); 154 259 155 156 260 /*A and B are the same point */ 157 261 if ( ( A->x == B->x) && (A->y == B->y) ) 158 return distance2d_pt_seg(A,C,D); 159 262 { 263 lw_dist2d_comp_pt_seg(A,C,D,dl); 264 return; 265 } 160 266 /*U and V are the same point */ 161 267 162 268 if ( ( C->x == D->x) && (C->y == D->y) ) 163 return distance2d_pt_seg(D,A,B); 164 269 { 270 dl->twisted= ((dl->twisted) * (-1)); 271 lw_dist2d_comp_pt_seg(D,A,B,dl); 272 return; 273 } 165 274 /* AB and CD are line segments */ 166 275 /* from comp.graphics.algo 167 276 … … 192 301 193 302 if ( (r_bot==0) || (s_bot == 0) ) 194 303 { 195 return ( 196 LW_MIN(distance2d_pt_seg(A,C,D), 197 LW_MIN(distance2d_pt_seg(B,C,D), 198 LW_MIN(distance2d_pt_seg(C,A,B), 199 distance2d_pt_seg(D,A,B)) 200 ) 201 ) 202 ); 304 lw_dist2d_comp_pt_seg(A,C,D,dl); 305 lw_dist2d_comp_pt_seg(B,C,D,dl); 306 307 dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/ 308 lw_dist2d_comp_pt_seg(C,A,B,dl); 309 lw_dist2d_comp_pt_seg(D,A,B,dl); 310 return; 203 311 } 312 204 313 s = s_top/s_bot; 205 314 r= r_top/r_bot; 206 315 207 if (( r<0) || (r>1) || (s<0) || (s>1))316 if (((r<0) || (r>1) || (s<0) || (s>1)) || (dl->thedir <0)) 208 317 { 209 /*no intersection */ 210 return ( 211 LW_MIN(distance2d_pt_seg(A,C,D), 212 LW_MIN(distance2d_pt_seg(B,C,D), 213 LW_MIN(distance2d_pt_seg(C,A,B), 214 distance2d_pt_seg(D,A,B)) 215 ) 216 ) 217 ); 218 318 lw_dist2d_comp_pt_seg(A,C,D,dl); 319 lw_dist2d_comp_pt_seg(B,C,D,dl); 320 321 dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/ 322 lw_dist2d_comp_pt_seg(C,A,B,dl); 323 lw_dist2d_comp_pt_seg(D,A,B,dl); 324 return; 219 325 } 220 326 else 221 return -0; /*intersection exists */ 327 { 328 if (dl->thedir >0) /*If there is intersection we identify the intersection point and return it but only if we are looking for mindistance*/ 329 { /*lwnotice("intersection lines: %e %e %e %e",A->x,A->y,B->x,B->y);*/ 330 POINT2D theP; 331 332 if(((A->x==C->x)&&(A->y==C->y))||((A->x==D->x)&&(A->y==D->y))) 333 { 334 theP.x = A->x; 335 theP.y = A->y; 336 } 337 else if(((B->x==C->x)&&(B->y==C->y))||((B->x==D->x)&&(B->y==D->y))) 338 { 339 theP.x = B->x; 340 theP.y = B->y; 341 } 342 else 343 { 344 theP.x = A->x+r*(B->x-A->x); 345 theP.y = A->y+r*(B->y-A->y); 346 } 347 /*lwnotice("intersection: %e %e ",theP.x, theP.y);*/ 348 dl->d=0.0; 349 dl->p1=theP; 350 dl->p2=theP; 351 } 352 return; 222 353 354 } 355 223 356 } 224 357 225 358 /* 226 359 * search all the segments of pointarray to see which one is closest to p1 227 360 * Returns minimum distance between point and pointarray 228 361 */ 229 double distance2d_pt_ptarray(POINT2D *p, POINTARRAY *pa)362 void lw_dist2d_comp_pt_ptarray(POINT2D *p, POINTARRAY *pa,DISTPTS *dl) 230 363 { 231 double result = 0; 364 365 232 366 int t; 233 367 POINT2D start, end; 234 368 int twist = dl->twisted; 235 369 getPoint2d_p(pa, 0, &start); 236 370 237 371 for (t=1; t<pa->npoints; t++) 238 372 { 239 d ouble dist;373 dl->twisted=twist; 240 374 getPoint2d_p(pa, t, &end); 241 dist = distance2d_pt_seg(p, &start, &end); 242 if (t==1) result = dist; 243 else result = LW_MIN(result, dist); 244 245 if ( result == 0 ) return 0; 246 375 lw_dist2d_comp_pt_seg(p, &start, &end,dl); 376 377 if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ 378 247 379 start = end; 248 380 } 249 381 250 return result;382 return; 251 383 } 252 384 253 385 /* test each segment of l1 against each segment of l2. Return min */ 254 double distance2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2)386 void lw_dist2d_comp_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl) 255 387 { 256 double result = 99999999999.9; 257 char result_okay = 0; /*result is a valid min */ 388 /*lwnotice("lw_dist2d_comp_ptarray_ptarray");*/ 258 389 int t,u; 259 390 POINT2D start, end; 260 391 POINT2D start2, end2; 261 262 LWDEBUGF(2, " distance2d_ptarray_ptarray called (points: %d-%d)",392 int twist = dl->twisted; 393 LWDEBUGF(2, "lw_dist2d_comp_ptarray_ptarray called (points: %d-%d)", 263 394 l1->npoints, l2->npoints); 264 395 265 396 getPoint2d_p(l1, 0, &start); 266 397 for (t=1; t<l1->npoints; t++) /*for each segment in L1 */ 267 398 { 268 399 getPoint2d_p(l1, t, &end); 269 270 400 getPoint2d_p(l2, 0, &start2); 271 401 for (u=1; u<l2->npoints; u++) /*for each segment in L2 */ 272 402 { 273 double dist;274 275 403 getPoint2d_p(l2, u, &end2); 404 dl->twisted=twist; 405 lw_dist2d_comp_seg_seg(&start, &end, &start2, &end2,dl); 276 406 277 dist = distance2d_seg_seg(&start, &end, &start2, &end2);278 407 408 279 409 LWDEBUGF(4, "line_line; seg %i * seg %i, dist = %g\n",t,u,dist); 280 281 if (result_okay)282 result = LW_MIN(result,dist);283 else284 {285 result_okay = 1;286 result = dist;287 }288 289 410 LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f", 290 411 t, u, dist, result); 291 292 if (result <= 0) return 0; /*intersection */ 293 412 if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ 294 413 start2 = end2; 295 414 } 296 415 start = end; 297 416 } 298 299 return result; 417 return ; 300 418 } 301 419 420 302 421 /* true if point is in poly (and not in its holes) */ 303 422 int pt_in_poly_2d(POINT2D *p, LWPOLY *poly) 304 423 { … … 317 436 return 1; /* In outer ring, not in holes */ 318 437 } 319 438 439 440 320 441 /* 321 442 * Brute force. 322 443 * Test line-ring distance against each ring. … … 327 448 * otherwise return min distance to a ring (could be outside 328 449 * polygon or inside a hole) 329 450 */ 330 double distance2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly) 451 452 void lw_dist2d_comp_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl) 331 453 { 332 454 POINT2D pt; 333 455 int i; 334 double mindist = 0;335 456 336 LWDEBUGF(2, "distance2d_ptarray_poly called (%d rings)", poly->nrings);337 457 458 LWDEBUGF(2, "lw_dist2d_comp_ptarray_poly called (%d rings)", poly->nrings); 459 460 338 461 for (i=0; i<poly->nrings; i++) 339 462 { 340 double dist = distance2d_ptarray_ptarray(pa, poly->rings[i]); 341 if (i) mindist = LW_MIN(mindist, dist); 342 else mindist = dist; 463 lw_dist2d_comp_ptarray_ptarray(pa, poly->rings[i], dl); 343 464 344 465 LWDEBUGF(3, " distance from ring %d: %f, mindist: %f", 345 466 i, dist, mindist); 346 347 if ( mindist <= 0 ) return 0.0; /* intersection */ 467 if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ 348 468 } 349 469 350 470 /* … … 357 477 * Outside outer ring, so min distance to a ring 358 478 * is the actual min distance 359 479 */ 360 if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) return mindist; 480 if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) 481 { 482 return ; 483 } 361 484 362 363 485 /* 364 486 * Its in the outer ring. 365 487 * Have to check if its inside a hole … … 372 494 * Its inside a hole, then the actual 373 495 * distance is the min ring distance 374 496 */ 375 return mindist;497 return; 376 498 } 377 499 } 500 if (dl->thedir >0) 501 { 502 dl->d=0.0; 503 dl->p1.x=pt.x; 504 dl->p1.y=pt.y; 505 dl->p2.x=pt.x; 506 dl->p2.y=pt.y; 507 } 508 return ; /* Not in hole, so inside polygon */ 378 509 379 return 0.0; /* Not in hole, so inside polygon */380 510 } 381 511 382 double distance2d_point_point(LWPOINT *point1, LWPOINT *point2) 512 513 514 515 516 void lw_dist2d_comp_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl) 383 517 { 384 518 POINT2D p1; 385 519 POINT2D p2; … … 387 521 getPoint2d_p(point1->point, 0, &p1); 388 522 getPoint2d_p(point2->point, 0, &p2); 389 523 390 return distance2d_pt_pt(&p1, &p2); 524 lw_dist2d_comp_pt_pt(&p1, &p2,dl); 525 return; 391 526 } 392 527 393 double distance2d_point_line(LWPOINT *point, LWLINE *line)528 void lw_dist2d_comp_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl) 394 529 { 395 530 POINT2D p; 396 531 POINTARRAY *pa = line->points; 397 532 getPoint2d_p(point->point, 0, &p); 398 return distance2d_pt_ptarray(&p, pa); 533 lw_dist2d_comp_pt_ptarray(&p, pa, dl); 534 return; 399 535 } 400 536 401 double distance2d_line_line(LWLINE *line1, LWLINE *line2)537 void lw_dist2d_comp_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl) 402 538 { 403 539 POINTARRAY *pa1 = line1->points; 404 540 POINTARRAY *pa2 = line2->points; 405 return distance2d_ptarray_ptarray(pa1, pa2); 541 lw_dist2d_comp_ptarray_ptarray(pa1, pa2, dl); 542 return; 406 543 } 407 544 408 545 /* … … 410 547 * 2. if in the boundary, test to see if its in a hole. 411 548 * if so, then return dist to hole, else return 0 (point in polygon) 412 549 */ 413 double distance2d_point_poly(LWPOINT *point, LWPOLY *poly) 550 551 void lw_dist2d_comp_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl) 414 552 { 415 553 POINT2D p; 416 554 int i; 417 555 418 556 getPoint2d_p(point->point, 0, &p); 419 557 420 LWDEBUG(2, " distance2d_point_poly called");558 LWDEBUG(2, "lw_dist2d_comp_point_poly called"); 421 559 560 561 if (dl->thedir < 0) 562 { 563 LWDEBUG(3, "looking for maxdistance"); 564 lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl); 565 return; 566 } 422 567 /* Return distance to outer ring if not inside it */ 423 568 if ( ! pt_in_ring_2d(&p, poly->rings[0]) ) 424 569 { 425 570 LWDEBUG(3, " not inside outer-ring"); 426 427 return distance2d_pt_ptarray(&p, poly->rings[0]);571 lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl); 572 return; 428 573 } 429 574 430 575 /* … … 439 584 if ( pt_in_ring_2d(&p, poly->rings[i]) ) 440 585 { 441 586 LWDEBUG(3, " inside an hole"); 442 443 return distance2d_pt_ptarray(&p, poly->rings[i]);587 lw_dist2d_comp_pt_ptarray(&p, poly->rings[i], dl); 588 return; 444 589 } 445 590 } 446 591 447 592 LWDEBUG(3, " inside the polygon"); 448 449 return 0.0; /* Is inside the polygon */ 593 if (dl->thedir >0) 594 { 595 dl->d=0.0; 596 dl->p1.x=p.x; 597 dl->p1.y=p.y; 598 dl->p2.x=p.x; 599 dl->p2.y=p.y; 600 } 601 return ; /* Is inside the polygon */ 450 602 } 451 452 603 /* 453 * Brute force.454 * Test to see if any rings intersect. 455 * If yes, dist=0. 456 * Test to see if one inside the other and if they are inside holes. 457 * Find min distance ring-to-ring.604 1 if we are looking for maxdistance, just check the outer rings. 605 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings 606 3 check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole of poly1 607 4 check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole of poly2 608 5 If we have come all the way here we know that the first point of one of them is inside the other ones outer ring and not in holes so we check wich one is inside. 458 609 */ 459 double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2)610 void lw_dist2d_comp_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl) 460 611 { 612 461 613 POINT2D pt; 462 double mindist = -1;463 614 int i; 615 LWDEBUG(2, "lw_dist2d_comp_poly_poly called"); 464 616 465 LWDEBUG(2, "distance2d_poly_poly called"); 617 /*1 if we are looking for maxdistance, just check the outer rings.*/ 618 if (dl->thedir <0) 619 {lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); 620 return; 621 } 466 622 467 /* if poly1 inside poly2 return 0 */468 getPoint2d_p(poly1->rings[0], 0, &pt);469 if ( pt_in_poly_2d(&pt, poly2) ) return 0.0;470 623 471 /* if poly2 inside poly1 return 0 */ 472 getPoint2d_p(poly2->rings[0], 0, &pt); 473 if ( pt_in_poly_2d(&pt, poly1) ) return 0.0; 474 475 LWDEBUG(3, " polys not inside each other"); 476 477 /* 478 * foreach ring in Poly1 479 * foreach ring in Poly2 480 * if intersect, return 0 481 */ 482 for (i=0; i<poly1->nrings; i++) 624 /* 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings 625 here it would be possible to handle the information about wich on is inside wichone and only search for the smaller ones in the bigger ones holes.*/ 626 getPoint2d_p(poly1->rings[0], 0, &pt); 627 if ( !pt_in_ring_2d(&pt, poly2->rings[0])) 483 628 { 484 int j; 485 for (j=0; j<poly2->nrings; j++) 629 getPoint2d_p(poly2->rings[0], 0, &pt); 630 if (!pt_in_ring_2d(&pt, poly1->rings[0])) 631 { 632 lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); 633 return; 634 } 635 } 636 637 /*3 check if first point of poly2 is in a hole of poly1. If so check outer ring of poly2 against that hole of poly1*/ 638 getPoint2d_p(poly2->rings[0], 0, &pt); 639 for (i=1; i<poly1->nrings; i++) 640 { 641 /* Inside a hole */ 642 if ( pt_in_ring_2d(&pt, poly1->rings[i]) ) 486 643 { 487 double d = distance2d_ptarray_ptarray(poly1->rings[i], 488 poly2->rings[j]); 489 if ( d <= 0 ) return 0.0; 644 lw_dist2d_comp_ptarray_ptarray(poly1->rings[i], poly2->rings[0],dl); 645 return; 646 } 647 648 } 649 650 /*4 check if first point of poly1 is in a hole of poly2. If so check outer ring of poly1 against that hole of poly2*/ 651 getPoint2d_p(poly1->rings[0], 0, &pt); 652 for (i=1; i<poly2->nrings; i++) 653 { 654 /* Inside a hole */ 655 if ( pt_in_ring_2d(&pt, poly2->rings[i]) ) 656 { 657 lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[i],dl); 658 return; 659 } 660 661 } 662 490 663 491 /* mindist is -1 when not yet set */ 492 if (mindist > -1) mindist = LW_MIN(mindist, d); 493 else mindist = d; 664 /*5 If we have come all the way here we know that the first point of one of them is inside the other ones outer ring and not in holes so we check wich one is inside.*/ 665 getPoint2d_p(poly1->rings[0], 0, &pt); 666 if ( pt_in_ring_2d(&pt, poly2->rings[0])) 667 { 668 dl->d=0.0; 669 dl->p1.x=pt.x; 670 dl->p1.y=pt.y; 671 dl->p2.x=pt.x; 672 dl->p2.y=pt.y; 673 return ; 674 } 494 675 495 LWDEBUGF(3, " ring%i-%i dist: %f, mindist: %f", i, j, d, mindist); 496 } 497 676 getPoint2d_p(poly2->rings[0], 0, &pt); 677 if (pt_in_ring_2d(&pt, poly1->rings[0])) 678 { 679 dl->d=0.0; 680 dl->p1.x=pt.x; 681 dl->p1.y=pt.y; 682 dl->p2.x=pt.x; 683 dl->p2.y=pt.y; 684 return ; 498 685 } 499 500 /* otherwise return closest approach of rings (no intersection) */501 return mindist;502 686 687 688 lwerror("You have found a bug. You could report that lw_dist2d_comp_poly_poly is not working as it should"); 689 503 690 } 504 691 505 double distance2d_line_poly(LWLINE *line, LWPOLY *poly)692 void lw_dist2d_comp_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl) 506 693 { 507 return distance2d_ptarray_poly(line->points, poly); 694 lw_dist2d_comp_ptarray_poly(line->points, poly, dl); 695 return; 508 696 } 509 697 510 511 698 /*find the 2d length of the given POINTARRAY (even if it's 3d) */ 512 699 double lwgeom_pointarray_length2d(POINTARRAY *pts) 513 700 { … … 552 739 ((frm.y - to.y)*(frm.y - to.y) ) + 553 740 ((frm.z - to.z)*(frm.z - to.z) ) ); 554 741 } 555 556 742 return dist; 557 743 } 558 744 … … 645 831 } 646 832 647 833 double 834 lwgeom_maxdistance2d_recursive(uchar *lw1, uchar *lw2) 835 { 836 return lwgeom_maxdistance2d_recursive_tolerance( lw1, lw2, 0.0 ); 837 } 838 839 840 double 841 lwgeom_maxdistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) 842 { 843 /*double thedist;*/ 844 DISTPTS thedl; 845 thedl.thedir = (-1); 846 thedl.d= -1.0; 847 thedl.tolerance = tolerance; 848 DISTPTS *dl = &thedl; 849 lw_dist2d_comp(lw1,lw2,dl); 850 return dl->d; 851 /*lwgeom_release((DISTPTS *)dl); 852 return thedist; 853 return 1212.00;*/ 854 } 855 856 double 648 857 lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) 649 858 { 650 859 return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 ); 651 860 } 652 861 862 653 863 double 654 864 lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) 655 865 { 866 /*double thedist;*/ 867 DISTPTS thedl; 868 thedl.thedir = (1); 869 thedl.d= 99999999.0; 870 thedl.tolerance = tolerance; 871 DISTPTS *dl = &thedl; 872 lw_dist2d_comp( lw1,lw2,dl); 873 return dl->d; 874 /*lwgeom_release((DISTPTS *)dl); 875 return thedist; 876 return 1212.00;*/ 877 } 878 879 void 880 lw_dist2d_comp(uchar *lw1, uchar *lw2, DISTPTS *dl) 881 { 656 882 LWGEOM_INSPECTED *in1, *in2; 657 883 int i, j; 658 double mindist = -1;659 884 660 885 in1 = lwgeom_inspect(lw1); 661 886 in2 = lwgeom_inspect(lw2); … … 664 889 { 665 890 uchar *g1 = lwgeom_getsubgeometry_inspected(in1, i); 666 891 int t1 = lwgeom_getType(g1[0]); 667 double dist=tolerance;668 892 669 /* Argument 1 is a multitype... recurse */893 /* it's a multitype... recurse */ 670 894 if ( lwgeom_contains_subgeoms(t1) ) 671 895 { 672 dist = lwgeom_mindistance2d_recursive_tolerance(g1, lw2, tolerance); 673 if ( dist <= tolerance ) return tolerance; /* can't be closer */ 674 if ( mindist == -1 ) mindist = dist; 675 else mindist = LW_MIN(dist, mindist); 896 lw_dist2d_comp(g1, lw2, dl); 897 898 if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ 676 899 continue; 677 900 } 678 901 … … 681 904 uchar *g2 = lwgeom_getsubgeometry_inspected(in2, j); 682 905 int t2 = lwgeom_getType(g2[0]); 683 906 684 /* Argument 2 is a multitype... recurse */685 907 if ( lwgeom_contains_subgeoms(t2) ) 686 908 { 687 dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); 688 if ( dist <= tolerance ) return tolerance; /* can't be closer */ 689 if ( mindist == -1 ) mindist = dist; 690 else mindist = LW_MIN(dist, mindist); 909 lw_dist2d_comp(g1, g2, dl); 910 if (dl->d<=dl->tolerance) 911 { 912 /*dl->d=tolerance;*/ 913 return; /* can't be closer */ 914 } 691 915 continue; 692 916 } 693 694 917 if ( t1 == POINTTYPE ) 695 918 { 696 919 if ( t2 == POINTTYPE ) 697 920 { 698 dist = distance2d_point_point( 699 lwpoint_deserialize(g1), 700 lwpoint_deserialize(g2) 701 ); 921 dl->twisted=1; 922 lw_dist2d_comp_point_point( 923 lwpoint_deserialize(g1), 924 lwpoint_deserialize(g2), 925 dl); 702 926 } 703 927 else if ( t2 == LINETYPE ) 704 928 { 705 dist = distance2d_point_line( 706 lwpoint_deserialize(g1), 707 lwline_deserialize(g2) 708 ); 929 dl->twisted=1; 930 lw_dist2d_comp_point_line( 931 lwpoint_deserialize(g1), 932 lwline_deserialize(g2), 933 dl); 709 934 } 710 935 else if ( t2 == POLYGONTYPE ) 711 936 { 712 dist = distance2d_point_poly( 713 lwpoint_deserialize(g1), 714 lwpoly_deserialize(g2) 715 ); 937 dl->twisted=1; 938 lw_dist2d_comp_point_poly( 939 lwpoint_deserialize(g1), 940 lwpoly_deserialize(g2), 941 dl); 716 942 } 717 943 else 718 944 { … … 723 949 { 724 950 if ( t2 == POINTTYPE ) 725 951 { 726 dist = distance2d_point_line( 727 lwpoint_deserialize(g2), 728 lwline_deserialize(g1) 729 ); 952 dl->twisted=(-1); 953 lw_dist2d_comp_point_line( 954 lwpoint_deserialize(g2), 955 lwline_deserialize(g1), 956 dl); 730 957 } 731 958 else if ( t2 == LINETYPE ) 732 959 { 733 dist = distance2d_line_line( 734 lwline_deserialize(g1), 735 lwline_deserialize(g2) 736 ); 960 dl->twisted=1; 961 /*lwnotice("start line");*/ 962 lw_dist2d_comp_line_line( 963 lwline_deserialize(g1), 964 lwline_deserialize(g2), 965 dl); 737 966 } 738 967 else if ( t2 == POLYGONTYPE ) 739 968 { 740 dist = distance2d_line_poly( 741 lwline_deserialize(g1), 742 lwpoly_deserialize(g2) 743 ); 969 dl->twisted=1; 970 lw_dist2d_comp_line_poly( 971 lwline_deserialize(g1), 972 lwpoly_deserialize(g2), 973 dl); 744 974 } 745 975 else 746 976 { … … 751 981 { 752 982 if ( t2 == POLYGONTYPE ) 753 983 { 754 dist = distance2d_poly_poly( 755 lwpoly_deserialize(g2), 756 lwpoly_deserialize(g1) 757 ); 984 dl->twisted=(-1); 985 lw_dist2d_comp_poly_poly( 986 lwpoly_deserialize(g2), 987 lwpoly_deserialize(g1), 988 dl); 758 989 } 759 990 else if ( t2 == POINTTYPE ) 760 991 { 761 dist = distance2d_point_poly( 762 lwpoint_deserialize(g2), 763 lwpoly_deserialize(g1) 764 ); 992 dl->twisted=(-1); 993 lw_dist2d_comp_point_poly( 994 lwpoint_deserialize(g2), 995 lwpoly_deserialize(g1), 996 dl); 765 997 } 766 998 else if ( t2 == LINETYPE ) 767 999 { 768 dist = distance2d_line_poly( 769 lwline_deserialize(g2), 770 lwpoly_deserialize(g1) 771 ); 1000 dl->twisted=(-1); 1001 lw_dist2d_comp_line_poly( 1002 lwline_deserialize(g2), 1003 lwpoly_deserialize(g1), 1004 dl); 772 1005 } 773 1006 else 774 1007 { 775 1008 lwerror("Unsupported geometry type: %s", lwgeom_typename(t2)); 776 1009 } 777 1010 } 778 // else if (lwgeom_contains_subgeoms(t1)) /* it's a multitype... recurse */ 779 // { 780 // dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); 781 // } 1011 /* else if (lwgeom_contains_subgeoms(t1)) 1012 { 1013 lw_dist2d_comp(g1, g2, tolerance,dl); 1014 1015 }*/ 782 1016 else 783 1017 { 784 1018 lwerror("Unsupported geometry type: %s", lwgeom_typename(t1)); 785 1019 } 786 1020 787 if (mindist == -1 ) mindist = dist;788 else mindist = LW_MIN(dist, mindist);789 1021 790 LWDEBUGF(3, "dist %d-%d: %f - mindist: %f",791 i, j, dist, mindist);792 1022 793 1023 794 if (mindist <= tolerance) return tolerance; /* can't be closer*/1024 if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ 795 1025 796 1026 } 797 1027 798 1028 } 799 800 if (mindist<0) mindist = 0; 801 802 return mindist; 1029 return ; 803 1030 } 804 1031 805 1032 -
postgis/lwgeom_functions_basic.c
33 33 Datum LWGEOM_nrings(PG_FUNCTION_ARGS); 34 34 Datum LWGEOM_area_polygon(PG_FUNCTION_ARGS); 35 35 Datum LWGEOM_dwithin(PG_FUNCTION_ARGS); 36 Datum LWGEOM_dfullywithin(PG_FUNCTION_ARGS); 36 37 Datum postgis_uses_stats(PG_FUNCTION_ARGS); 37 38 Datum postgis_autocache_bbox(PG_FUNCTION_ARGS); 38 39 Datum postgis_scripts_released(PG_FUNCTION_ARGS); … … 43 44 Datum LWGEOM_length_linestring(PG_FUNCTION_ARGS); 44 45 Datum LWGEOM_perimeter2d_poly(PG_FUNCTION_ARGS); 45 46 Datum LWGEOM_perimeter_poly(PG_FUNCTION_ARGS); 47 Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS); 46 48 Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS); 47 Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS); 49 Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS); 50 Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS); 48 51 Datum LWGEOM_inside_circle_point(PG_FUNCTION_ARGS); 49 52 Datum LWGEOM_collect(PG_FUNCTION_ARGS); 50 53 Datum LWGEOM_accum(PG_FUNCTION_ARGS); … … 1522 1525 PG_RETURN_POINTER(result); 1523 1526 } 1524 1527 1528 1529 PG_FUNCTION_INFO_V1(LWGEOM_shortestline2d); 1530 Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS) 1531 { 1532 1533 PG_LWGEOM *geom1; 1534 PG_LWGEOM *geom2; 1535 LWPOINT *point1; 1536 LWPOINT *point2; 1537 int SRID; 1538 double x1,x2,y1,y2; 1539 PG_LWGEOM *result=NULL; 1540 LWPOINT *lwpoints[2]; 1541 LWLINE *outline; 1542 1543 DISTPTS thedl; 1544 DISTPTS *dl = &thedl; 1545 dl->d=99999999.0; 1546 dl->thedir=(1); 1547 dl->tolerance = 0.0; 1548 1549 geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1550 geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1551 1552 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1553 { 1554 elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); 1555 PG_RETURN_NULL(); 1556 } 1557 1558 SRID = pglwgeom_getSRID(geom1); 1559 1560 lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2), dl ); 1561 1562 x1=dl->p1.x; 1563 y1=dl->p1.y; 1564 x2=dl->p2.x; 1565 y2=dl->p2.y; 1566 1567 point1 = make_lwpoint2d(SRID, x1, y1); 1568 point2 = make_lwpoint2d(SRID, x2, y2); 1569 1570 lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1))); 1571 lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2))); 1572 1573 outline = lwline_from_lwpointarray(SRID, 2, lwpoints); 1574 1575 result = pglwgeom_serialize((LWGEOM *)outline); 1576 PG_RETURN_POINTER(result); 1577 1578 } 1579 1580 PG_FUNCTION_INFO_V1(LWGEOM_longestline2d); 1581 Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS) 1582 { 1583 1584 PG_LWGEOM *geom1; 1585 PG_LWGEOM *geom2; 1586 LWPOINT *point1; 1587 LWPOINT *point2; 1588 int SRID; 1589 1590 double x1,x2,y1,y2; 1591 PG_LWGEOM *result=NULL; 1592 LWPOINT *lwpoints[2]; 1593 LWLINE *outline; 1594 DISTPTS thedl; 1595 DISTPTS *dl = &thedl; 1596 dl->tolerance = 0.0; 1597 dl->d=(-1.0); 1598 dl->thedir=(-1); 1599 1600 geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1601 geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1602 1603 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1604 { 1605 elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); 1606 PG_RETURN_NULL(); 1607 } 1608 1609 SRID = pglwgeom_getSRID(geom1); 1610 1611 lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2), dl ); 1612 1613 x1=dl->p1.x; 1614 y1=dl->p1.y; 1615 x2=dl->p2.x; 1616 y2=dl->p2.y; 1617 1618 point1 = make_lwpoint2d(SRID, x1, y1); 1619 point2 = make_lwpoint2d(SRID, x2, y2); 1620 1621 lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1))); 1622 lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2))); 1623 1624 outline = lwline_from_lwpointarray(SRID, 2, lwpoints); 1625 1626 result = pglwgeom_serialize((LWGEOM *)outline); 1627 1628 PG_RETURN_POINTER(result); 1629 } 1630 1525 1631 /* Minimum 2d distance between objects in geom1 and geom2. */ 1526 1632 PG_FUNCTION_INFO_V1(LWGEOM_mindistance2d); 1527 1633 Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS) … … 1594 1700 PG_RETURN_BOOL(tolerance >= mindist); 1595 1701 } 1596 1702 1597 /* 1598 * Maximum 2d distance between linestrings. 1599 * Returns NULL if given geoms are not linestrings. 1600 * This is very bogus (or I'm missing its meaning) 1601 */ 1703 PG_FUNCTION_INFO_V1(LWGEOM_dfullywithin); 1704 Datum LWGEOM_dfullywithin(PG_FUNCTION_ARGS) 1705 { 1706 PG_LWGEOM *geom1; 1707 PG_LWGEOM *geom2; 1708 double maxdist, tolerance; 1709 1710 PROFSTART(PROF_QRUN); 1711 1712 geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1713 geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1714 tolerance = PG_GETARG_FLOAT8(2); 1715 1716 if ( tolerance < 0 ) 1717 { 1718 elog(ERROR,"Tolerance cannot be less than zero\n"); 1719 PG_RETURN_NULL(); 1720 } 1721 1722 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1723 { 1724 elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); 1725 PG_RETURN_NULL(); 1726 } 1727 maxdist = lwgeom_maxdistance2d_recursive_tolerance( 1728 SERIALIZED_FORM(geom1), 1729 SERIALIZED_FORM(geom2), 1730 tolerance 1731 ); 1732 PROFSTOP(PROF_QRUN); 1733 PROFREPORT("dist",geom1, geom2, NULL); 1734 1735 PG_FREE_IF_COPY(geom1, 0); 1736 PG_FREE_IF_COPY(geom2, 1); 1737 1738 PG_RETURN_BOOL(tolerance >= maxdist); 1739 } 1740 1602 1741 PG_FUNCTION_INFO_V1(LWGEOM_maxdistance2d_linestring); 1603 1742 Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS) 1604 1743 { 1605 1606 1744 PG_LWGEOM *geom1; 1607 1745 PG_LWGEOM *geom2; 1608 LWLINE *line1; 1609 LWLINE *line2; 1610 double maxdist = 0; 1611 int i; 1746 double maxdist; 1612 1747 1613 elog(ERROR, "This function is unimplemented yet"); 1614 PG_RETURN_NULL(); 1748 PROFSTART(PROF_QRUN); 1615 1749 1616 1750 geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1617 line1 = lwline_deserialize(SERIALIZED_FORM(geom1));1618 if ( line1 == NULL ) PG_RETURN_NULL(); /* not a linestring */1619 1620 1751 geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1621 line2 = lwline_deserialize(SERIALIZED_FORM(geom2));1622 if ( line2 == NULL ) PG_RETURN_NULL(); /* not a linestring */1623 1752 1624 1753 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1625 1754 { … … 1627 1756 PG_RETURN_NULL(); 1628 1757 } 1629 1758 1630 for (i=0; i<line1->points->npoints; i++)1631 {1632 POINT2D p;1633 double dist;1634 1759 1635 getPoint2d_p(line1->points, i, &p); 1636 dist = distance2d_pt_ptarray(&p, line2->points); 1760 maxdist = lwgeom_maxdistance2d_recursive(SERIALIZED_FORM(geom1), 1761 SERIALIZED_FORM(geom2)); 1762 /* 1763 mindist = lwgeom_mindistance2d_recursive_tolerance( 1764 SERIALIZED_FORM(geom1), 1765 SERIALIZED_FORM(geom2), 1766 tolerance 1767 ); 1768 */ 1637 1769 1638 if (dist > maxdist) maxdist = dist;1639 }1770 PROFSTOP(PROF_QRUN); 1771 PROFREPORT("maxdist",geom1, geom2, NULL); 1640 1772 1641 1773 PG_FREE_IF_COPY(geom1, 0); 1642 1774 PG_FREE_IF_COPY(geom2, 1); … … 1644 1776 PG_RETURN_FLOAT8(maxdist); 1645 1777 } 1646 1778 1647 /** 1648 * @brief Longitude shift: 1649 * Y remains the same 1650 * X is converted: 1651 * from -180..180 to 0..360 1652 * from 0..360 to -180..180 1653 */ 1779 1654 1780 PG_FUNCTION_INFO_V1(LWGEOM_longitude_shift); 1655 1781 Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS) 1656 1782 { -
regress/measures.sql
47 47 ) as foo; 48 48 49 49 50 51 52 --st_max_distance 53 54 select 'st_max_distance_134', st_max_distance('POINT(1 2)', 'POINT(1 2)'); 55 select 'st_max_distance_135', st_max_distance('POINT(5 0)', 'POINT(10 12)'); 56 57 select 'st_max_distance_136', st_max_distance('POINT(0 0)', ST_translate('POINT(0 0)', 5, 12, 0)); 58 59 -- postgis-users/2006-May/012174.html 60 select 'st_max_distance_dist', st_max_distance(a,b), st_max_distance(b,a) from ( 61 select 'POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))'::geometry as a, 62 'POLYGON((11 0, 11 10, 20 10, 20 0, 11 0), 63 (15 5, 15 8, 17 8, 17 5, 15 5))'::geometry as b 64 ) as foo; 65 66 67 68 select 'distancetest1', 69 st_distance(a, b), 70 st_max_distance(a, b) 71 from ( 72 select geomfromtext('MULTILINESTRING((17 16, 16 17, 17 18, 18 17, 17 16), (28 35,29 39, 30 35))') as a, 73 geomfromtext('MULTIPOLYGON(((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14)), 74 ((33 35,33 40, 35 40, 35 35, 33 35)))') as b 75 ) as foo; 76 77 78 79 select 'distancetest2', 80 st_distance(a, b), 81 st_max_distance(a, b) 82 from ( 83 select geomfromtext('LINESTRING(-40 -20 , 4 2)') as a, 84 geomfromtext('LINESTRING(-10 20, 1 -2)') as b 85 ) as foo; 86 87 88 89 select 'distancepoly1', 90 st_distance(a, b), 91 st_max_distance(a, b) 92 from ( 93 select geomfromtext('MULTIPOLYGON(((17 16, 16 17, 17 18, 18 17, 17 16)), ((28 35,29 39, 30 35, 28 35)))') as a, 94 geomfromtext('MULTIPOLYGON(((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14)), 95 ((33 35,33 40, 35 40, 35 35, 33 35)))') as b 96 ) as foo; 97 98 99 100 select 'distancepoly2', 101 st_distance(a, b), 102 st_max_distance(a, b) 103 from ( 104 select geomfromtext('POLYGON((17 14, 16 17, 17 18, 18 17, 17 14))') as a, 105 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 106 ) as foo; 107 108 109 110 select 'distancepoly3', 111 st_distance(a, b), 112 st_max_distance(a, b) 113 from ( 114 select geomfromtext('POLYGON((17 16, 16 17, 17 19, 18 17, 17 16))') as a, 115 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 116 ) as foo; 117 118 119 select 'distancepoly4', 120 st_distance(a, b), 121 st_max_distance(a, b) 122 from ( 123 select geomfromtext('POLYGON((17 16, 16 17, 16 20, 18 20, 18 17, 17 16))') as a, 124 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 125 ) as foo; 126 127 128 129 select 'distancepoly5', 130 st_distance(a, b), 131 st_max_distance(a, b) 132 from ( 133 select geomfromtext('POLYGON((17 12, 16 17, 17 18, 18 17, 17 12))') as a, 134 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 135 ) as foo; 136 137 138 139 140 select 'distancepoly6', 141 st_distance(a, b), 142 st_max_distance(a, b) 143 from ( 144 select geomfromtext('POLYGON((2 2, 2 3, 3 3, 3 2, 2 2))') as a, 145 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 146 ) as foo; 147 148 149 150 50 151 -- Area of an empty polygon 51 152 select 'emptyPolyArea', st_area('POLYGON EMPTY'); 52 153 -
regress/measures_expected
18 18 135|13 19 19 136|13 20 20 dist|1|1 21 st_max_distance_134|0 22 st_max_distance_135|13 23 st_max_distance_136|13 24 st_max_distance_dist|22.3606797749979|22.3606797749979 25 distancetest1|1|50 26 distancetest2|0|50 27 distancepoly1|1|30 28 distancepoly2|0|26.1725046566048 29 distancepoly3|0|26.9072480941474 30 distancepoly4|0|28.3196045170126 31 distancepoly5|0|26.1725046566048 32 distancepoly6|0|32.5269119345812 21 33 emptyPolyArea|0 22 34 emptyLineArea|0 23 35 emptyPointArea|0
