[pal] Use GEOS for calculating minimum distances between geometries

Also more cleanups, fixes for geos handling

Author: nyalldawson

src/core/pal/costcalculator.cpp

@@ -227,9 +227,7 @@ namespace pal
 
   void PolygonCostCalculator::update( PointSet *pset )
   {
-    double rx, ry;
-    pset->getDist( px, py, &rx, &ry );
-    double d = dist_euc2d_sq( px, py, rx, ry );
+    double d = pset->minDistanceToPoint( px, py );
     if ( d < dist )
     {
       dist = d;

src/core/pal/feature.cpp

@@ -1365,9 +1365,11 @@ namespace pal
 
     if ( mOwnsGeom ) // delete old geometry if we own it
       GEOSGeom_destroy_r( ctxt, mGeos );
+    GEOSPreparedGeom_destroy_r( ctxt, mPreparedGeom );
     // set up new geometry
     mGeos = gTmp;
     mOwnsGeom = true;
+    mPreparedGeom = 0;
 
     deleteCoords();
     qDeleteAll( mHoles );

src/core/pal/geomfunction.cpp

@@ -145,24 +145,6 @@ namespace pal
              && cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
   }
 
-  /*
-   * \brief Intersection bw a line and a segment
-   * \return true if the point exist false otherwise
-   */
-  bool computeLineSegIntersection( double x1, double y1, double x2, double y2,  // 1st line
-                                   double x3, double y3, double x4, double y4,  // 2nd segment
-                                   double *x, double *y )
-  {
-    double cp1, cp2;
-    cp1 = cross_product( x1, y1, x2, y2, x3, y3 );
-    cp2 = cross_product( x1, y1, x2, y2, x4, y4 );
-
-    if ( cp1 * cp2 <= 0 )
-      return computeLineIntersection( x1, y1, x2, y2, x3, y3, x4, y4, x, y );
-
-    return false;
-  }
-
   /*
    * \brief compute the point wherre two lines intersects
    * \return true if the ok false if line are parallel

src/core/pal/geomfunction.h

@@ -69,14 +69,6 @@ namespace pal
   bool isSegIntersects( double x1, double y1, double x2, double y2,  // 1st segment
                         double x3, double y3, double x4, double y4 ); // 2nd segment
 
-  /*
-   * \brief Intersection bw a line and a segment
-   * \return true if the point exist false otherwise
-   */
-  bool computeLineSegIntersection( double x1, double y1, double x2, double y2,  // 1st line
-                                   double x3, double y3, double x4, double y4,  // 2nd segment
-                                   double *x, double *y );
-
   /*
    * \brief compute the point wherre two lines intersects
    * \return true if the ok false if line are parallel

src/core/pal/layer.cpp

@@ -378,7 +378,7 @@ namespace pal
     while ( !featureParts->isEmpty() )
     {
       FeaturePart* fpart = featureParts->takeFirst();
-      const GEOSGeometry* geom = fpart->getGeometry();
+      const GEOSGeometry* geom = fpart->geos();
       double chopInterval = fpart->getFeature()->repeatDistance();
       if ( chopInterval != 0. && GEOSGeomTypeId_r( geosctxt, geom ) == GEOS_LINESTRING )
       {

src/core/pal/pointset.cpp

@@ -154,13 +154,29 @@ namespace pal
   void PointSet::createGeosGeom() const
   {
     GEOSContextHandle_t geosctxt = geosContext();
-    GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints, 2 );
+
+    bool needClose = false;
+    if ( x[0] != x[ nbPoints - 1] || y[0] != y[ nbPoints - 1 ] )
+    {
+      needClose = true;
+    }
+
+    GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints + ( needClose ? 1 : 0 ), 2 );
     for ( int i = 0; i < nbPoints; ++i )
     {
       GEOSCoordSeq_setX_r( geosctxt, coord, i, x[i] );
       GEOSCoordSeq_setY_r( geosctxt, coord, i, y[i] );
     }
+
+    //close ring if needed
+    if ( needClose )
+    {
+      GEOSCoordSeq_setX_r( geosctxt, coord, nbPoints, x[0] );
+      GEOSCoordSeq_setY_r( geosctxt, coord, nbPoints, y[0] );
+    }
+
     mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), 0, 0 );
+
     mOwnsGeom = true;
   }
 
@@ -180,7 +196,7 @@ namespace pal
   {
     GEOSContextHandle_t geosctxt = geosContext();
 
-    if ( mOwnsGeom )
+    if ( mGeos && mOwnsGeom )
     {
       GEOSGeom_destroy_r( geosctxt, mGeos );
       mGeos = NULL;
@@ -625,44 +641,6 @@ namespace pal
     }
   }
 
-
-  PointSet* PointSet::createProblemSpecificPointSet( double bbmin[2], double bbmax[2], bool *inside )
-  {
-#ifdef _DEBUG_FULL_
-    std::cout << "CreateProblemSpecific:" << std::endl;
-#endif
-    PointSet *shape = new PointSet();
-    shape->x = new double[nbPoints];
-    shape->y = new double[nbPoints];
-    shape->nbPoints = nbPoints;
-    shape->type = type;
-
-    shape->xmin = xmin;
-    shape->xmax = xmax;
-    shape->ymin = ymin;
-    shape->ymax = ymax;
-
-    *inside = true;
-
-    for ( int i = 0; i < nbPoints; i++ )
-    {
-      shape->x[i] = this->x[i];
-      shape->y[i] = this->y[i];
-
-      // check whether it's not outside
-      if ( x[i] < bbmin[0] || x[i] > bbmax[0] || y[i] < bbmin[1] || y[i] > bbmax[1] )
-        *inside = false;
-    }
-
-    shape->holeOf = NULL;
-    shape->parent = NULL;
-
-    return shape;
-  }
-
-
-
-
   CHullBox * PointSet::compute_chull_bbox()
   {
     int i;
@@ -830,74 +808,47 @@ namespace pal
     return finalBb;
   }
 
-  double PointSet::getDist( double px, double py, double *rx, double *ry )
+  double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry )
   {
-    if ( nbPoints == 1 || type == GEOS_POINT )
-    {
-      if ( rx && ry )
-      {
-        *rx = x[0];
-        *ry = y[0];
-      }
-      return dist_euc2d_sq( x[0], y[0], px, py );
-    }
-
-    int a, b;
-    int nbP = ( type == GEOS_POLYGON ? nbPoints : nbPoints - 1 );
+    if ( !mGeos )
+      createGeosGeom();
 
-    double best_dist = DBL_MAX;
-    double d;
+    if ( !mGeos )
+      return 0;
 
-    for ( a = 0; a < nbP; a++ )
-    {
-      b = ( a + 1 ) % nbPoints;
-
-      double px2, py2;
-      px2 = px - y[b] + y[a];
-      py2 = py + x[b] - x[a];
-      double ix, iy;
-
-      // (px,py)->(px2,py2) is a line perpendicular to a->b
-      // Check the line p->p2 cross the segment a->b
-      if ( computeLineSegIntersection( px, py, px2, py2,
-                                       x[a], y[a], x[b], y[b],
-                                       &ix, &iy ) )
-      {
-        d = dist_euc2d_sq( px, py, ix, iy );
-      }
-      else
-      {
-        double d1 = dist_euc2d_sq( x[a], y[a], px, py );
-        double d2 = dist_euc2d_sq( x[b], y[b], px, py );
-        if ( d1 < d2 )
-        {
-          d = d1;
-          ix = x[a];
-          iy = y[a];
-        }
-        else
-        {
-          d = d2;
-          ix = x[b];
-          iy = y[b];
-        }
-      }
+    GEOSContextHandle_t geosctxt = geosContext();
 
-      if ( d < best_dist )
-      {
-        best_dist = d;
-        if ( rx && ry )
-        {
-          *rx = ix;
-          *ry = iy;
-        }
-      }
-    } // end for (a in nbPoints)
+    GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
+    GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
+    GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
+    GEOSGeometry* geosPt = GEOSGeom_createPoint_r( geosctxt, coord );
 
-    return best_dist;
+    int type = GEOSGeomTypeId_r( geosctxt, mGeos );
+    const GEOSGeometry* extRing = 0;
+    if ( type != GEOS_POLYGON )
+    {
+      extRing = mGeos;
+    }
+    else
+    {
+      //for polygons, we want distance to exterior ring (not an interior point)
+      extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
+    }
+    GEOSCoordSequence *nearestCoord = GEOSNearestPoints_r( geosctxt, extRing, geosPt );
+    double nx;
+    double ny;
+    ( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord, 0, &nx );
+    ( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord, 0, &ny );
+    GEOSGeom_destroy_r( geosctxt, geosPt );
+
+    if ( rx )
+      *rx = nx;
+    if ( ry )
+      *ry = ny;
+
+    return dist_euc2d_sq( px, py, nx, ny );
   }
 
-
   void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
   {
     if ( !mGeos )
@@ -932,6 +883,15 @@ namespace pal
     GEOSGeom_destroy_r( geosctxt, centroidGeom );
   }
 
+  const GEOSGeometry *PointSet::geos() const
+  {
+    if ( !mGeos )
+      createGeosGeom();
+
+    return mGeos;
+  }
+
+
 } // end namespace
 
 #endif

src/core/pal/pointset.h

@@ -44,18 +44,6 @@ namespace pal
   class Projection;
   class LabelPosition;
 
-  typedef struct _cross
-  {
-    int pt;
-    double d;
-    double x;
-    double y;
-    int seg;        // seg{0,1,2,3}
-    int nextCorner; // pt{0,1,2,3}
-    int way;
-
-  } Crossing;
-
   class PointSet;
 
   typedef struct _cHullBox
@@ -70,18 +58,6 @@ namespace pal
   } CHullBox;
 
 
-
-  inline bool ptrCrossingCompare( Crossing * a, Crossing * b )
-  {
-    return a == b;
-  }
-
-  inline bool crossingDist( void *a, void *b )
-  {
-    return (( Crossing* ) a )->d > (( Crossing* ) b )->d;
-  }
-
-
   class CORE_EXPORT PointSet
   {
       friend class FeaturePart;
@@ -95,10 +71,6 @@ namespace pal
       PointSet( int nbPoints, double *x, double *y );
       virtual ~PointSet();
 
-      /** Returns the point set's GEOS geometry.
-       */
-      const GEOSGeometry* getGeometry() const { return mGeos; }
-
       PointSet* extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty );
 
       /** Tests whether point set contains a specified point.
@@ -108,27 +80,23 @@ namespace pal
        */
       bool containsPoint( double x, double y ) const;
 
-      PointSet* createProblemSpecificPointSet( double bbmin[2], double bbmax[2], bool *inside );
-
       CHullBox * compute_chull_bbox();
 
-
       /** Split a concave shape into several convex shapes.
        */
       static void splitPolygons( QLinkedList<PointSet *> &shapes_toProcess,
                                  QLinkedList<PointSet *> &shapes_final,
                                  double xrm, double yrm, const QString &uid );
 
-
-
-      /** Return the minimum distance bw this and the point (px,py). Optionally, store the nearest point in (rx,ry).
+      /** Returns the minimum distance between the point set geometry and the point (px,py)
+       * Optionally, the nearest point is stored in (rx,ry).
        * @param px x coordinate of the point
        * @param py y coordinate of the points
        * @param rx pointer to x coorinates of the nearest point (can be NULL)
        * @param ry pointer to y coorinates of the nearest point (can be NULL)
        * @returns minimum distance
        */
-      double getDist( double px, double py, double *rx, double *ry );
+      double minDistanceToPoint( double px, double py, double *rx = 0, double *ry = 0 );
 
       void getCentroid( double &px, double &py, bool forceInside = false ) const;
 
@@ -151,7 +119,7 @@ namespace pal
        * @param dl distance to traverse along line
        * @param px final x coord on line
        * @param py final y coord on line
-       */
+      */
       inline void getPoint( double *d, double *ad, double dl,
                             double *px, double *py )
       {
@@ -192,6 +160,10 @@ namespace pal
         }
       }
 
+      /** Returns the point set's GEOS geometry.
+      */
+      const GEOSGeometry* geos() const;
+
     protected:
       mutable GEOSGeometry *mGeos;
       mutable bool mOwnsGeom;