[pal] Use GEOS for linear referencing

src/core/pal/feature.cpp

@@ -500,9 +500,6 @@ namespace pal
     double xrm = mFeature->label_x;
     double yrm = mFeature->label_y;
 
-    double *d; // segments lengths distance bw pt[i] && pt[i+1]
-    double *ad;  // absolute distance bw pt[0] and pt[i] along the line
-    double ll; // line length
     double dist;
     double bx, by, ex, ey;
     int nbls;
@@ -528,24 +525,7 @@ namespace pal
     x = line->x;
     y = line->y;
 
-    d = new double[nbPoints-1];
-    ad = new double[nbPoints];
-
-    ll = 0.0; // line length
-    for ( i = 0; i < line->nbPoints - 1; i++ )
-    {
-      if ( i == 0 )
-        ad[i] = 0;
-      else
-        ad[i] = ad[i-1] + d[i-1];
-
-      d[i] = dist_euc2d( x[i], y[i], x[i+1], y[i+1] );
-      ll += d[i];
-    }
-
-    ad[line->nbPoints-1] = ll;
-
-
+    double ll = line->length(); // line length
     nbls = ( int )( ll / xrm ); // ratio bw line length and label width
 
 #ifdef _DEBUG_FULL_
@@ -580,9 +560,9 @@ namespace pal
     while ( l < ll - xrm )
     {
       // => bx, by
-      line->getPoint( d, ad, l, &bx, &by );
+      line->getPointByDistance( l, &bx, &by );
       // same but l = l+xrm
-      line->getPoint( d, ad, l + xrm, &ex, &ey );
+      line->getPointByDistance( l + xrm, &ex, &ey );
 
       // Label is bigger than line ...
       if ( l < 0 )
@@ -649,11 +629,6 @@ namespace pal
         break;
     }
 
-    //delete line;
-
-    delete[] d;
-    delete[] ad;
-
     int nbp = positions.size();
     *lPos = new LabelPosition *[nbp];
     i = 0;

src/core/pal/pointset.cpp

@@ -942,6 +942,25 @@ namespace pal
     GEOSGeom_destroy_r( geosctxt, centroidGeom );
   }
 
+  void PointSet::getPointByDistance( double distance, double *px, double *py ) const
+  {
+    if ( !mGeos )
+      createGeosGeom();
+
+    if ( !mGeos )
+      return;
+
+    GEOSContextHandle_t geosctxt = geosContext();
+    GEOSGeometry *point = GEOSInterpolate_r( geosctxt, mGeos, distance );
+    if ( point )
+    {
+      const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, point );
+      GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, px );
+      GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, py );
+    }
+    GEOSGeom_destroy_r( geosctxt, point );
+  }
+
   const GEOSGeometry *PointSet::geos() const
   {
     if ( !mGeos )
@@ -950,6 +969,21 @@ namespace pal
     return mGeos;
   }
 
+  double PointSet::length() const
+  {
+    if ( !mGeos )
+      createGeosGeom();
+
+    if ( !mGeos )
+      return -1;
+
+    GEOSContextHandle_t geosctxt = geosContext();
+
+    double len = 0;
+    ( void )GEOSLength_r( geosctxt, mGeos, &len );
+    return len;
+  }
+
 
 } // end namespace
 

src/core/pal/pointset.h

@@ -123,57 +123,21 @@ namespace pal
 
       int getNumPoints() const { return nbPoints; }
 
-      /** Iterate on line by real step of dl on x,y points.
-       * @param d array of distances between points
-       * @param ad absolute distance from pt0 to each point (ad0 = pt0->pt0)
-       * @param dl distance to traverse along line
+      /** Get a point a set distance along a line geometry.
+       * @param distance 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 )
-      {
-        int i;
-        double dx, dy, di;
-        double distr;
-
-        i = 0;
-        if ( dl >= 0 )
-        {
-          while ( i < nbPoints && ad[i] <= dl ) i++;
-          i--;
-        }
-
-        if ( i < nbPoints - 1 )
-        {
-          if ( dl < 0 )
-          {
-            dx = x[nbPoints-1] - x[0];
-            dy = y[nbPoints-1] - y[0];
-            di = sqrt( dx * dx + dy * dy );
-          }
-          else
-          {
-            dx = x[i+1] - x[i];
-            dy = y[i+1] - y[i];
-            di = d[i];
-          }
-
-          distr = dl - ad[i];
-          *px = x[i] + dx * distr / di;
-          *py = y[i] + dy * distr / di;
-        }
-        else    // just select last point...
-        {
-          *px = x[i];
-          *py = y[i];
-        }
-      }
+      void getPointByDistance( double distance, double *px, double *py ) const;
 
       /** Returns the point set's GEOS geometry.
       */
       const GEOSGeometry* geos() const;
 
+      /** Returns length of line geometry.
+       */
+      double length() const;
+
     protected:
       mutable GEOSGeometry *mGeos;
       mutable bool mOwnsGeom;