Ticket #63: shortest_longest_fullywithin_new_poly_poly_20090717.patch
| File shortest_longest_fullywithin_new_poly_poly_20090717.patch, 51.8 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"/><xref linkend="ST_LongestLine"/></para> 9881 </refsection> 9882 </refentry> 9883 9884 <refentry id="ST_ShortestLine"> 9885 <refnamediv> 9886 <refname>ST_ShortestLine</refname> 9887 9888 <refpurpose>Returns the 2-dimensional shortest line between two geometries</refpurpose> 9889 </refnamediv> 9890 9891 <refsynopsisdiv> 9892 <funcsynopsis> 9893 <funcprototype> 9894 <funcdef>geometry <function>ST_ShortestLine</function></funcdef> 9895 9896 <paramdef><type>geometry </type> 9897 <parameter>g1</parameter></paramdef> 9898 9899 <paramdef><type>geometry </type> 9900 <parameter>g2</parameter></paramdef> 9901 </funcprototype> 9902 </funcsynopsis> 9903 </refsynopsisdiv> 9904 9905 <refsection> 9906 <title>Description</title> 9907 9908 <para>Returns the 2-dimensional shortest line between two geometries. The function will 9909 only return the first shortest line if more than one, that the function finds. 9910 If g1 and g2 intersects in just one point the function will return a line with both start 9911 and end in that intersection-point. 9912 If g1 and g2 are intersecting with more than one point the function will return a line with start 9913 and end in the same point but it can be any of the intersecting points. 9914 The line returned will always start in g1 and end in g2. 9915 The length of the line this function returns will always be the same as st_distance returns for g1 and g2. 9916 </para> 9917 9918 </refsection> 9919 9920 <refsection> 9921 <title>Examples</title> 9922 9923 <programlisting>postgis=# SELECT astext(st_shortestline('POINT(0 0)'::geometry, 'LINESTRING ( 2 0, 0 2 )'::geometry)); 9924 st_shortestline 9925 ----------------- 9926 "LINESTRING(0 0,1 1)" 9927 (1 row)</programlisting> 9928 </refsection> 9929 9930 <refsection> 9931 <title>See Also</title> 9932 9933 <para><xref linkend="ST_Distance"/><xref linkend="ST_LongestLine"/></para> 9934 </refsection> 9935 </refentry> 9936 9937 <refentry id="ST_LongestLine"> 9938 <refnamediv> 9939 <refname>ST_LongestLine</refname> 9940 9941 <refpurpose>Returns the 2-dimensional longest line points of two geometries. 9942 The function will only return the first longest line if more than one, that the function finds. 9943 The line returned will always start in g1 and end in g2. 9944 The length of the line this function returns will always be the same as st_max_distance returns for g1 and g2.</refpurpose> 9945 </refnamediv> 9946 9947 <refsynopsisdiv> 9948 <funcsynopsis> 9949 <funcprototype> 9950 <funcdef>geometry <function>ST_LongestLine</function></funcdef> 9951 9952 <paramdef><type>geometry </type> 9953 <parameter>g1</parameter></paramdef> 9954 9955 <paramdef><type>geometry </type> 9956 <parameter>g2</parameter></paramdef> 9957 </funcprototype> 9958 </funcsynopsis> 9959 </refsynopsisdiv> 9960 9961 <refsection> 9962 <title>Description</title> 9963 9964 <para>Returns the 2-dimensional longest line between the points of two geometries. 9965 </para> 9966 9967 </refsection> 9968 9969 <refsection> 9970 <title>Examples</title> 9971 9972 <programlisting>postgis=# SELECT astext(st_longestline('POINT(0 0)'::geometry, 'LINESTRING ( 2 0, 0 2 )'::geometry)); 9973 st_shortestline 9974 ----------------- 9975 "LINESTRING(0 0,2 0)" 9976 (1 row)</programlisting> 9977 </refsection> 9978 9979 <refsection> 9980 <title>See Also</title> 9981 9982 <para><xref linkend="ST_Max_Distance"/><xref linkend="ST_ShortestLine"/></para> 9983 </refsection> 9984 </refentry> 9985 9836 9986 <refentry id="ST_HausdorffDistance"> 9837 9987 <refnamediv> 9838 9988 <refname>ST_HausdorffDistance</refname> … … 10059 10209 <para>Returns true if the geometries are within the specified distance 10060 10210 of one another. The distance is specified in units defined by the 10061 10211 spatial reference system of the geometries. For this function to make 10062 sense, the source geometries must both be of the same coor indate projection,10212 sense, the source geometries must both be of the same coordinate projection, 10063 10213 having the same SRID.</para> 10064 10214 10065 10215 <note> … … 10109 10259 <refsection> 10110 10260 <title>See Also</title> 10111 10261 10112 <para><xref linkend="ST_Distance"/>, <xref linkend="ST_Expand"/> </para>10262 <para><xref linkend="ST_Distance"/>, <xref linkend="ST_Expand"/>, <xref linkend="ST_DFullyWithin"/></para> 10113 10263 </refsection> 10114 10264 </refentry> 10115 10265 10266 <refentry id="ST_DFullyWithin"> 10267 <refnamediv> 10268 <refname>ST_DFullyWithin</refname> 10269 10270 <refpurpose>Returns true if all of the geometries are within the specified 10271 distance of one another</refpurpose> 10272 </refnamediv> 10273 10274 <refsynopsisdiv> 10275 <funcsynopsis> 10276 <funcprototype> 10277 <funcdef>boolean <function>ST_DFullyWithin</function></funcdef> 10278 10279 <paramdef><type>geometry </type> 10280 <parameter>g1</parameter></paramdef> 10281 10282 <paramdef><type>geometry </type> 10283 <parameter>g2</parameter></paramdef> 10284 10285 <paramdef><type>double precision </type> 10286 <parameter>distance</parameter></paramdef> 10287 </funcprototype> 10288 </funcsynopsis> 10289 </refsynopsisdiv> 10290 10291 <refsection> 10292 <title>Description</title> 10293 10294 <para>Returns true if the geometries fully within the specified distance 10295 of one another. The distance is specified in units defined by the 10296 spatial reference system of the geometries. For this function to make 10297 sense, the source geometries must both be of the same coordinate projection, 10298 having the same SRID.</para> 10299 10300 <note> 10301 <para>This function call will automatically include a bounding box 10302 comparison that will make use of any indexes that are available on 10303 the geometries.</para> 10304 </note> 10305 10306 10307 </refsection> 10308 10309 <refsection> 10310 <title>Examples</title> 10311 <programlisting>postgis=# SELECT ST_DFullyWithin(geom_a, geom_b, 10) as FullyWithin10, ST_DWithin(geom_a, geom_b, 10) as Within10, ST_DFullyWithin(geom_a, geom_b, 20) as FullyWithin20 from 10312 (select ST_GeomFromText('POINT(1 1)') as geom_a,ST_GeomFromText('LINESTRING(1 5, 2 7, 1 9, 14 12)') as geom_b) t1; 10313 10314 ----------------- 10315 FullyWithin10 | Within10 | FullyWithin20 | 10316 ---------------+----------+---------------+ 10317 f | t | t | </programlisting> 10318 </refsection> 10319 10320 <refsection> 10321 <title>See Also</title> 10322 10323 <para><xref linkend="ST_Max_Distance"/>, <xref linkend="ST_DWithin"/></para> 10324 </refsection> 10325 </refentry> 10326 10116 10327 <refentry id="ST_Equals"> 10117 10328 <refnamediv> 10118 10329 <refname>ST_Equals</refname> … … 10684 10895 </refsection> 10685 10896 </refentry> 10686 10897 10687 10688 10689 10690 <refentry id="ST_Max_Distance">10691 <refnamediv>10692 <refname>ST_Max_Distance</refname>10693 10694 <refpurpose>Returns the 2-dimensional largest distance between two geometries in10695 projected units.</refpurpose>10696 </refnamediv>10697 10698 <refsynopsisdiv>10699 <funcsynopsis>10700 <funcprototype>10701 <funcdef>float <function>ST_Max_Distance</function></funcdef>10702 10703 <paramdef><type>geometry </type>10704 <parameter>g1</parameter></paramdef>10705 10706 <paramdef><type>geometry </type>10707 <parameter>g2</parameter></paramdef>10708 </funcprototype>10709 </funcsynopsis>10710 </refsynopsisdiv>10711 10712 <refsection>10713 <title>Description</title>10714 10715 <para>Returns the 2-dimensional maximum cartesian distance between two linestrings in10716 projected units.</para>10717 10718 </refsection>10719 10720 <refsection>10721 <title>Examples</title>10722 10723 <programlisting>--ALL EXAMPLES current throw NOT YET IMPLEMENTED</programlisting>10724 </refsection>10725 10726 <refsection>10727 <title>See Also</title>10728 10729 <para><xref linkend="ST_Distance"/></para>10730 </refsection>10731 </refentry>10732 10733 10734 10898 <refentry id="ST_OrderingEquals"> 10735 10899 <refnamediv> 10736 10900 <refname>ST_OrderingEquals</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 } DISTPTS; 293 303 /******************************************************************/ 294 304 295 305 /* … … 1184 1194 extern void lwgeom_force3dz_recursive(uchar *serialized, uchar *optr, size_t *retsize); 1185 1195 extern void lwgeom_force3dm_recursive(uchar *serialized, uchar *optr, size_t *retsize); 1186 1196 extern void lwgeom_force4d_recursive(uchar *serialized, uchar *optr, size_t *retsize); 1197 1187 1198 extern double distance2d_pt_pt(POINT2D *p1, POINT2D *p2); 1199 extern void lw_dist2d_comp_pt_pt(POINT2D *p1, POINT2D *p2, DISTPTS *dl); 1188 1200 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); 1201 extern void lw_dist2d_comp_pt_seg(POINT2D *p, POINT2D *A, POINT2D *B, DISTPTS *dl); 1202 extern void lw_dist2d_comp_seg_seg(POINT2D *A, POINT2D *B, POINT2D *C, POINT2D *D, DISTPTS *dl); 1203 extern void lw_dist2d_comp_pt_ptarray(POINT2D *p, POINTARRAY *pa, DISTPTS *dl); 1204 extern void lw_dist2d_comp_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2, DISTPTS *dl); 1192 1205 extern int pt_in_ring_2d(POINT2D *p, POINTARRAY *ring); 1193 1206 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);1207 extern void lw_dist2d_comp_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl); 1208 extern void lw_dist2d_comp_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl); 1209 extern void lw_dist2d_comp_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl); 1210 extern void lw_dist2d_comp_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl); 1211 extern void lw_dist2d_comp_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl); 1212 extern void lw_dist2d_comp_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl); 1213 extern void lw_dist2d_comp_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl); 1201 1214 extern int azimuth_pt_pt(POINT2D *p1, POINT2D *p2, double *ret); 1202 1215 extern double lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2); 1203 1216 extern double lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance); 1217 extern double lwgeom_maxdistance2d_recursive(uchar *lw1, uchar *lw2); 1218 extern double lwgeom_maxdistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance); 1219 extern void lw_dist2d_comp(uchar *lw1, uchar *lw2, double tolerance, DISTPTS *dl); 1204 1220 extern int lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad); 1205 1221 extern int32 lwgeom_npoints(uchar *serialized); 1206 1222 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; 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); 375 lw_dist2d_comp_pt_seg(p, &start, &end,dl); 376 start = end; 244 377 245 if ( result == 0 ) return 0;246 247 start = end;248 378 } 249 379 250 return result;380 return; 251 381 } 252 382 253 383 /* test each segment of l1 against each segment of l2. Return min */ 254 double distance2d_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2)384 void lw_dist2d_comp_ptarray_ptarray(POINTARRAY *l1, POINTARRAY *l2,DISTPTS *dl) 255 385 { 256 double result = 99999999999.9; 257 char result_okay = 0; /*result is a valid min */ 386 /*lwnotice("lw_dist2d_comp_ptarray_ptarray");*/ 258 387 int t,u; 259 388 POINT2D start, end; 260 389 POINT2D start2, end2; 261 262 LWDEBUGF(2, " distance2d_ptarray_ptarray called (points: %d-%d)",390 int twist = dl->twisted; 391 LWDEBUGF(2, "lw_dist2d_comp_ptarray_ptarray called (points: %d-%d)", 263 392 l1->npoints, l2->npoints); 264 393 265 394 getPoint2d_p(l1, 0, &start); 266 395 for (t=1; t<l1->npoints; t++) /*for each segment in L1 */ 267 396 { 268 397 getPoint2d_p(l1, t, &end); 269 270 398 getPoint2d_p(l2, 0, &start2); 271 399 for (u=1; u<l2->npoints; u++) /*for each segment in L2 */ 272 400 { 273 double dist;274 275 401 getPoint2d_p(l2, u, &end2); 402 dl->twisted=twist; 403 lw_dist2d_comp_seg_seg(&start, &end, &start2, &end2,dl); 276 404 277 dist = distance2d_seg_seg(&start, &end, &start2, &end2);278 405 406 279 407 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 408 LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f", 290 409 t, u, dist, result); 291 410 292 if (result <= 0) return 0; /*intersection */293 294 411 start2 = end2; 295 412 } 296 413 start = end; 297 414 } 298 299 return result; 415 return ; 300 416 } 301 417 418 302 419 /* true if point is in poly (and not in its holes) */ 303 420 int pt_in_poly_2d(POINT2D *p, LWPOLY *poly) 304 421 { … … 317 434 return 1; /* In outer ring, not in holes */ 318 435 } 319 436 437 438 320 439 /* 321 440 * Brute force. 322 441 * Test line-ring distance against each ring. … … 327 446 * otherwise return min distance to a ring (could be outside 328 447 * polygon or inside a hole) 329 448 */ 330 double distance2d_ptarray_poly(POINTARRAY *pa, LWPOLY *poly) 449 450 void lw_dist2d_comp_ptarray_poly(POINTARRAY *pa, LWPOLY *poly, DISTPTS *dl) 331 451 { 332 452 POINT2D pt; 333 453 int i; 334 double mindist = 0;335 454 336 LWDEBUGF(2, "distance2d_ptarray_poly called (%d rings)", poly->nrings);337 455 456 LWDEBUGF(2, "lw_dist2d_comp_ptarray_poly called (%d rings)", poly->nrings); 457 458 338 459 for (i=0; i<poly->nrings; i++) 339 460 { 340 double dist = distance2d_ptarray_ptarray(pa, poly->rings[i]); 341 if (i) mindist = LW_MIN(mindist, dist); 342 else mindist = dist; 461 lw_dist2d_comp_ptarray_ptarray(pa, poly->rings[i], dl); 343 462 344 463 LWDEBUGF(3, " distance from ring %d: %f, mindist: %f", 345 464 i, dist, mindist); 346 347 if ( mindist <= 0 ) return 0.0; /* intersection */348 465 } 349 466 350 467 /* … … 357 474 * Outside outer ring, so min distance to a ring 358 475 * is the actual min distance 359 476 */ 360 if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) return mindist; 477 if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) 478 { 479 return ; 480 } 361 481 362 363 482 /* 364 483 * Its in the outer ring. 365 484 * Have to check if its inside a hole … … 372 491 * Its inside a hole, then the actual 373 492 * distance is the min ring distance 374 493 */ 375 return mindist;494 return; 376 495 } 377 496 } 497 if (dl->thedir >0) 498 { 499 dl->d=0; 500 dl->p1.x=pt.x; 501 dl->p1.y=pt.y; 502 dl->p2.x=pt.x; 503 dl->p2.y=pt.y; 504 } 505 return ; /* Not in hole, so inside polygon */ 378 506 379 return 0.0; /* Not in hole, so inside polygon */380 507 } 381 508 382 double distance2d_point_point(LWPOINT *point1, LWPOINT *point2) 509 510 511 512 513 void lw_dist2d_comp_point_point(LWPOINT *point1, LWPOINT *point2, DISTPTS *dl) 383 514 { 384 515 POINT2D p1; 385 516 POINT2D p2; … … 387 518 getPoint2d_p(point1->point, 0, &p1); 388 519 getPoint2d_p(point2->point, 0, &p2); 389 520 390 return distance2d_pt_pt(&p1, &p2); 521 lw_dist2d_comp_pt_pt(&p1, &p2,dl); 522 return; 391 523 } 392 524 393 double distance2d_point_line(LWPOINT *point, LWLINE *line)525 void lw_dist2d_comp_point_line(LWPOINT *point, LWLINE *line, DISTPTS *dl) 394 526 { 395 527 POINT2D p; 396 528 POINTARRAY *pa = line->points; 397 529 getPoint2d_p(point->point, 0, &p); 398 return distance2d_pt_ptarray(&p, pa); 530 lw_dist2d_comp_pt_ptarray(&p, pa, dl); 531 return; 399 532 } 400 533 401 double distance2d_line_line(LWLINE *line1, LWLINE *line2)534 void lw_dist2d_comp_line_line(LWLINE *line1, LWLINE *line2, DISTPTS *dl) 402 535 { 403 536 POINTARRAY *pa1 = line1->points; 404 537 POINTARRAY *pa2 = line2->points; 405 return distance2d_ptarray_ptarray(pa1, pa2); 538 lw_dist2d_comp_ptarray_ptarray(pa1, pa2, dl); 539 return; 406 540 } 407 541 408 542 /* … … 410 544 * 2. if in the boundary, test to see if its in a hole. 411 545 * if so, then return dist to hole, else return 0 (point in polygon) 412 546 */ 413 double distance2d_point_poly(LWPOINT *point, LWPOLY *poly) 547 548 void lw_dist2d_comp_point_poly(LWPOINT *point, LWPOLY *poly, DISTPTS *dl) 414 549 { 415 550 POINT2D p; 416 551 int i; 417 552 418 553 getPoint2d_p(point->point, 0, &p); 419 554 420 LWDEBUG(2, " distance2d_point_poly called");555 LWDEBUG(2, "lw_dist2d_comp_point_poly called"); 421 556 557 558 if (dl->thedir < 0) 559 { 560 LWDEBUG(3, "looking for maxdistance"); 561 lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl); 562 return; 563 } 422 564 /* Return distance to outer ring if not inside it */ 423 565 if ( ! pt_in_ring_2d(&p, poly->rings[0]) ) 424 566 { 425 567 LWDEBUG(3, " not inside outer-ring"); 426 427 return distance2d_pt_ptarray(&p, poly->rings[0]);568 lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl); 569 return; 428 570 } 429 571 430 572 /* … … 439 581 if ( pt_in_ring_2d(&p, poly->rings[i]) ) 440 582 { 441 583 LWDEBUG(3, " inside an hole"); 442 443 return distance2d_pt_ptarray(&p, poly->rings[i]);584 lw_dist2d_comp_pt_ptarray(&p, poly->rings[i], dl); 585 return; 444 586 } 445 587 } 446 588 447 589 LWDEBUG(3, " inside the polygon"); 448 449 return 0.0; /* Is inside the polygon */ 590 if (dl->thedir >0) 591 { 592 dl->d=0.0; 593 dl->p1.x=p.x; 594 dl->p1.y=p.y; 595 dl->p2.x=p.x; 596 dl->p2.y=p.y; 597 } 598 return ; /* Is inside the polygon */ 450 599 } 451 452 600 /* 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.601 1 if we are looking for maxdistance, just check the outer rings. 602 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings 603 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 604 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 605 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 606 */ 459 double distance2d_poly_poly(LWPOLY *poly1, LWPOLY *poly2)607 void lw_dist2d_comp_poly_poly(LWPOLY *poly1, LWPOLY *poly2, DISTPTS *dl) 460 608 { 609 461 610 POINT2D pt; 462 double mindist = -1;611 BOX2DFLOAT4 box1, box2; 463 612 int i; 613 LWDEBUG(2, "lw_dist2d_comp_poly_poly called"); 464 614 465 LWDEBUG(2, "distance2d_poly_poly called"); 615 /*1 if we are looking for maxdistance, just check the outer rings.*/ 616 if (dl->thedir <0) 617 {lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); 618 return; 619 } 466 620 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 621 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++) 622 /* 2 check if poly1 has first point outside poly2 and vice versa, if so, just check outer rings*/ 623 getPoint2d_p(poly1->rings[0], 0, &pt); 624 if ( !pt_in_ring_2d(&pt, poly2->rings[0])) 483 625 { 484 int j; 485 for (j=0; j<poly2->nrings; j++) 626 getPoint2d_p(poly2->rings[0], 0, &pt); 627 if (!pt_in_ring_2d(&pt, poly1->rings[0])) 628 { 629 lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[0],dl); 630 return; 631 } 632 } 633 634 /*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*/ 635 getPoint2d_p(poly2->rings[0], 0, &pt); 636 for (i=1; i<poly1->nrings; i++) 637 { 638 /* Inside a hole */ 639 if ( pt_in_ring_2d(&pt, poly1->rings[i]) ) 486 640 { 487 double d = distance2d_ptarray_ptarray(poly1->rings[i], 488 poly2->rings[j]); 489 if ( d <= 0 ) return 0.0; 641 lw_dist2d_comp_ptarray_ptarray(poly1->rings[i], poly2->rings[0],dl); 642 return; 643 } 644 645 } 646 647 /*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*/ 648 getPoint2d_p(poly1->rings[0], 0, &pt); 649 for (i=1; i<poly2->nrings; i++) 650 { 651 /* Inside a hole */ 652 if ( pt_in_ring_2d(&pt, poly2->rings[i]) ) 653 { 654 lw_dist2d_comp_ptarray_ptarray(poly1->rings[0], poly2->rings[i],dl); 655 return; 656 } 657 658 } 659 490 660 491 /* mindist is -1 when not yet set */ 492 if (mindist > -1) mindist = LW_MIN(mindist, d); 493 else mindist = d; 661 /*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.*/ 662 getPoint2d_p(poly1->rings[0], 0, &pt); 663 if ( pt_in_ring_2d(&pt, poly2->rings[0])) 664 { 665 dl->d=0.0; 666 dl->p1.x=pt.x; 667 dl->p1.y=pt.y; 668 dl->p2.x=pt.x; 669 dl->p2.y=pt.y; 670 lwnotice("nr 1 p1.x: %f",pt.x ); 671 return ; 672 } 494 673 495 LWDEBUGF(3, " ring%i-%i dist: %f, mindist: %f", i, j, d, mindist); 496 } 497 674 getPoint2d_p(poly2->rings[0], 0, &pt); 675 if (pt_in_ring_2d(&pt, poly1->rings[0])) 676 { 677 dl->d=0.0; 678 dl->p1.x=pt.x; 679 dl->p1.y=pt.y; 680 dl->p2.x=pt.x; 681 dl->p2.y=pt.y; 682 lwnotice("nr 2 p1.x: %f",pt.x ); 683 return ; 498 684 } 499 500 /* otherwise return closest approach of rings (no intersection) */501 return mindist;502 685 686 687 lwerror("You have found a bug. You could report that lw_dist2d_comp_poly_poly is not working as it should"); 688 503 689 } 504 690 505 double distance2d_line_poly(LWLINE *line, LWPOLY *poly)691 void lw_dist2d_comp_line_poly(LWLINE *line, LWPOLY *poly, DISTPTS *dl) 506 692 { 507 return distance2d_ptarray_poly(line->points, poly); 693 lw_dist2d_comp_ptarray_poly(line->points, poly, dl); 694 return; 508 695 } 509 696 510 511 697 /*find the 2d length of the given POINTARRAY (even if it's 3d) */ 512 698 double lwgeom_pointarray_length2d(POINTARRAY *pts) 513 699 { … … 552 738 ((frm.y - to.y)*(frm.y - to.y) ) + 553 739 ((frm.z - to.z)*(frm.z - to.z) ) ); 554 740 } 555 556 741 return dist; 557 742 } 558 743 … … 645 830 } 646 831 647 832 double 833 lwgeom_maxdistance2d_recursive(uchar *lw1, uchar *lw2) 834 { 835 return lwgeom_maxdistance2d_recursive_tolerance( lw1, lw2, 99999999.0 ); 836 } 837 838 839 double 840 lwgeom_maxdistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) 841 { 842 /*double thedist;*/ 843 DISTPTS thedl; 844 thedl.thedir = (-1); 845 thedl.d= -1.0; 846 DISTPTS *dl = &thedl; 847 lw_dist2d_comp( lw1, lw2,tolerance, dl ); 848 return dl->d; 849 /*lwgeom_release((DISTPTS *)dl); 850 return thedist; 851 return 1212.00;*/ 852 } 853 854 double 648 855 lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2) 649 856 { 650 857 return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 ); 651 858 } 652 859 860 653 861 double 654 862 lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance) 655 863 { 864 /*double thedist;*/ 865 DISTPTS thedl; 866 thedl.thedir = (1); 867 thedl.d= 99999999.0; 868 DISTPTS *dl = &thedl; 869 lw_dist2d_comp( lw1, lw2,tolerance, dl ); 870 return dl->d; 871 /*lwgeom_release((DISTPTS *)dl); 872 return thedist; 873 return 1212.00;*/ 874 } 875 876 void 877 lw_dist2d_comp(uchar *lw1, uchar *lw2, double tolerance, DISTPTS *dl) 878 { 656 879 LWGEOM_INSPECTED *in1, *in2; 657 880 int i, j; 658 double mindist = -1;659 881 660 882 in1 = lwgeom_inspect(lw1); 661 883 in2 = lwgeom_inspect(lw2); … … 664 886 { 665 887 uchar *g1 = lwgeom_getsubgeometry_inspected(in1, i); 666 888 int t1 = lwgeom_getType(g1[0]); 667 double dist=tolerance;668 889 669 /* Argument 1 is a multitype... recurse */890 /* it's a multitype... recurse */ 670 891 if ( lwgeom_contains_subgeoms(t1) ) 671 892 { 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); 893 lw_dist2d_comp(g1, lw2, tolerance, dl); 894 if ( ((tolerance - dl->d)*dl->thedir) >0 ) 895 { 896 /*dl->d=tolerance;*/ 897 return; /* can't be closer */ 898 } 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, tolerance, dl); 910 if ( ((tolerance - dl->d)*dl->thedir) >0 ) 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 ( ((tolerance - dl->d)*dl->thedir) >0 ) 1025 { 1026 /*dl->d=tolerance;*/ 1027 return; /* can't be closer */ 1028 } 795 1029 796 1030 } 797 1031 798 1032 } 799 800 if (mindist<0) mindist = 0; 801 802 return mindist; 1033 return ; 803 1034 } 804 1035 805 1036 -
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 1548 geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1549 geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1550 1551 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1552 { 1553 elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); 1554 PG_RETURN_NULL(); 1555 } 1556 1557 SRID = pglwgeom_getSRID(geom1); 1558 1559 lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2),0.0, dl ); 1560 1561 x1=dl->p1.x; 1562 y1=dl->p1.y; 1563 x2=dl->p2.x; 1564 y2=dl->p2.y; 1565 1566 point1 = make_lwpoint2d(SRID, x1, y1); 1567 point2 = make_lwpoint2d(SRID, x2, y2); 1568 1569 lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1))); 1570 lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2))); 1571 1572 outline = lwline_from_lwpointarray(SRID, 2, lwpoints); 1573 1574 result = pglwgeom_serialize((LWGEOM *)outline); 1575 PG_RETURN_POINTER(result); 1576 1577 } 1578 1579 PG_FUNCTION_INFO_V1(LWGEOM_longestline2d); 1580 Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS) 1581 { 1582 1583 PG_LWGEOM *geom1; 1584 PG_LWGEOM *geom2; 1585 LWPOINT *point1; 1586 LWPOINT *point2; 1587 int SRID; 1588 1589 double x1,x2,y1,y2; 1590 PG_LWGEOM *result=NULL; 1591 LWPOINT *lwpoints[2]; 1592 LWLINE *outline; 1593 DISTPTS thedl; 1594 DISTPTS *dl = &thedl; 1595 dl->d=(-1.0); 1596 dl->thedir=(-1); 1597 1598 geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1599 geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1600 1601 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1602 { 1603 elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); 1604 PG_RETURN_NULL(); 1605 } 1606 1607 SRID = pglwgeom_getSRID(geom1); 1608 1609 lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2),99999999.0, dl ); 1610 1611 x1=dl->p1.x; 1612 y1=dl->p1.y; 1613 x2=dl->p2.x; 1614 y2=dl->p2.y; 1615 1616 point1 = make_lwpoint2d(SRID, x1, y1); 1617 point2 = make_lwpoint2d(SRID, x2, y2); 1618 1619 lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1))); 1620 lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2))); 1621 1622 outline = lwline_from_lwpointarray(SRID, 2, lwpoints); 1623 1624 result = pglwgeom_serialize((LWGEOM *)outline); 1625 1626 PG_RETURN_POINTER(result); 1627 } 1628 1525 1629 /* Minimum 2d distance between objects in geom1 and geom2. */ 1526 1630 PG_FUNCTION_INFO_V1(LWGEOM_mindistance2d); 1527 1631 Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS) … … 1594 1698 PG_RETURN_BOOL(tolerance >= mindist); 1595 1699 } 1596 1700 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 */ 1701 PG_FUNCTION_INFO_V1(LWGEOM_dfullywithin); 1702 Datum LWGEOM_dfullywithin(PG_FUNCTION_ARGS) 1703 { 1704 PG_LWGEOM *geom1; 1705 PG_LWGEOM *geom2; 1706 double maxdist, tolerance; 1707 1708 PROFSTART(PROF_QRUN); 1709 1710 geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); 1711 geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); 1712 tolerance = PG_GETARG_FLOAT8(2); 1713 1714 if ( tolerance < 0 ) 1715 { 1716 elog(ERROR,"Tolerance cannot be less than zero\n"); 1717 PG_RETURN_NULL(); 1718 } 1719 1720 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1721 { 1722 elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n"); 1723 PG_RETURN_NULL(); 1724 } 1725 maxdist = lwgeom_maxdistance2d_recursive_tolerance( 1726 SERIALIZED_FORM(geom1), 1727 SERIALIZED_FORM(geom2), 1728 tolerance 1729 ); 1730 PROFSTOP(PROF_QRUN); 1731 PROFREPORT("dist",geom1, geom2, NULL); 1732 1733 PG_FREE_IF_COPY(geom1, 0); 1734 PG_FREE_IF_COPY(geom2, 1); 1735 1736 PG_RETURN_BOOL(tolerance >= maxdist); 1737 } 1738 1602 1739 PG_FUNCTION_INFO_V1(LWGEOM_maxdistance2d_linestring); 1603 1740 Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS) 1604 1741 { 1605 1606 1742 PG_LWGEOM *geom1; 1607 1743 PG_LWGEOM *geom2; 1608 LWLINE *line1; 1609 LWLINE *line2; 1610 double maxdist = 0; 1611 int i; 1744 double maxdist; 1612 1745 1613 elog(ERROR, "This function is unimplemented yet"); 1614 PG_RETURN_NULL(); 1746 PROFSTART(PROF_QRUN); 1615 1747 1616 1748 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 1749 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 1750 1624 1751 if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2)) 1625 1752 { … … 1627 1754 PG_RETURN_NULL(); 1628 1755 } 1629 1756 1630 for (i=0; i<line1->points->npoints; i++)1631 {1632 POINT2D p;1633 double dist;1634 1757 1635 getPoint2d_p(line1->points, i, &p); 1636 dist = distance2d_pt_ptarray(&p, line2->points); 1758 maxdist = lwgeom_maxdistance2d_recursive(SERIALIZED_FORM(geom1), 1759 SERIALIZED_FORM(geom2)); 1760 /* 1761 mindist = lwgeom_mindistance2d_recursive_tolerance( 1762 SERIALIZED_FORM(geom1), 1763 SERIALIZED_FORM(geom2), 1764 tolerance 1765 ); 1766 */ 1637 1767 1638 if (dist > maxdist) maxdist = dist;1639 }1768 PROFSTOP(PROF_QRUN); 1769 PROFREPORT("maxdist",geom1, geom2, NULL); 1640 1770 1641 1771 PG_FREE_IF_COPY(geom1, 0); 1642 1772 PG_FREE_IF_COPY(geom2, 1); … … 1644 1774 PG_RETURN_FLOAT8(maxdist); 1645 1775 } 1646 1776 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 */ 1777 1654 1778 PG_FUNCTION_INFO_V1(LWGEOM_longitude_shift); 1655 1779 Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS) 1656 1780 { -
postgis/postgis.sql.in.c
1446 1446 AS 'MODULE_PATHNAME', 'LWGEOM_maxdistance2d_linestring' 1447 1447 LANGUAGE 'C' IMMUTABLE STRICT; 1448 1448 1449 1450 --Two new functions 1451 1452 CREATE OR REPLACE FUNCTION ST_ShortestLine(geometry,geometry) 1453 RETURNS geometry 1454 AS 'MODULE_PATHNAME', 'LWGEOM_shortestline2d' 1455 LANGUAGE 'C' IMMUTABLE STRICT; 1456 1457 CREATE OR REPLACE FUNCTION ST_LongestLine(geometry,geometry) 1458 RETURNS geometry 1459 AS 'MODULE_PATHNAME', 'LWGEOM_longestline2d' 1460 LANGUAGE 'C' IMMUTABLE STRICT; 1461 1449 1462 -- Deprecation in 1.2.3 1450 1463 CREATE OR REPLACE FUNCTION point_inside_circle(geometry,float8,float8,float8) 1451 1464 RETURNS bool … … 4317 4330 AS 'SELECT $1 && ST_Expand($2,$3) AND $2 && ST_Expand($1,$3) AND _ST_DWithin($1, $2, $3)' 4318 4331 LANGUAGE 'SQL' IMMUTABLE; 4319 4332 4333 CREATE OR REPLACE FUNCTION _ST_DFullyWithin(geometry,geometry,float8) 4334 RETURNS boolean 4335 AS 'MODULE_PATHNAME', 'LWGEOM_dfullywithin' 4336 LANGUAGE 'C' IMMUTABLE STRICT; 4337 4338 CREATE OR REPLACE FUNCTION ST_DFullyWithin(geometry, geometry, float8) 4339 RETURNS boolean 4340 AS 'SELECT $1 && ST_Expand($2,$3) AND $2 && ST_Expand($1,$3) AND _ST_DFullyWithin($1, $2, $3)' 4341 LANGUAGE 'SQL' IMMUTABLE; 4342 4320 4343 -- Deprecation in 1.2.3 4321 4344 CREATE OR REPLACE FUNCTION intersects(geometry,geometry) 4322 4345 RETURNS boolean -
postgis/uninstall_postgis.sql.in.c
26 26 27 27 DROP FUNCTION ST_MinimumBoundingCircle(geometry); 28 28 DROP FUNCTION ST_MinimumBoundingCircle(inputgeom geometry, segs_per_quarter integer); 29 29 DROP FUNCTION ST_ShortestLine(geometry, geometry); 30 DROP FUNCTION ST_LongestLine(geometry, geometry); 31 DROP FUNCTION ST_DFullyWithin(geometry, geometry, float8); 30 32 --------------------------------------------------------------- 31 33 -- SQL-MM 32 34 --------------------------------------------------------------- -
regress/measures.sql
47 47 ) as foo; 48 48 49 49 50 --st_shortestline 51 52 select 'st_shortestline_134', st_astext(st_shortestline('POINT(1 2)', 'POINT(1 2)')); 53 select 'st_shortestline_135', st_astext(st_shortestline('POINT(5 0)', 'POINT(10 12)')); 54 55 select 'st_shortestline_136', st_astext(st_shortestline('POINT(0 0)', ST_translate('POINT(0 0)', 5, 12, 0))); 56 57 -- postgis-users/2006-May/012174.html 58 select 'st_shortestline_dist', st_astext(st_shortestline(a,b)), st_astext(st_shortestline(b,a)) from ( 59 select 'POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))'::geometry as a, 60 'POLYGON((11 0, 11 10, 20 10, 20 0, 11 0), 61 (15 5, 15 8, 17 8, 17 5, 15 5))'::geometry as b 62 ) as foo; 63 64 65 --st_max_distance 66 67 select 'st_max_distance_134', st_max_distance('POINT(1 2)', 'POINT(1 2)'); 68 select 'st_max_distance_135', st_max_distance('POINT(5 0)', 'POINT(10 12)'); 69 70 select 'st_max_distance_136', st_max_distance('POINT(0 0)', ST_translate('POINT(0 0)', 5, 12, 0)); 71 72 -- postgis-users/2006-May/012174.html 73 select 'st_max_distance_dist', st_max_distance(a,b), st_max_distance(b,a) from ( 74 select 'POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))'::geometry as a, 75 'POLYGON((11 0, 11 10, 20 10, 20 0, 11 0), 76 (15 5, 15 8, 17 8, 17 5, 15 5))'::geometry as b 77 ) as foo; 78 79 80 81 --st_longestline 82 83 select 'st_longestline_134', st_astext(st_longestline('POINT(1 2)', 'POINT(1 2)')); 84 select 'st_longestline_135', st_astext(st_longestline('POINT(5 0)', 'POINT(10 12)')); 85 86 select 'st_longestline_136', st_astext(st_longestline('POINT(0 0)', ST_translate('POINT(0 0)', 5, 12, 0))); 87 88 -- postgis-users/2006-May/012174.html 89 select 'st_longestline_dist', st_astext(st_longestline(a,b)), st_astext(st_longestline(b,a)) from ( 90 select 'POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))'::geometry as a, 91 'POLYGON((11 0, 11 10, 20 10, 20 0, 11 0), 92 (15 5, 15 8, 17 8, 17 5, 15 5))'::geometry as b 93 ) as foo; 94 95 96 97 98 99 100 select 'distancetest1', 101 st_distance(a, b), 102 st_max_distance(a, b), 103 st_astext(st_shortestline(a,b)), 104 st_astext(st_shortestline(b,a)), 105 st_astext(st_longestline(a,b)), 106 st_astext(st_longestline(b,a)) from ( 107 select geomfromtext('MULTILINESTRING((17 16, 16 17, 17 18, 18 17, 17 16), (28 35,29 39, 30 35))') as a, 108 geomfromtext('MULTIPOLYGON(((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14)), 109 ((33 35,33 40, 35 40, 35 35, 33 35)))') as b 110 ) as foo; 111 112 113 114 select 'distancetest2', 115 st_distance(a, b), 116 st_max_distance(a, b), 117 st_astext(st_shortestline(a,b)), 118 st_astext(st_shortestline(b,a)), 119 st_astext(st_longestline(a,b)), 120 st_astext(st_longestline(b,a)) from ( 121 select geomfromtext('LINESTRING(-40 -20 , 4 2)') as a, 122 geomfromtext('LINESTRING(-10 20, 1 -2)') as b 123 ) as foo; 124 125 126 127 select 'distancepoly1', 128 st_distance(a, b), 129 st_max_distance(a, b), 130 st_astext(st_shortestline(a,b)), 131 st_astext(st_shortestline(b,a)), 132 st_astext(st_longestline(a,b)), 133 st_astext(st_longestline(b,a)) from ( 134 select geomfromtext('MULTIPOLYGON(((17 16, 16 17, 17 18, 18 17, 17 16), (28 35,29 39, 30 35, 28 35)))') as a, 135 geomfromtext('MULTIPOLYGON(((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14)), 136 ((33 35,33 40, 35 40, 35 35, 33 35)))') as b 137 ) as foo; 138 139 140 141 select 'distancepoly2', 142 st_distance(a, b), 143 st_max_distance(a, b), 144 st_astext(st_shortestline(a,b)), 145 st_astext(st_shortestline(b,a)), 146 st_astext(st_longestline(a,b)), 147 st_astext(st_longestline(b,a)) from ( 148 select geomfromtext('POLYGON((17 14, 16 17, 17 18, 18 17, 17 14))') as a, 149 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 150 ) as foo; 151 152 153 154 select 'distancepoly3', 155 st_distance(a, b), 156 st_max_distance(a, b), 157 st_astext(st_shortestline(a,b)), 158 st_astext(st_shortestline(b,a)), 159 st_astext(st_longestline(a,b)), 160 st_astext(st_longestline(b,a)) from ( 161 select geomfromtext('POLYGON((17 16, 16 17, 17 19, 18 17, 17 16))') as a, 162 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 163 ) as foo; 164 165 166 select 'distancepoly4', 167 st_distance(a, b), 168 st_max_distance(a, b), 169 st_astext(st_shortestline(a,b)), 170 st_astext(st_shortestline(b,a)), 171 st_astext(st_longestline(a,b)), 172 st_astext(st_longestline(b,a)) from ( 173 select geomfromtext('POLYGON((17 16, 16 17, 16 20, 18 20, 18 17, 17 16))') as a, 174 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 175 ) as foo; 176 177 178 179 select 'distancepoly5', 180 st_distance(a, b), 181 st_max_distance(a, b), 182 st_astext(st_shortestline(a,b)), 183 st_astext(st_shortestline(b,a)), 184 st_astext(st_longestline(a,b)), 185 st_astext(st_longestline(b,a)) from ( 186 select geomfromtext('POLYGON((17 12, 16 17, 17 18, 18 17, 17 12))') as a, 187 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 188 ) as foo; 189 190 191 192 193 select 'distancepoly6', 194 st_distance(a, b), 195 st_max_distance(a, b), 196 st_astext(st_shortestline(a,b)), 197 st_astext(st_shortestline(b,a)), 198 st_astext(st_longestline(a,b)), 199 st_astext(st_longestline(b,a)) from ( 200 select geomfromtext('POLYGON((2 2, 2 3, 3 3, 3 2, 2 2))') as a, 201 geomfromtext('POLYGON((-1 -1, -1 25, 25 25, 25 -1, -1 -1), (14 14,14 19,19 19,19 14,14 14))') as b 202 ) as foo; 203 204 205 206 50 207 -- Area of an empty polygon 51 208 select 'emptyPolyArea', st_area('POLYGON EMPTY'); 52 209 -
regress/measures_expected
18 18 135|13 19 19 136|13 20 20 dist|1|1 21 st_shortestline_134|LINESTRING(1 2,1 2) 22 st_shortestline_135|LINESTRING(5 0,10 12) 23 st_shortestline_136|LINESTRING(0 0,5 12) 24 st_shortestline_dist|LINESTRING(10 10,11 10)|LINESTRING(11 10,10 10) 25 st_max_distance_134|0 26 st_max_distance_135|13 27 st_max_distance_136|13 28 st_max_distance_dist|22.3606797749979|22.3606797749979 29 st_longestline_134|LINESTRING(1 2,1 2) 30 st_longestline_135|LINESTRING(5 0,10 12) 31 st_longestline_136|LINESTRING(0 0,5 12) 32 st_longestline_dist|LINESTRING(0 0,20 10)|LINESTRING(20 10,0 0) 33 distancetest1|1|50|LINESTRING(17 19,17 18)|LINESTRING(17 18,17 19)|LINESTRING(29 39,-1 -1)|LINESTRING(-1 -1,29 39) 34 distancetest2|0|50|LINESTRING(0 0,0 0)|LINESTRING(0 0,0 0)|LINESTRING(-40 -20,-10 20)|LINESTRING(-10 20,-40 -20) 35 distancepoly1|1|30|LINESTRING(17 18,17 19)|LINESTRING(17 19,17 18)|LINESTRING(17 16,35 40)|LINESTRING(35 40,17 16) 36 distancepoly2|0|26.1725046566048|LINESTRING(17 14,17 14)|LINESTRING(17 14,17 14)|LINESTRING(17 18,-1 -1)|LINESTRING(-1 -1,17 18) 37 distancepoly3|0|26.9072480941474|LINESTRING(17 19,17 19)|LINESTRING(17 19,17 19)|LINESTRING(17 19,-1 -1)|LINESTRING(-1 -1,17 19) 38 distancepoly4|0|28.3196045170126|LINESTRING(18 19,18 19)|LINESTRING(18 19,18 19)|LINESTRING(18 20,-1 -1)|LINESTRING(-1 -1,18 20) 39 distancepoly5|0|26.1725046566048|LINESTRING(17 12,17 12)|LINESTRING(17 12,17 12)|LINESTRING(17 18,-1 -1)|LINESTRING(-1 -1,17 18) 40 distancepoly6|0|32.5269119345812|LINESTRING(2 2,2 2)|LINESTRING(2 2,2 2)|LINESTRING(2 2,25 25)|LINESTRING(25 25,2 2) 21 41 emptyPolyArea|0 22 42 emptyLineArea|0 23 43 emptyPointArea|0
