[pal] Use GEOS for point in polygon test, also use prepared GEOS

geometries where possible

Author: nyalldawson

src/core/pal/costcalculator.cpp

@@ -59,7 +59,7 @@ namespace pal
         {
           case PolygonInterior:
             // n ranges from 0 -> 12
-            n = lp->getNumPointsInPolygon( obstacle->getNumPoints(), obstacle->x, obstacle->y );
+            n = lp->getNumPointsInPolygon( obstacle );
             break;
           case PolygonBoundary:
             // penalty may need tweaking, given that interior mode ranges up to 12

src/core/pal/feature.cpp

@@ -1061,8 +1061,7 @@ namespace pal
             {
               for ( ry = py, j = 0; j < 2; ry = ry + 2 * yrm, j++ )
               {
-                // TODO should test with the polyone insteand of the bbox
-                if ( !isPointInPolygon( 4, box->x, box->y, rx, ry ) )
+                if ( !mapShape->containsPoint( rx, ry ) )
                 {
                   enoughPlace = false;
                   break;
@@ -1130,7 +1129,7 @@ namespace pal
               ry += box->y[0];
 
               // Only accept candidate that center is in the polygon
-              if ( isPointInPolygon( mapShape->nbPoints, mapShape->x, mapShape->y, rx, ry ) )
+              if ( mapShape->containsPoint( rx, ry ) )
               {
                 // cost is set to minimal value, evaluated later
                 positions.append( new LabelPosition( id++, rx - dlx, ry - dly, xrm, yrm, alpha, 0.0001, this ) ); // Polygon
@@ -1336,13 +1335,10 @@ namespace pal
 
   bool FeaturePart::isConnected( FeaturePart* p2 )
   {
-    if ( !mGeos )
-      createGeosGeom();
-
     if ( !p2->mGeos )
       p2->createGeosGeom();
 
-    return ( GEOSTouches_r( geosContext(), mGeos, p2->mGeos ) == 1 );
+    return ( GEOSPreparedTouches_r( geosContext(), preparedGeom(), p2->mGeos ) == 1 );
   }
 
   bool FeaturePart::mergeWithFeaturePart( FeaturePart* other )

src/core/pal/geomfunction.cpp

@@ -378,25 +378,6 @@ namespace pal
 
   }
 
-
-  bool isPointInPolygon( int npol, double *xp, double *yp, double x, double y )
-  {
-    // code from Randolph Franklin (found at http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/)
-    int i, j;
-    bool c = false;
-
-    for ( i = 0, j = npol - 1; i < npol; j = i++ )
-    {
-      if (((( yp[i] <= y ) && ( y < yp[j] ) ) ||
-           (( yp[j] <= y ) && ( y < yp[i] ) ) )
-          && ( x < ( xp[j] - xp[i] ) *( y - yp[i] ) / ( yp[j] - yp[i] ) + xp[i] ) )
-      {
-        c = !c;
-      }
-    }
-    return c;
-  }
-
   void findLineCircleIntersection( double cx, double cy, double radius,
                                    double x1, double y1, double x2, double y2,
                                    double& xRes, double& yRes )

src/core/pal/geomfunction.h

@@ -59,8 +59,6 @@ namespace pal
     return ( x2 - x1 ) *( x2 - x1 ) + ( y2 - y1 ) *( y2 - y1 );
   }
 
-  bool isPointInPolygon( int npol, double *xp, double *yp, double x, double y );
-
   void findLineCircleIntersection( double cx, double cy, double radius,
                                    double x1, double y1, double x2, double y2,
                                    double& xRes, double& yRes );

src/core/pal/labelposition.cpp

@@ -586,20 +586,20 @@ namespace pal
     return false;
   }
 
-  int LabelPosition::getNumPointsInPolygon( int npol, double *xp, double *yp )
+  int LabelPosition::getNumPointsInPolygon( PointSet *polygon ) const
   {
     int a, k, count = 0;
     double px, py;
 
-    // cheack each corner
+    // check each corner
     for ( k = 0; k < 4; k++ )
     {
       px = x[k];
       py = y[k];
 
       for ( a = 0; a < 2; a++ ) // and each middle of segment
       {
-        if ( isPointInPolygon( npol, xp, yp, px, py ) )
+        if ( polygon->containsPoint( px, py ) )
           count++;
         px = ( x[k] + x[( k+1 ) %4] ) / 2.0;
         py = ( y[k] + y[( k+1 ) %4] ) / 2.0;
@@ -610,7 +610,7 @@ namespace pal
     py = ( y[0] + y[2] ) / 2.0;
 
     // and the label center
-    if ( isPointInPolygon( npol, xp, yp, px, py ) )
+    if ( polygon->containsPoint( px, py ) )
       count += 4; // virtually 4 points
 
     // TODO: count with nextFeature

src/core/pal/labelposition.h

@@ -132,7 +132,7 @@ namespace pal
       bool isBorderCrossingLine( PointSet* feat );
 
       /** Returns number of intersections with polygon (testing border and center) */
-      int getNumPointsInPolygon( int npol, double *xp, double *yp );
+      int getNumPointsInPolygon( PointSet* polygon ) const;
 
       /** Shift the label by specified offset */
       void offsetPosition( double xOffset, double yOffset );

src/core/pal/pointset.cpp

@@ -43,7 +43,6 @@
 namespace pal
 {
 
-
   PointSet::PointSet()
       : mGeos( 0 )
       , mOwnsGeom( false )
@@ -53,6 +52,7 @@ namespace pal
       , xmax( -DBL_MAX )
       , ymin( DBL_MAX )
       , ymax( -DBL_MAX )
+      , mPreparedGeom( 0 )
   {
     nbPoints = cHullSize =  0;
     x = NULL;
@@ -71,6 +71,7 @@ namespace pal
       , xmax( -DBL_MAX )
       , ymin( DBL_MAX )
       , ymax( -DBL_MAX )
+      , mPreparedGeom( 0 )
   {
     this->nbPoints = nbPoints;
     this->x = new double[nbPoints];
@@ -93,6 +94,7 @@ namespace pal
       , xmax( aY )
       , ymin( aX )
       , ymax( aY )
+      , mPreparedGeom( 0 )
   {
     nbPoints = cHullSize = 1;
     x = new double[1];
@@ -115,6 +117,7 @@ namespace pal
       , xmax( -DBL_MAX )
       , ymin( DBL_MAX )
       , ymax( -DBL_MAX )
+      , mPreparedGeom( 0 )
   {
     int i;
 
@@ -148,7 +151,7 @@ namespace pal
     holeOf = ps.holeOf;
   }
 
-  void PointSet::createGeosGeom()
+  void PointSet::createGeosGeom() const
   {
     GEOSContextHandle_t geosctxt = geosContext();
     GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints, 2 );
@@ -161,13 +164,28 @@ namespace pal
     mOwnsGeom = true;
   }
 
+  const GEOSPreparedGeometry *PointSet::preparedGeom() const
+  {
+    if ( !mGeos )
+      createGeosGeom();
+
+    if ( !mPreparedGeom )
+    {
+      mPreparedGeom = GEOSPrepare_r( geosContext(), mGeos );
+    }
+    return mPreparedGeom;
+  }
+
   PointSet::~PointSet()
   {
+    GEOSContextHandle_t geosctxt = geosContext();
+
     if ( mOwnsGeom )
     {
-      GEOSGeom_destroy_r( geosContext(), mGeos );
+      GEOSGeom_destroy_r( geosctxt, mGeos );
       mGeos = NULL;
     }
+    GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
 
     deleteCoords();
 
@@ -230,6 +248,18 @@ namespace pal
     return newShape;
   }
 
+  bool PointSet::containsPoint( double x, double y ) const
+  {
+    GEOSContextHandle_t geosctxt = geosContext();
+    GEOSCoordSequence* seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
+    GEOSCoordSeq_setX_r( geosctxt, seq, 0, x );
+    GEOSCoordSeq_setY_r( geosctxt, seq, 0, y );
+    GEOSGeometry* point = GEOSGeom_createPoint_r( geosctxt, seq );
+
+    bool result = ( GEOSPreparedContains_r( geosctxt, preparedGeom(), point ) == 1 );
+    GEOSGeom_destroy_r( geosctxt, point );
+    return result;
+  }
 
   void PointSet::splitPolygons( QLinkedList<PointSet*> &shapes_toProcess,
                                 QLinkedList<PointSet*> &shapes_final,
@@ -868,7 +898,7 @@ namespace pal
   }
 
 
-  void PointSet::getCentroid( double &px, double &py, bool forceInside )
+  void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
   {
     if ( !mGeos )
       createGeosGeom();
@@ -886,7 +916,7 @@ namespace pal
     }
 
     // check if centroid inside in polygon
-    if ( forceInside && !isPointInPolygon( nbPoints, x, y, px, py ) )
+    if ( forceInside && !containsPoint( px, py ) )
     {
       GEOSGeometry *pointGeom = GEOSPointOnSurface_r( geosctxt, mGeos );
 

src/core/pal/pointset.h

@@ -101,6 +101,13 @@ namespace pal
 
       PointSet* extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty );
 
+      /** Tests whether point set contains a specified point.
+       * @param x x-coordinate of point
+       * @param y y-coordinate of point
+       * @returns true if point set contains a specified point
+       */
+      bool containsPoint( double x, double y ) const;
+
       PointSet* createProblemSpecificPointSet( double bbmin[2], double bbmax[2], bool *inside );
 
       CHullBox * compute_chull_bbox();
@@ -123,7 +130,7 @@ namespace pal
        */
       double getDist( double px, double py, double *rx, double *ry );
 
-      void getCentroid( double &px, double &py, bool forceInside = false );
+      void getCentroid( double &px, double &py, bool forceInside = false ) const;
 
       int getGeosType() const { return type; }
 
@@ -189,8 +196,8 @@ namespace pal
       }
 
     protected:
-      GEOSGeometry *mGeos;
-      bool mOwnsGeom;
+      mutable GEOSGeometry *mGeos;
+      mutable bool mOwnsGeom;
 
       int nbPoints;
       double *x;
@@ -209,13 +216,18 @@ namespace pal
       PointSet( PointSet &ps );
 
       void deleteCoords();
-      void createGeosGeom();
+      void createGeosGeom() const;
+      const GEOSPreparedGeometry* preparedGeom() const;
 
       double xmin;
       double xmax;
       double ymin;
       double ymax;
 
+    private:
+
+      mutable const GEOSPreparedGeometry* mPreparedGeom;
+
   };
 
 } // namespace pal