[pal] Move storage of GEOS geometry from Feature up to PointSet. Use

GEOS to calculate feature centroids rather than custom code.
qgis · Jul 19, 2015 · cb249c3 · cb249c3
1 parent 9a23cec
commit cb249c3
Show file tree

Hide file tree

Showing 4 changed files with 79 additions and 80 deletions.
diff --git a/src/core/pal/feature.cpp b/src/core/pal/feature.cpp
@@ -116,11 +116,6 @@ namespace pal
     qDeleteAll( mHoles );
     mHoles.clear();
 
-    if ( mOwnsGeom )
-    {
-      GEOSGeom_destroy_r( geosContext(), mGeos );
-      mGeos = NULL;
-    }
   }
 
   void FeaturePart::extractCoords( const GEOSGeometry* geom )
@@ -1298,6 +1293,9 @@ namespace pal
 
   void FeaturePart::addSizePenalty( int nbp, LabelPosition** lPos, double bbx[4], double bby[4] )
   {
+    if ( !mGeos )
+      createGeosGeom();
+
     GEOSContextHandle_t ctxt = geosContext();
     int geomType = GEOSGeomTypeId_r( ctxt, mGeos );
 
@@ -1338,11 +1336,22 @@ namespace pal
 
   bool FeaturePart::isConnected( FeaturePart* p2 )
   {
+    if ( !mGeos )
+      createGeosGeom();
+
+    if ( !p2->mGeos )
+      p2->createGeosGeom();
+
     return ( GEOSTouches_r( geosContext(), mGeos, p2->mGeos ) == 1 );
   }
 
   bool FeaturePart::mergeWithFeaturePart( FeaturePart* other )
   {
+    if ( !mGeos )
+      createGeosGeom();
+    if ( !other->mGeos )
+      other->createGeosGeom();
+
     GEOSContextHandle_t ctxt = geosContext();
     GEOSGeometry* g1 = GEOSGeom_clone_r( ctxt, mGeos );
     GEOSGeometry* g2 = GEOSGeom_clone_r( ctxt, other->mGeos );

diff --git a/src/core/pal/feature.h b/src/core/pal/feature.h
@@ -230,10 +230,6 @@ namespace pal
        */
       Feature* getFeature() { return mFeature; }
 
-      /** Returns the part's GEOS geometry.
-       */
-      const GEOSGeometry* getGeometry() const { return mGeos; }
-
       /** Returns the layer that feature belongs to.
        */
       Layer* layer();
@@ -292,8 +288,6 @@ namespace pal
 
       Feature* mFeature;
       QList<FeaturePart*> mHoles;
-      GEOSGeometry *mGeos;
-      bool mOwnsGeom;
 
       /** \brief read coordinates from a GEOS geom */
       void extractCoords( const GEOSGeometry* geom );

diff --git a/src/core/pal/pointset.cpp b/src/core/pal/pointset.cpp
@@ -45,7 +45,9 @@ namespace pal
 
 
   PointSet::PointSet()
-      : holeOf( NULL )
+      : mGeos( 0 )
+      , mOwnsGeom( false )
+      , holeOf( NULL )
       , parent( NULL )
       , xmin( DBL_MAX )
       , xmax( -DBL_MAX )
@@ -60,7 +62,9 @@ namespace pal
   }
 
   PointSet::PointSet( int nbPoints, double *x, double *y )
-      : cHullSize( 0 )
+      : mGeos( 0 )
+      , mOwnsGeom( false )
+      , cHullSize( 0 )
       , holeOf( NULL )
       , parent( NULL )
       , xmin( DBL_MAX )
@@ -83,8 +87,12 @@ namespace pal
   }
 
   PointSet::PointSet( double aX, double aY )
-      : xmin( aX ), xmax( aY )
-      , ymin( aX ), ymax( aY )
+      : mGeos( 0 )
+      , mOwnsGeom( false )
+      , xmin( aX )
+      , xmax( aY )
+      , ymin( aX )
+      , ymax( aY )
   {
     nbPoints = cHullSize = 1;
     x = new double[1];
@@ -100,7 +108,9 @@ namespace pal
   }
 
   PointSet::PointSet( PointSet &ps )
-      : parent( 0 )
+      : mGeos( 0 )
+      , mOwnsGeom( false )
+      , parent( 0 )
       , xmin( DBL_MAX )
       , xmax( -DBL_MAX )
       , ymin( DBL_MAX )
@@ -138,8 +148,27 @@ namespace pal
     holeOf = ps.holeOf;
   }
 
+  void PointSet::createGeosGeom()
+  {
+    GEOSContextHandle_t geosctxt = geosContext();
+    GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints, 2 );
+    for ( int i = 0; i < nbPoints; ++i )
+    {
+      GEOSCoordSeq_setX_r( geosctxt, coord, i, x[i] );
+      GEOSCoordSeq_setY_r( geosctxt, coord, i, y[i] );
+    }
+    mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), 0, 0 );
+    mOwnsGeom = true;
+  }
+
   PointSet::~PointSet()
   {
+    if ( mOwnsGeom )
+    {
+      GEOSGeom_destroy_r( geosContext(), mGeos );
+      mGeos = NULL;
+    }
+
     deleteCoords();
 
     if ( cHull )
@@ -156,9 +185,6 @@ namespace pal
     y = NULL;
   }
 
-
-
-
   PointSet* PointSet::extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty )
   {
 
@@ -842,65 +868,38 @@ namespace pal
   }
 
 
-
   void PointSet::getCentroid( double &px, double &py, bool forceInside )
   {
-    // for explanation see this page:
-    // http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/
+    if ( !mGeos )
+      createGeosGeom();
 
-    int i, j;
-    double cx = 0, cy = 0, A = 0, tmp, sumx = 0, sumy = 0;
-    for ( i = 0; i < nbPoints; i++ )
-    {
-      j = i + 1; if ( j == nbPoints ) j = 0;
-      tmp = (( x[i] - x[0] ) * ( y[j] - y[0] ) - ( x[j] - x[0] ) * ( y[i] - y[0] ) );
-      cx += ( x[i] + x[j] - 2 * x[0] ) * tmp;
-      cy += ( y[i] + y[j] - 2 * y[0] ) * tmp;
-      A += tmp;
-      sumx += x[i];
-      sumy += y[i];
-    }
+    if ( !mGeos )
+      return;
 
-    if ( A == 0 )
+    GEOSContextHandle_t geosctxt = geosContext();
+    GEOSGeometry *centroidGeom = GEOSGetCentroid_r( geosctxt, mGeos );
+    if ( centroidGeom )
     {
-      px = sumx / nbPoints;
-      py = sumy / nbPoints;
-    }
-    else
-    {
-      px = cx / ( 3 * A ) + x[0];
-      py = cy / ( 3 * A ) + y[0];
+      const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom );
+      GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
+      GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
     }
 
     // check if centroid inside in polygon
     if ( forceInside && !isPointInPolygon( nbPoints, x, y, px, py ) )
     {
-      GEOSContextHandle_t geosctxt = geosContext();
-      GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints, 2 );
-
-      for ( int i = 0; i < nbPoints; ++i )
-      {
-        GEOSCoordSeq_setX_r( geosctxt, coord, i, x[i] );
-        GEOSCoordSeq_setY_r( geosctxt, coord, i, y[i] );
-      }
-
-      GEOSGeometry *geom = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), 0, 0 );
+      GEOSGeometry *pointGeom = GEOSPointOnSurface_r( geosctxt, mGeos );
 
-      if ( geom )
+      if ( pointGeom )
       {
-        GEOSGeometry *pointGeom = GEOSPointOnSurface_r( geosctxt, geom );
-
-        if ( pointGeom )
-        {
-          const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom );
-          GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
-          GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
-
-          GEOSGeom_destroy_r( geosctxt, pointGeom );
-        }
-        GEOSGeom_destroy_r( geosctxt, geom );
+        const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom );
+        GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
+        GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
+        GEOSGeom_destroy_r( geosctxt, pointGeom );
       }
     }
+
+    GEOSGeom_destroy_r( geosctxt, centroidGeom );
   }
 
 } // end namespace

diff --git a/src/core/pal/pointset.h b/src/core/pal/pointset.h
@@ -95,43 +95,36 @@ 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 );
 
       PointSet* createProblemSpecificPointSet( double bbmin[2], double bbmax[2], bool *inside );
 
       CHullBox * compute_chull_bbox();
 
 
-      /*
-       * split a concave shape into several convex shapes
-       *
+      /** 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 );
 
 
 
-      /**
-       * \brief return the minimum distance bw this and the point (px,py)
-       *
-       * compute the minimum distance bw the point (px,py) and this.
-       * Optionnaly, store the nearest point in (rx,ry)
-       *
+      /** Return the minimum distance bw this and the point (px,py). Optionally, store the nearest point 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 getDistInside(double px, double py);
-
       void getCentroid( double &px, double &py, bool forceInside = false );
 
-
       int getGeosType() const { return type; }
 
       void getBoundingBox( double min[2], double max[2] ) const
@@ -145,8 +138,7 @@ namespace pal
 
       int getNumPoints() const { return nbPoints; }
 
-      /*
-       * Iterate on line by real step of dl on x,y points
+      /** Iterate on line by real step of dl on x,y points.
        * @param nbPoint # point in line
        * @param x x coord
        * @param y y coord
@@ -197,6 +189,9 @@ namespace pal
       }
 
     protected:
+      GEOSGeometry *mGeos;
+      bool mOwnsGeom;
+
       int nbPoints;
       double *x;
       double *y;   // points order is counterclockwise
@@ -214,11 +209,13 @@ namespace pal
       PointSet( PointSet &ps );
 
       void deleteCoords();
+      void createGeosGeom();
 
       double xmin;
       double xmax;
       double ymin;
       double ymax;
+
   };
 
 } // namespace pal