Index: doc/reference.xml =================================================================== --- doc/reference.xml (revision 4310) +++ doc/reference.xml (working copy) @@ -9833,6 +9833,54 @@ + + + ST_Max_Distance + + Returns the 2-dimensional largest distance between two geometries in + projected units. + + + + + + float ST_Max_Distance + + geometry + g1 + + geometry + g2 + + + + + + Description + + Returns the 2-dimensional maximum cartesian distance between two linestrings in + projected units. If g1 and g2 is the same geometry the function will return the distance between + the two vertices most far from eachother in that geometry. + + + + + Examples + + postgis=# SELECT ST_Max_Distance('POINT(0 0)'::geometry, 'LINESTRING ( 2 0, 0 2 )'::geometry); + st_max_distance +----------------- + 2 +(1 row) + + + + See Also + + + + + ST_Distance_Sphere Index: liblwgeom/liblwgeom.h =================================================================== --- liblwgeom/liblwgeom.h (revision 4310) +++ liblwgeom/liblwgeom.h (working copy) @@ -1,3 +1,4 @@ + /********************************************************************** * $Id$ * @@ -19,6 +20,7 @@ #include #include + /** * @file liblwgeom.h * @@ -290,6 +292,15 @@ } POINT4D; +typedef struct +{ + double d; /*the distance between p1 and p2*/ + POINT2D p1; + POINT2D p2; + int thedir; /*the direction of looking, if thedir = -1 then we look for maxdistance and if it is 1 then we look for mindistance*/ + int twisted; /*To preserve the order of incoming points to match the first and secon point in shortest and longest line*/ + double tolerance; /*the tolerance for dwithin and dfullywithin*/ +} DISTPTS; /******************************************************************/ /* @@ -1184,23 +1195,29 @@ extern void lwgeom_force3dz_recursive(uchar *serialized, uchar *optr, size_t *retsize); extern void lwgeom_force3dm_recursive(uchar *serialized, uchar *optr, size_t *retsize); extern void lwgeom_force4d_recursive(uchar *serialized, uchar *optr, size_t *retsize); + extern double distance2d_pt_pt(POINT2D *p1, POINT2D *p2); +extern void lw_dist2d_comp_pt_pt(POINT2D *p1, POINT2D *p2, DISTPTS *dl); extern double distance2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B); -extern double distance2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D); -extern double distance2d_pt_ptarray(POINT2D *p, POINTARRAY *pa); -extern double distance2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2); +extern void lw_dist2d_comp_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl); +extern void lw_dist2d_comp_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl); +extern void lw_dist2d_comp_pt_ptarray(POINT2D *p, POINTARRAY *pa, DISTPTS *dl); +extern void lw_dist2d_comp_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl); extern int pt_in_ring_2d(POINT2D *p, POINTARRAY *ring); extern int pt_in_poly_2d(POINT2D *p, LWPOLY *poly); -extern double distance2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly); -extern double distance2d_point_point(LWPOINT *point1, LWPOINT *point2); -extern double distance2d_point_line(LWPOINT *point, LWLINE *line); -extern double distance2d_line_line(LWLINE *line1, LWLINE *line2); -extern double distance2d_point_poly(LWPOINT *point, LWPOLY *poly); -extern double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2); -extern double distance2d_line_poly(LWLINE *line, LWPOLY *poly); +extern void lw_dist2d_comp_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl); +extern void lw_dist2d_comp_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl); +extern void lw_dist2d_comp_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl); +extern void lw_dist2d_comp_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl); +extern void lw_dist2d_comp_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl); +extern void lw_dist2d_comp_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl); +extern void lw_dist2d_comp_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl); extern int azimuth_pt_pt(POINT2D *p1, POINT2D *p2, double *ret); extern double lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2); extern double lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance); +extern double lwgeom_maxdistance2d_recursive(uchar *lw1, uchar *lw2); +extern double lwgeom_maxdistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance); +extern void lw_dist2d_comp(uchar *lw1, uchar *lw2, DISTPTS *dl); extern int lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad); extern int32 lwgeom_npoints(uchar *serialized); extern char ptarray_isccw(const POINTARRAY *pa); Index: liblwgeom/measures.c =================================================================== --- liblwgeom/measures.c (revision 4310) +++ liblwgeom/measures.c (working copy) @@ -81,6 +81,31 @@ return (cn&1); /* 0 if even (out), and 1 if odd (in) */ } +/*Compares incomming points and stores the points closest to each other or most far away from each other depending on dl->thedir (thedirection) */ +void lw_dist2d_comp_pt_pt(POINT2D *thep1, POINT2D *thep2,DISTPTS *dl) +{ + double hside = thep2->x - thep1->x; + double vside = thep2->y - thep1->y; + double dist = sqrt ( hside*hside + vside*vside ); + if (((dl->d - dist)*(dl->thedir))>0) /*multiplication with thedir to handle mindistance (thedir=1) and maxdistance (thedir = (-1)*/ + { + dl->d = dist; + + 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*/ + { + dl->p1 = *thep1; + dl->p2 = *thep2; + } + else + { + dl->p1 = *thep2; + dl->p2 = *thep1; + } + } + return; +} + +/*The old function nessecary for ptarray_segmentize2d in ptarray.c*/ double distance2d_pt_pt(POINT2D *p1, POINT2D *p2) { double hside = p2->x - p1->x; @@ -94,7 +119,81 @@ ); */ } -/*distance2d from p to line A->B */ +/*lw_dist2d_comp from p to line A->B +This one is now sending every occation to lw_dist2d_comp_pt_pt +Berfore it was handling occations where r was between 0 and 1 internally and just returning the distance without identifying the points. +To get this points it was nessecary to change and it also showed to be about 10%faster. */ +void lw_dist2d_comp_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl) +{ + /*lwnotice("lw_dist2d_comp_pt_seg %e",dl->d);*/ + double r; + + /*if start==end, then use pt distance */ + if ( ( A->x == B->x) && (A->y == B->y) ) + { + lw_dist2d_comp_pt_pt(p,A,dl); + return; + } + /* + * otherwise, we use comp.graphics.algorithms + * Frequently Asked Questions method + * + * (1) AC dot AB + * r = --------- + * ||AB||^2 + * r has the following meaning: + * r=0 P = A + * r=1 P = B + * r<0 P is on the backward extension of AB + * r>1 P is on the forward extension of AB + * 0x-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) ); + + /*This is for finding the maxdistance. + the maxdistance have to be between two vertexes as I understand it, + compared to mindistance which can be between some point between the vertices and a vertex.*/ + if (dl->thedir < 1) + { + if (r>=0.5) + { + lw_dist2d_comp_pt_pt(p,A,dl); + return; + } + if (r<0.5) + { + lw_dist2d_comp_pt_pt(p,B,dl); + return; + } + } + + if (r<0) /*If the first vertex A is closest to the point p*/ + { + lw_dist2d_comp_pt_pt(p,A,dl); + return; + } + if (r>1) /*If the second vertex B is closest to the point p*/ + { + + lw_dist2d_comp_pt_pt(p,B,dl); + return; + } + + + /*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*/ + POINT2D c; + c.x=A->x + r * (B->x-A->x); + c.y=A->y + r * (B->y-A->y); + + lw_dist2d_comp_pt_pt(p,&c,dl); + return; +} + + + + +/*The old function nessecary for ptarray_segmentize2d in ptarray.c*/ double distance2d_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B) { double r,s; @@ -142,26 +241,36 @@ ); } -/* find the minimum 2d distance from AB to CD */ -double distance2d_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D) + + + + +/*This function is changed so it is not doing any comparasion of distance +but just sending every possible combination further to lw_dist2d_comp_pt_seg*/ +void lw_dist2d_comp_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl) { double s_top, s_bot,s; double r_top, r_bot,r; - LWDEBUGF(2, "distance2d_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]", + + LWDEBUGF(2, "lw_dist2d_comp_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]", A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y); - /*A and B are the same point */ if ( ( A->x == B->x) && (A->y == B->y) ) - return distance2d_pt_seg(A,C,D); - + { + lw_dist2d_comp_pt_seg(A,C,D,dl); + return; + } /*U and V are the same point */ if ( ( C->x == D->x) && (C->y == D->y) ) - return distance2d_pt_seg(D,A,B); - + { + dl->twisted= ((dl->twisted) * (-1)); + lw_dist2d_comp_pt_seg(D,A,B,dl); + return; + } /* AB and CD are line segments */ /* from comp.graphics.algo @@ -192,113 +301,123 @@ if ( (r_bot==0) || (s_bot == 0) ) { - return ( - LW_MIN(distance2d_pt_seg(A,C,D), - LW_MIN(distance2d_pt_seg(B,C,D), - LW_MIN(distance2d_pt_seg(C,A,B), - distance2d_pt_seg(D,A,B)) - ) - ) - ); + lw_dist2d_comp_pt_seg(A,C,D,dl); + lw_dist2d_comp_pt_seg(B,C,D,dl); + + dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/ + lw_dist2d_comp_pt_seg(C,A,B,dl); + lw_dist2d_comp_pt_seg(D,A,B,dl); + return; } + s = s_top/s_bot; r= r_top/r_bot; - if ((r<0) || (r>1) || (s<0) || (s>1) ) + if (((r<0) || (r>1) || (s<0) || (s>1)) || (dl->thedir <0)) { - /*no intersection */ - return ( - LW_MIN(distance2d_pt_seg(A,C,D), - LW_MIN(distance2d_pt_seg(B,C,D), - LW_MIN(distance2d_pt_seg(C,A,B), - distance2d_pt_seg(D,A,B)) - ) - ) - ); - + lw_dist2d_comp_pt_seg(A,C,D,dl); + lw_dist2d_comp_pt_seg(B,C,D,dl); + + dl->twisted= ((dl->twisted) * (-1)); /*here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted*/ + lw_dist2d_comp_pt_seg(C,A,B,dl); + lw_dist2d_comp_pt_seg(D,A,B,dl); + return; } else - return -0; /*intersection exists */ + { + if (dl->thedir >0) /*If there is intersection we identify the intersection point and return it but only if we are looking for mindistance*/ + { /*lwnotice("intersection lines: %e %e %e %e",A->x,A->y,B->x,B->y);*/ + POINT2D theP; + + if(((A->x==C->x)&&(A->y==C->y))||((A->x==D->x)&&(A->y==D->y))) + { + theP.x = A->x; + theP.y = A->y; + } + else if(((B->x==C->x)&&(B->y==C->y))||((B->x==D->x)&&(B->y==D->y))) + { + theP.x = B->x; + theP.y = B->y; + } + else + { + theP.x = A->x+r*(B->x-A->x); + theP.y = A->y+r*(B->y-A->y); + } + /*lwnotice("intersection: %e %e ",theP.x, theP.y);*/ + dl->d=0.0; + dl->p1=theP; + dl->p2=theP; + } + return; + } + } /* * search all the segments of pointarray to see which one is closest to p1 * Returns minimum distance between point and pointarray */ -double distance2d_pt_ptarray(POINT2D *p, POINTARRAY *pa) +void lw_dist2d_comp_pt_ptarray(POINT2D *p, POINTARRAY *pa,DISTPTS *dl) { - double result = 0; + + int t; POINT2D start, end; - + int twist = dl->twisted; getPoint2d_p(pa, 0, &start); for (t=1; tnpoints; t++) { - double dist; + dl->twisted=twist; getPoint2d_p(pa, t, &end); - dist = distance2d_pt_seg(p, &start, &end); - if (t==1) result = dist; - else result = LW_MIN(result, dist); - - if ( result == 0 ) return 0; - + lw_dist2d_comp_pt_seg(p, &start, &end,dl); + + if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ + start = end; } - return result; + return; } /* test each segment of l1 against each segment of l2. Return min */ -double distance2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2) +void lw_dist2d_comp_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl) { - double result = 99999999999.9; - char result_okay = 0; /*result is a valid min */ + /*lwnotice("lw_dist2d_comp_ptarray_ptarray");*/ int t,u; POINT2D start, end; POINT2D start2, end2; - - LWDEBUGF(2, "distance2d_ptarray_ptarray called (points: %d-%d)", + int twist = dl->twisted; + LWDEBUGF(2, "lw_dist2d_comp_ptarray_ptarray called (points: %d-%d)", l1->npoints, l2->npoints); getPoint2d_p(l1, 0, &start); for (t=1; tnpoints; t++) /*for each segment in L1 */ { getPoint2d_p(l1, t, &end); - getPoint2d_p(l2, 0, &start2); for (u=1; unpoints; u++) /*for each segment in L2 */ { - double dist; - getPoint2d_p(l2, u, &end2); + dl->twisted=twist; + lw_dist2d_comp_seg_seg(&start, &end, &start2, &end2,dl); - dist = distance2d_seg_seg(&start, &end, &start2, &end2); + LWDEBUGF(4, "line_line; seg %i * seg %i, dist = %g\n",t,u,dist); - - if (result_okay) - result = LW_MIN(result,dist); - else - { - result_okay = 1; - result = dist; - } - LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f", t, u, dist, result); - - if (result <= 0) return 0; /*intersection */ - + if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ start2 = end2; } start = end; } - - return result; + return ; } + /* true if point is in poly (and not in its holes) */ int pt_in_poly_2d(POINT2D *p, LWPOLY *poly) { @@ -317,6 +436,8 @@ return 1; /* In outer ring, not in holes */ } + + /* * Brute force. * Test line-ring distance against each ring. @@ -327,24 +448,23 @@ * otherwise return min distance to a ring (could be outside * polygon or inside a hole) */ -double distance2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly) + +void lw_dist2d_comp_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl) { POINT2D pt; int i; - double mindist = 0; - LWDEBUGF(2, "distance2d_ptarray_poly called (%d rings)", poly->nrings); + LWDEBUGF(2, "lw_dist2d_comp_ptarray_poly called (%d rings)", poly->nrings); + + for (i=0; inrings; i++) { - double dist = distance2d_ptarray_ptarray(pa, poly->rings[i]); - if (i) mindist = LW_MIN(mindist, dist); - else mindist = dist; + lw_dist2d_comp_ptarray_ptarray(pa, poly->rings[i], dl); LWDEBUGF(3, " distance from ring %d: %f, mindist: %f", i, dist, mindist); - - if ( mindist <= 0 ) return 0.0; /* intersection */ + if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ } /* @@ -357,9 +477,11 @@ * Outside outer ring, so min distance to a ring * is the actual min distance */ - if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) return mindist; + if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) + { + return ; + } - /* * Its in the outer ring. * Have to check if its inside a hole @@ -372,14 +494,26 @@ * Its inside a hole, then the actual * distance is the min ring distance */ - return mindist; + return; } } + if (dl->thedir >0) + { + dl->d=0.0; + dl->p1.x=pt.x; + dl->p1.y=pt.y; + dl->p2.x=pt.x; + dl->p2.y=pt.y; + } + return ; /* Not in hole, so inside polygon */ - return 0.0; /* Not in hole, so inside polygon */ } -double distance2d_point_point(LWPOINT *point1, LWPOINT *point2) + + + + +void lw_dist2d_comp_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl) { POINT2D p1; POINT2D p2; @@ -387,22 +521,25 @@ getPoint2d_p(point1->point, 0, &p1); getPoint2d_p(point2->point, 0, &p2); - return distance2d_pt_pt(&p1, &p2); + lw_dist2d_comp_pt_pt(&p1, &p2,dl); + return; } -double distance2d_point_line(LWPOINT *point, LWLINE *line) +void lw_dist2d_comp_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl) { POINT2D p; POINTARRAY *pa = line->points; getPoint2d_p(point->point, 0, &p); - return distance2d_pt_ptarray(&p, pa); + lw_dist2d_comp_pt_ptarray(&p, pa, dl); + return; } -double distance2d_line_line(LWLINE *line1, LWLINE *line2) +void lw_dist2d_comp_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl) { POINTARRAY *pa1 = line1->points; POINTARRAY *pa2 = line2->points; - return distance2d_ptarray_ptarray(pa1, pa2); + lw_dist2d_comp_ptarray_ptarray(pa1, pa2, dl); + return; } /* @@ -410,21 +547,29 @@ * 2. if in the boundary, test to see if its in a hole. * if so, then return dist to hole, else return 0 (point in polygon) */ -double distance2d_point_poly(LWPOINT *point, LWPOLY *poly) + +void lw_dist2d_comp_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl) { POINT2D p; int i; getPoint2d_p(point->point, 0, &p); - LWDEBUG(2, "distance2d_point_poly called"); + LWDEBUG(2, "lw_dist2d_comp_point_poly called"); + + if (dl->thedir < 0) + { + LWDEBUG(3, "looking for maxdistance"); + lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl); + return; + } /* Return distance to outer ring if not inside it */ if ( ! pt_in_ring_2d(&p, poly->rings[0]) ) { LWDEBUG(3, " not inside outer-ring"); - - return distance2d_pt_ptarray(&p, poly->rings[0]); + lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl); + return; } /* @@ -439,75 +584,117 @@ if ( pt_in_ring_2d(&p, poly->rings[i]) ) { LWDEBUG(3, " inside an hole"); - - return distance2d_pt_ptarray(&p, poly->rings[i]); + lw_dist2d_comp_pt_ptarray(&p, poly->rings[i], dl); + return; } } LWDEBUG(3, " inside the polygon"); - - return 0.0; /* Is inside the polygon */ + if (dl->thedir >0) + { + dl->d=0.0; + dl->p1.x=p.x; + dl->p1.y=p.y; + dl->p2.x=p.x; + dl->p2.y=p.y; + } + return ; /* Is inside the polygon */ } - /* - * Brute force. - * Test to see if any rings intersect. - * If yes, dist=0. - * Test to see if one inside the other and if they are inside holes. - * Find min distance ring-to-ring. +1 if we are looking for maxdistance, just check the outer rings. +2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings +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 +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 +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. */ -double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2) +void lw_dist2d_comp_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl) { + POINT2D pt; - double mindist = -1; int i; + LWDEBUG(2, "lw_dist2d_comp_poly_poly called"); - LWDEBUG(2, "distance2d_poly_poly called"); + /*1 if we are looking for maxdistance, just check the outer rings.*/ + if (dl->thedir <0) + {lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); + return; + } - /* if poly1 inside poly2 return 0 */ - getPoint2d_p(poly1->rings[0], 0, &pt); - if ( pt_in_poly_2d(&pt, poly2) ) return 0.0; - /* if poly2 inside poly1 return 0 */ - getPoint2d_p(poly2->rings[0], 0, &pt); - if ( pt_in_poly_2d(&pt, poly1) ) return 0.0; - - LWDEBUG(3, " polys not inside each other"); - - /* - * foreach ring in Poly1 - * foreach ring in Poly2 - * if intersect, return 0 - */ - for (i=0; inrings; i++) +/* 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings +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.*/ + getPoint2d_p(poly1->rings[0], 0, &pt); + if ( !pt_in_ring_2d(&pt, poly2->rings[0])) { - int j; - for (j=0; jnrings; j++) + getPoint2d_p(poly2->rings[0], 0, &pt); + if (!pt_in_ring_2d(&pt, poly1->rings[0])) + { + lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); + return; + } + } + +/*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*/ + getPoint2d_p(poly2->rings[0], 0, &pt); + for (i=1; inrings; i++) + { + /* Inside a hole */ + if ( pt_in_ring_2d(&pt, poly1->rings[i]) ) { - double d = distance2d_ptarray_ptarray(poly1->rings[i], - poly2->rings[j]); - if ( d <= 0 ) return 0.0; + lw_dist2d_comp_ptarray_ptarray(poly1->rings[i], poly2->rings[0],dl); + return; + } + + } + +/*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*/ + getPoint2d_p(poly1->rings[0], 0, &pt); + for (i=1; inrings; i++) + { + /* Inside a hole */ + if ( pt_in_ring_2d(&pt, poly2->rings[i]) ) + { + lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[i],dl); + return; + } + + } + - /* mindist is -1 when not yet set */ - if (mindist > -1) mindist = LW_MIN(mindist, d); - else mindist = d; +/*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.*/ + getPoint2d_p(poly1->rings[0], 0, &pt); + if ( pt_in_ring_2d(&pt, poly2->rings[0])) + { + dl->d=0.0; + dl->p1.x=pt.x; + dl->p1.y=pt.y; + dl->p2.x=pt.x; + dl->p2.y=pt.y; + return ; + } - LWDEBUGF(3, " ring%i-%i dist: %f, mindist: %f", i, j, d, mindist); - } - + getPoint2d_p(poly2->rings[0], 0, &pt); + if (pt_in_ring_2d(&pt, poly1->rings[0])) + { + dl->d=0.0; + dl->p1.x=pt.x; + dl->p1.y=pt.y; + dl->p2.x=pt.x; + dl->p2.y=pt.y; + return ; } - - /* otherwise return closest approach of rings (no intersection) */ - return mindist; - + + + lwerror("You have found a bug. You could report that lw_dist2d_comp_poly_poly is not working as it should"); + } -double distance2d_line_poly(LWLINE *line, LWPOLY *poly) +void lw_dist2d_comp_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl) { - return distance2d_ptarray_poly(line->points, poly); + lw_dist2d_comp_ptarray_poly(line->points, poly, dl); + return; } - /*find the 2d length of the given POINTARRAY (even if it's 3d) */ double lwgeom_pointarray_length2d(POINTARRAY *pts) { @@ -552,7 +739,6 @@ ((frm.y - to.y)*(frm.y - to.y) ) + ((frm.z - to.z)*(frm.z - to.z) ) ); } - return dist; } @@ -645,17 +831,56 @@ } double +lwgeom_maxdistance2d_recursive(uchar *lw1, uchar *lw2) +{ + return lwgeom_maxdistance2d_recursive_tolerance( lw1, lw2, 0.0 ); +} + + +double +lwgeom_maxdistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) +{ + /*double thedist;*/ + DISTPTS thedl; + thedl.thedir = (-1); + thedl.d= -1.0; + thedl.tolerance = tolerance; + DISTPTS *dl = &thedl; + lw_dist2d_comp(lw1,lw2,dl); + return dl->d; + /*lwgeom_release((DISTPTS *)dl); + return thedist; + return 1212.00;*/ +} + +double lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) { return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 ); } + double lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) { + /*double thedist;*/ + DISTPTS thedl; + thedl.thedir = (1); + thedl.d= 99999999.0; + thedl.tolerance = tolerance; + DISTPTS *dl = &thedl; + lw_dist2d_comp( lw1,lw2,dl); + return dl->d; + /*lwgeom_release((DISTPTS *)dl); + return thedist; + return 1212.00;*/ +} + +void +lw_dist2d_comp(uchar *lw1, uchar *lw2, DISTPTS *dl) +{ LWGEOM_INSPECTED *in1, *in2; int i, j; - double mindist = -1; in1 = lwgeom_inspect(lw1); in2 = lwgeom_inspect(lw2); @@ -664,15 +889,13 @@ { uchar *g1 = lwgeom_getsubgeometry_inspected(in1, i); int t1 = lwgeom_getType(g1[0]); - double dist=tolerance; - /* Argument 1 is a multitype... recurse */ + /* it's a multitype... recurse */ if ( lwgeom_contains_subgeoms(t1) ) { - dist = lwgeom_mindistance2d_recursive_tolerance(g1, lw2, tolerance); - if ( dist <= tolerance ) return tolerance; /* can't be closer */ - if ( mindist == -1 ) mindist = dist; - else mindist = LW_MIN(dist, mindist); + lw_dist2d_comp(g1, lw2, dl); + + if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ continue; } @@ -681,38 +904,41 @@ uchar *g2 = lwgeom_getsubgeometry_inspected(in2, j); int t2 = lwgeom_getType(g2[0]); - /* Argument 2 is a multitype... recurse */ if ( lwgeom_contains_subgeoms(t2) ) { - dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); - if ( dist <= tolerance ) return tolerance; /* can't be closer */ - if ( mindist == -1 ) mindist = dist; - else mindist = LW_MIN(dist, mindist); + lw_dist2d_comp(g1, g2, dl); + if (dl->d<=dl->tolerance) + { + /*dl->d=tolerance;*/ + return; /* can't be closer */ + } continue; } - if ( t1 == POINTTYPE ) { if ( t2 == POINTTYPE ) { - dist = distance2d_point_point( - lwpoint_deserialize(g1), - lwpoint_deserialize(g2) - ); + dl->twisted=1; + lw_dist2d_comp_point_point( + lwpoint_deserialize(g1), + lwpoint_deserialize(g2), + dl); } else if ( t2 == LINETYPE ) { - dist = distance2d_point_line( - lwpoint_deserialize(g1), - lwline_deserialize(g2) - ); + dl->twisted=1; + lw_dist2d_comp_point_line( + lwpoint_deserialize(g1), + lwline_deserialize(g2), + dl); } else if ( t2 == POLYGONTYPE ) { - dist = distance2d_point_poly( - lwpoint_deserialize(g1), - lwpoly_deserialize(g2) - ); + dl->twisted=1; + lw_dist2d_comp_point_poly( + lwpoint_deserialize(g1), + lwpoly_deserialize(g2), + dl); } else { @@ -723,24 +949,28 @@ { if ( t2 == POINTTYPE ) { - dist = distance2d_point_line( - lwpoint_deserialize(g2), - lwline_deserialize(g1) - ); + dl->twisted=(-1); + lw_dist2d_comp_point_line( + lwpoint_deserialize(g2), + lwline_deserialize(g1), + dl); } else if ( t2 == LINETYPE ) { - dist = distance2d_line_line( - lwline_deserialize(g1), - lwline_deserialize(g2) - ); + dl->twisted=1; + /*lwnotice("start line");*/ + lw_dist2d_comp_line_line( + lwline_deserialize(g1), + lwline_deserialize(g2), + dl); } else if ( t2 == POLYGONTYPE ) { - dist = distance2d_line_poly( - lwline_deserialize(g1), - lwpoly_deserialize(g2) - ); + dl->twisted=1; + lw_dist2d_comp_line_poly( + lwline_deserialize(g1), + lwpoly_deserialize(g2), + dl); } else { @@ -751,55 +981,52 @@ { if ( t2 == POLYGONTYPE ) { - dist = distance2d_poly_poly( - lwpoly_deserialize(g2), - lwpoly_deserialize(g1) - ); + dl->twisted=(-1); + lw_dist2d_comp_poly_poly( + lwpoly_deserialize(g2), + lwpoly_deserialize(g1), + dl); } else if ( t2 == POINTTYPE ) { - dist = distance2d_point_poly( - lwpoint_deserialize(g2), - lwpoly_deserialize(g1) - ); + dl->twisted=(-1); + lw_dist2d_comp_point_poly( + lwpoint_deserialize(g2), + lwpoly_deserialize(g1), + dl); } else if ( t2 == LINETYPE ) { - dist = distance2d_line_poly( - lwline_deserialize(g2), - lwpoly_deserialize(g1) - ); + dl->twisted=(-1); + lw_dist2d_comp_line_poly( + lwline_deserialize(g2), + lwpoly_deserialize(g1), + dl); } else { lwerror("Unsupported geometry type: %s", lwgeom_typename(t2)); } } -// else if (lwgeom_contains_subgeoms(t1)) /* it's a multitype... recurse */ -// { -// dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); -// } + /* else if (lwgeom_contains_subgeoms(t1)) + { + lw_dist2d_comp(g1, g2, tolerance,dl); + + }*/ else { lwerror("Unsupported geometry type: %s", lwgeom_typename(t1)); } - if (mindist == -1 ) mindist = dist; - else mindist = LW_MIN(dist, mindist); - LWDEBUGF(3, "dist %d-%d: %f - mindist: %f", - i, j, dist, mindist); - if (mindist <= tolerance) return tolerance; /* can't be closer */ + if (dl->d<=dl->tolerance && dl->thedir > 0) return; /*just a check if the answer is already given*/ } } - - if (mindist<0) mindist = 0; - - return mindist; + return ; } Index: postgis/lwgeom_functions_basic.c =================================================================== --- postgis/lwgeom_functions_basic.c (revision 4310) +++ postgis/lwgeom_functions_basic.c (working copy) @@ -33,6 +33,7 @@ Datum LWGEOM_nrings(PG_FUNCTION_ARGS); Datum LWGEOM_area_polygon(PG_FUNCTION_ARGS); Datum LWGEOM_dwithin(PG_FUNCTION_ARGS); +Datum LWGEOM_dfullywithin(PG_FUNCTION_ARGS); Datum postgis_uses_stats(PG_FUNCTION_ARGS); Datum postgis_autocache_bbox(PG_FUNCTION_ARGS); Datum postgis_scripts_released(PG_FUNCTION_ARGS); @@ -43,8 +44,10 @@ Datum LWGEOM_length_linestring(PG_FUNCTION_ARGS); Datum LWGEOM_perimeter2d_poly(PG_FUNCTION_ARGS); Datum LWGEOM_perimeter_poly(PG_FUNCTION_ARGS); +Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS); Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS); -Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS); +Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS); +Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS); Datum LWGEOM_inside_circle_point(PG_FUNCTION_ARGS); Datum LWGEOM_collect(PG_FUNCTION_ARGS); Datum LWGEOM_accum(PG_FUNCTION_ARGS); @@ -1522,6 +1525,109 @@ PG_RETURN_POINTER(result); } + +PG_FUNCTION_INFO_V1(LWGEOM_shortestline2d); +Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS) +{ + + PG_LWGEOM *geom1; + PG_LWGEOM *geom2; + LWPOINT *point1; + LWPOINT *point2; + int SRID; + double x1,x2,y1,y2; + PG_LWGEOM *result=NULL; + LWPOINT *lwpoints[2]; + LWLINE *outline; + + DISTPTS thedl; + DISTPTS *dl = &thedl; + dl->d=99999999.0; + dl->thedir=(1); + dl->tolerance = 0.0; + + geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) + { + elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); + PG_RETURN_NULL(); + } + + SRID = pglwgeom_getSRID(geom1); + + lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2), dl ); + + x1=dl->p1.x; + y1=dl->p1.y; + x2=dl->p2.x; + y2=dl->p2.y; + + point1 = make_lwpoint2d(SRID, x1, y1); + point2 = make_lwpoint2d(SRID, x2, y2); + + lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1))); + lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2))); + + outline = lwline_from_lwpointarray(SRID, 2, lwpoints); + + result = pglwgeom_serialize((LWGEOM *)outline); + PG_RETURN_POINTER(result); + +} + +PG_FUNCTION_INFO_V1(LWGEOM_longestline2d); +Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS) +{ + + PG_LWGEOM *geom1; + PG_LWGEOM *geom2; + LWPOINT *point1; + LWPOINT *point2; + int SRID; + + double x1,x2,y1,y2; + PG_LWGEOM *result=NULL; + LWPOINT *lwpoints[2]; + LWLINE *outline; + DISTPTS thedl; + DISTPTS *dl = &thedl; + dl->tolerance = 0.0; + dl->d=(-1.0); + dl->thedir=(-1); + + geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) + { + elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); + PG_RETURN_NULL(); + } + + SRID = pglwgeom_getSRID(geom1); + + lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2), dl ); + + x1=dl->p1.x; + y1=dl->p1.y; + x2=dl->p2.x; + y2=dl->p2.y; + + point1 = make_lwpoint2d(SRID, x1, y1); + point2 = make_lwpoint2d(SRID, x2, y2); + + lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1))); + lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2))); + + outline = lwline_from_lwpointarray(SRID, 2, lwpoints); + + result = pglwgeom_serialize((LWGEOM *)outline); + + PG_RETURN_POINTER(result); +} + /* Minimum 2d distance between objects in geom1 and geom2. */ PG_FUNCTION_INFO_V1(LWGEOM_mindistance2d); Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS) @@ -1594,32 +1700,55 @@ PG_RETURN_BOOL(tolerance >= mindist); } -/* - * Maximum 2d distance between linestrings. - * Returns NULL if given geoms are not linestrings. - * This is very bogus (or I'm missing its meaning) - */ +PG_FUNCTION_INFO_V1(LWGEOM_dfullywithin); +Datum LWGEOM_dfullywithin(PG_FUNCTION_ARGS) +{ + PG_LWGEOM *geom1; + PG_LWGEOM *geom2; + double maxdist, tolerance; + + PROFSTART(PROF_QRUN); + + geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + tolerance = PG_GETARG_FLOAT8(2); + + if ( tolerance < 0 ) + { + elog(ERROR,"Tolerance cannot be less than zero\n"); + PG_RETURN_NULL(); + } + + if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) + { + elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); + PG_RETURN_NULL(); + } + maxdist = lwgeom_maxdistance2d_recursive_tolerance( + SERIALIZED_FORM(geom1), + SERIALIZED_FORM(geom2), + tolerance + ); + PROFSTOP(PROF_QRUN); + PROFREPORT("dist",geom1, geom2, NULL); + + PG_FREE_IF_COPY(geom1, 0); + PG_FREE_IF_COPY(geom2, 1); + + PG_RETURN_BOOL(tolerance >= maxdist); +} + PG_FUNCTION_INFO_V1(LWGEOM_maxdistance2d_linestring); Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS) { - PG_LWGEOM *geom1; PG_LWGEOM *geom2; - LWLINE *line1; - LWLINE *line2; - double maxdist = 0; - int i; + double maxdist; - elog(ERROR, "This function is unimplemented yet"); - PG_RETURN_NULL(); + PROFSTART(PROF_QRUN); geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - line1 = lwline_deserialize(SERIALIZED_FORM(geom1)); - if ( line1 == NULL ) PG_RETURN_NULL(); /* not a linestring */ - geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); - line2 = lwline_deserialize(SERIALIZED_FORM(geom2)); - if ( line2 == NULL ) PG_RETURN_NULL(); /* not a linestring */ if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) { @@ -1627,16 +1756,19 @@ PG_RETURN_NULL(); } - for (i=0; ipoints->npoints; i++) - { - POINT2D p; - double dist; - getPoint2d_p(line1->points, i, &p); - dist = distance2d_pt_ptarray(&p, line2->points); + maxdist = lwgeom_maxdistance2d_recursive(SERIALIZED_FORM(geom1), + SERIALIZED_FORM(geom2)); +/* + mindist = lwgeom_mindistance2d_recursive_tolerance( + SERIALIZED_FORM(geom1), + SERIALIZED_FORM(geom2), + tolerance + ); +*/ - if (dist > maxdist) maxdist = dist; - } + PROFSTOP(PROF_QRUN); + PROFREPORT("maxdist",geom1, geom2, NULL); PG_FREE_IF_COPY(geom1, 0); PG_FREE_IF_COPY(geom2, 1); @@ -1644,13 +1776,7 @@ PG_RETURN_FLOAT8(maxdist); } -/** - * @brief Longitude shift: - * Y remains the same - * X is converted: - * from -180..180 to 0..360 - * from 0..360 to -180..180 - */ + PG_FUNCTION_INFO_V1(LWGEOM_longitude_shift); Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS) { Index: regress/measures.sql =================================================================== --- regress/measures.sql (revision 4310) +++ regress/measures.sql (working copy) @@ -47,6 +47,107 @@ ) as foo; + + +--st_max_distance + +select 'st_max_distance_134', st_max_distance('POINT(1 2)', 'POINT(1 2)'); +select 'st_max_distance_135', st_max_distance('POINT(5 0)', 'POINT(10 12)'); + +select 'st_max_distance_136', st_max_distance('POINT(0 0)', ST_translate('POINT(0 0)', 5, 12, 0)); + +-- postgis-users/2006-May/012174.html +select 'st_max_distance_dist', st_max_distance(a,b), st_max_distance(b,a) from ( + select 'POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))'::geometry as a, + 'POLYGON((11 0, 11 10, 20 10, 20 0, 11 0), + (15 5, 15 8, 17 8, 17 5, 15 5))'::geometry as b + ) as foo; + + + +select 'distancetest1', + st_distance(a, b), + st_max_distance(a, b) + from ( + select geomfromtext('MULTILINESTRING((17 16, 16 17, 17 18, 18 17, 17 16), (28 35,29 39, 30 35))') as a, + geomfromtext('MULTIPOLYGON(((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14)), + ((33 35,33 40, 35 40, 35 35, 33 35)))') as b + ) as foo; + + + +select 'distancetest2', + st_distance(a, b), + st_max_distance(a, b) + from ( + select geomfromtext('LINESTRING(-40 -20 , 4 2)') as a, + geomfromtext('LINESTRING(-10 20, 1 -2)') as b + ) as foo; + + + +select 'distancepoly1', + st_distance(a, b), + st_max_distance(a, b) +from ( + select geomfromtext('MULTIPOLYGON(((17 16, 16 17, 17 18, 18 17, 17 16)), ((28 35,29 39, 30 35, 28 35)))') as a, + geomfromtext('MULTIPOLYGON(((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14)), + ((33 35,33 40, 35 40, 35 35, 33 35)))') as b + ) as foo; + + + +select 'distancepoly2', + st_distance(a, b), + st_max_distance(a, b) +from ( + select geomfromtext('POLYGON((17 14, 16 17, 17 18, 18 17, 17 14))') as a, + geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b + ) as foo; + + + +select 'distancepoly3', + st_distance(a, b), + st_max_distance(a, b) +from ( + select geomfromtext('POLYGON((17 16, 16 17, 17 19, 18 17, 17 16))') as a, + geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b + ) as foo; + + +select 'distancepoly4', + st_distance(a, b), + st_max_distance(a, b) +from ( + select geomfromtext('POLYGON((17 16, 16 17, 16 20, 18 20, 18 17, 17 16))') as a, + geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b + ) as foo; + + + +select 'distancepoly5', + st_distance(a, b), + st_max_distance(a, b) +from ( + select geomfromtext('POLYGON((17 12, 16 17, 17 18, 18 17, 17 12))') as a, + geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b + ) as foo; + + + + +select 'distancepoly6', + st_distance(a, b), + st_max_distance(a, b) +from ( + select geomfromtext('POLYGON((2 2, 2 3, 3 3, 3 2, 2 2))') as a, + geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b + ) as foo; + + + + -- Area of an empty polygon select 'emptyPolyArea', st_area('POLYGON EMPTY'); Index: regress/measures_expected =================================================================== --- regress/measures_expected (revision 4310) +++ regress/measures_expected (working copy) @@ -18,6 +18,18 @@ 135|13 136|13 dist|1|1 +st_max_distance_134|0 +st_max_distance_135|13 +st_max_distance_136|13 +st_max_distance_dist|22.3606797749979|22.3606797749979 +distancetest1|1|50 +distancetest2|0|50 +distancepoly1|1|30 +distancepoly2|0|26.1725046566048 +distancepoly3|0|26.9072480941474 +distancepoly4|0|28.3196045170126 +distancepoly5|0|26.1725046566048 +distancepoly6|0|32.5269119345812 emptyPolyArea|0 emptyLineArea|0 emptyPointArea|0