[PostGIS] #6130: Geography distance is 0 for two lines where an end point of one is the antipode of an end point of the other

[PostGIS]

#6130: Geography distance is 0 for two lines where an end point of one is the
antipode of an end point of the other
----------------------+---------------------------
 Reporter:  ezimanyi  |      Owner:  pramsey
     Type:  defect    |     Status:  new
 Priority:  medium    |  Milestone:  PostGIS 3.6.5
Component:  postgis   |    Version:  3.6.x
 Keywords:            |
----------------------+---------------------------
 ST_Distance on geography answers 0, and ST_Intersects and ST_DWithin
 answer true, for two lines
 that are 6654 km apart, whenever an end point of one line is the antipode
 of an end point of the
 other. Moving that end point by 0.001 degree gives the expected distance.

 Reproduction, PostGIS 3.6.4 (94d984b), PostgreSQL 18, GEOS 3.12.1:

 {{{
 SELECT ST_Distance(geography 'LINESTRING(-90 0,0 0)', geography
 'LINESTRING(90 0,0 60)');
 -- 0
 SELECT ST_Distance(geography 'LINESTRING(-90 0,0 0)', geography
 'LINESTRING(90 0,0 60)', false);
 -- 0
 SELECT ST_Intersects(geography 'LINESTRING(-90 0,0 0)', geography
 'LINESTRING(90 0,0 60)');
 -- t
 SELECT ST_DWithin(geography 'LINESTRING(-90 0,0 0)', geography
 'LINESTRING(90 0,0 60)', 1000);
 -- t
 SELECT ST_AsText(ST_ShortestLine(geography 'LINESTRING(-90 0,0 0)',
 geography 'LINESTRING(90 0,0 60)'));
 -- LINESTRING(-90 0,-90 0)

 -- The same with the first line starting at -89.999:
 SELECT ST_Distance(geography 'LINESTRING(-89.999 0,0 0)', geography
 'LINESTRING(90 0,0 60)');
 -- 6654072.819...
 -- The closest pair is (0 0) and (0 60):
 SELECT ST_Distance(geography 'POINT(0 0)', geography 'POINT(0 60)');
 -- 6654072.819...
 }}}

 Cause: edge_intersects() in liblwgeom/lwgeodetic.c. After ruling out the
 "no intersection" and
 "straddle" cases, it returns PIR_INTERSECTS for every remaining case, that
 is, whenever an end
 point of one edge lies on the plane of the other edge (side == 0). Such a
 point lies on the line
 where the two great circles meet, and the other edge contains either that
 point or its antipode.
 Here (-90 0), the start of the first edge, lies on the great circle of the
 edge (90 0)-(0 60) but
 is the antipode of that edge's start, so the edges do not meet, yet the
 touch is reported. The
 straddle case checks point_in_cone() before reporting an intersection; the
 touch cases do not.

 Proposed fix: count a touch only where the touching end point lies on the
 side of the centre of
 the sphere that holds the other edge. An edge is shorter than half a great
 circle, so that side
 contains the meeting point and not its antipode. The test is the sign of
 dot(P, B1 + B2); unlike
 point_in_cone(), whose strict comparison can reject a vertex that nearly
 coincides with the other
 edge's end point, it stays robust at a shared or nearly shared vertex.

 Measured with the patch applied to a vendored copy of liblwgeom
 (MobilityDB), whose
 edge_intersects() and point_in_cone() are identical to master's: the
 geography distance of the
 reproduction goes from 0 to 6654072.819 m, the value of the moved line and
 of the closest point
 pair, while edges that share an end point, cross, or touch at an end point
 still answer 0. I
 have not run the PostGIS regression suite with it.

 Patch against master d28cf14b5b:

 {{{
 --- a/liblwgeom/lwgeodetic.c
 +++ b/liblwgeom/lwgeodetic.c
 @@ -3414,6 +3414,18 @@
  }

  /**
 +* Utility function for edge_intersects(), true if P lies on the side of
 the
 +* center of the sphere that holds the edge A1/A2.
 +*/
 +static int
 +edge_faces_point(const POINT3D *A1, const POINT3D *A2, const POINT3D *P)
 +{
 +       POINT3D AC;
 +       vector_sum(A1, A2, &AC);
 +       return dot_product(P, &AC) > 0.0;
 +}
 +
 +/**
  * Returns non-zero if edges A and B interact. The type of interaction is
 given in the
  * return value with the bitmask elements defined above.
  */
 @@ -3493,30 +3505,36 @@
                 return PIR_NO_INTERACT;
         }

 -       /* The rest are all intersects variants... */
 -       rv |= PIR_INTERSECTS;
 +       /* The rest are all touch variants. An end point on the plane of
 the */
 +       /* other edge lies where the two great circles meet, and the other
 */
 +       /* edge, shorter than half a circle, holds either that point or
 its */
 +       /* antipode: it touches only in the first case. */

         /* A touches B */
 -       if ( a1_side == 0 )
 +       if ( a1_side == 0 && edge_faces_point(B1, B2, A1) )
         {
                 /* Touches at A1, A2 is on what side? */
 +               rv |= PIR_INTERSECTS;
                 rv |= (a2_side < 0 ? PIR_A_TOUCH_RIGHT :
 PIR_A_TOUCH_LEFT);
         }
 -       else if ( a2_side == 0 )
 +       else if ( a2_side == 0 && edge_faces_point(B1, B2, A2) )
         {
                 /* Touches at A2, A1 is on what side? */
 +               rv |= PIR_INTERSECTS;
                 rv |= (a1_side < 0 ? PIR_A_TOUCH_RIGHT :
 PIR_A_TOUCH_LEFT);
         }

         /* B touches A */
 -       if ( b1_side == 0 )
 +       if ( b1_side == 0 && edge_faces_point(A1, A2, B1) )
         {
                 /* Touches at B1, B2 is on what side? */
 +               rv |= PIR_INTERSECTS;
                 rv |= (b2_side < 0 ? PIR_B_TOUCH_RIGHT :
 PIR_B_TOUCH_LEFT);
         }
 -       else if ( b2_side == 0 )
 +       else if ( b2_side == 0 && edge_faces_point(A1, A2, B2) )
         {
                 /* Touches at B2, B1 is on what side? */
 +               rv |= PIR_INTERSECTS;
                 rv |= (b1_side < 0 ? PIR_B_TOUCH_RIGHT :
 PIR_B_TOUCH_LEFT);
         }

 }}}
-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/6130>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.