@@ -955,21 +955,44 @@ namespace pal
955
955
956
956
void PointSet::getPointByDistance ( double distance, double *px, double *py ) const
957
957
{
958
- if ( !mGeos )
959
- createGeosGeom ();
958
+ // if anything fails, return the first point
959
+ *px = x[0 ];
960
+ *py = y[0 ];
960
961
961
- if ( !mGeos )
962
- return ;
962
+ if ( distance >= 0 )
963
+ {
964
+ // positive distance, use GEOS for interpolation
965
+ if ( !mGeos )
966
+ createGeosGeom ();
963
967
964
- GEOSContextHandle_t geosctxt = geosContext ();
965
- GEOSGeometry *point = GEOSInterpolate_r ( geosctxt, mGeos , distance );
966
- if ( point )
968
+ if ( !mGeos )
969
+ return ;
970
+
971
+ GEOSContextHandle_t geosctxt = geosContext ();
972
+
973
+ GEOSGeometry *point = GEOSInterpolate_r ( geosctxt, mGeos , distance );
974
+ if ( point )
975
+ {
976
+ const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r ( geosctxt, point );
977
+ GEOSCoordSeq_getX_r ( geosctxt, coordSeq, 0 , px );
978
+ GEOSCoordSeq_getY_r ( geosctxt, coordSeq, 0 , py );
979
+ }
980
+ GEOSGeom_destroy_r ( geosctxt, point );
981
+ }
982
+ else
967
983
{
968
- const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r ( geosctxt, point );
969
- GEOSCoordSeq_getX_r ( geosctxt, coordSeq, 0 , px );
970
- GEOSCoordSeq_getY_r ( geosctxt, coordSeq, 0 , py );
984
+ // negative distance. In this case we extrapolate backward from the first point,
985
+ // using the gradient of the entire linestring
986
+ double dx = x[nbPoints-1 ] - x[0 ];
987
+ double dy = y[nbPoints-1 ] - y[0 ];
988
+ double di = sqrt ( dx * dx + dy * dy );
989
+
990
+ if ( di != 0.0 )
991
+ {
992
+ *px = x[0 ] + distance * dx / di;
993
+ *py = y[0 ] + distance * dy / di;
994
+ }
971
995
}
972
- GEOSGeom_destroy_r ( geosctxt, point );
973
996
}
974
997
975
998
const GEOSGeometry *PointSet::geos () const
0 commit comments