Skip to content

Commit 42bef4e

Browse files
nyalldawsonnyalldawson
committedJul 20, 2015
[pal] Use GEOS for calculating minimum distances between geometries
Also more cleanups, fixes for geos handling
·
ltr-3_40archive/master_2
1 parent 323fa2c commit 42bef4e

File tree

7 files changed

+73
-167
lines changed

7 files changed

+73
-167
lines changed
 

‎src/core/pal/costcalculator.cpp

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -227,9 +227,7 @@ namespace pal
227227

228228
void PolygonCostCalculator::update( PointSet *pset )
229229
{
230-
double rx, ry;
231-
pset->getDist( px, py, &rx, &ry );
232-
double d = dist_euc2d_sq( px, py, rx, ry );
230+
double d = pset->minDistanceToPoint( px, py );
233231
if ( d < dist )
234232
{
235233
dist = d;

‎src/core/pal/feature.cpp

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1365,9 +1365,11 @@ namespace pal
13651365

13661366
if ( mOwnsGeom ) // delete old geometry if we own it
13671367
GEOSGeom_destroy_r( ctxt, mGeos );
1368+
GEOSPreparedGeom_destroy_r( ctxt, mPreparedGeom );
13681369
// set up new geometry
13691370
mGeos = gTmp;
13701371
mOwnsGeom = true;
1372+
mPreparedGeom = 0;
13711373

13721374
deleteCoords();
13731375
qDeleteAll( mHoles );

‎src/core/pal/geomfunction.cpp

Lines changed: 0 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -145,24 +145,6 @@ namespace pal
145145
&& cross_product( x3, y3, x4, y4, x1, y1 ) * cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
146146
}
147147

148-
/*
149-
* \brief Intersection bw a line and a segment
150-
* \return true if the point exist false otherwise
151-
*/
152-
bool computeLineSegIntersection( double x1, double y1, double x2, double y2, // 1st line
153-
double x3, double y3, double x4, double y4, // 2nd segment
154-
double *x, double *y )
155-
{
156-
double cp1, cp2;
157-
cp1 = cross_product( x1, y1, x2, y2, x3, y3 );
158-
cp2 = cross_product( x1, y1, x2, y2, x4, y4 );
159-
160-
if ( cp1 * cp2 <= 0 )
161-
return computeLineIntersection( x1, y1, x2, y2, x3, y3, x4, y4, x, y );
162-
163-
return false;
164-
}
165-
166148
/*
167149
* \brief compute the point wherre two lines intersects
168150
* \return true if the ok false if line are parallel

‎src/core/pal/geomfunction.h

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -69,14 +69,6 @@ namespace pal
6969
bool isSegIntersects( double x1, double y1, double x2, double y2, // 1st segment
7070
double x3, double y3, double x4, double y4 ); // 2nd segment
7171

72-
/*
73-
* \brief Intersection bw a line and a segment
74-
* \return true if the point exist false otherwise
75-
*/
76-
bool computeLineSegIntersection( double x1, double y1, double x2, double y2, // 1st line
77-
double x3, double y3, double x4, double y4, // 2nd segment
78-
double *x, double *y );
79-
8072
/*
8173
* \brief compute the point wherre two lines intersects
8274
* \return true if the ok false if line are parallel

‎src/core/pal/layer.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -378,7 +378,7 @@ namespace pal
378378
while ( !featureParts->isEmpty() )
379379
{
380380
FeaturePart* fpart = featureParts->takeFirst();
381-
const GEOSGeometry* geom = fpart->getGeometry();
381+
const GEOSGeometry* geom = fpart->geos();
382382
double chopInterval = fpart->getFeature()->repeatDistance();
383383
if ( chopInterval != 0. && GEOSGeomTypeId_r( geosctxt, geom ) == GEOS_LINESTRING )
384384
{

‎src/core/pal/pointset.cpp

Lines changed: 61 additions & 101 deletions
Original file line numberDiff line numberDiff line change
@@ -154,13 +154,29 @@ namespace pal
154154
void PointSet::createGeosGeom() const
155155
{
156156
GEOSContextHandle_t geosctxt = geosContext();
157-
GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints, 2 );
157+
158+
bool needClose = false;
159+
if ( x[0] != x[ nbPoints - 1] || y[0] != y[ nbPoints - 1 ] )
160+
{
161+
needClose = true;
162+
}
163+
164+
GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, nbPoints + ( needClose ? 1 : 0 ), 2 );
158165
for ( int i = 0; i < nbPoints; ++i )
159166
{
160167
GEOSCoordSeq_setX_r( geosctxt, coord, i, x[i] );
161168
GEOSCoordSeq_setY_r( geosctxt, coord, i, y[i] );
162169
}
170+
171+
//close ring if needed
172+
if ( needClose )
173+
{
174+
GEOSCoordSeq_setX_r( geosctxt, coord, nbPoints, x[0] );
175+
GEOSCoordSeq_setY_r( geosctxt, coord, nbPoints, y[0] );
176+
}
177+
163178
mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), 0, 0 );
179+
164180
mOwnsGeom = true;
165181
}
166182

@@ -180,7 +196,7 @@ namespace pal
180196
{
181197
GEOSContextHandle_t geosctxt = geosContext();
182198

183-
if ( mOwnsGeom )
199+
if ( mGeos && mOwnsGeom )
184200
{
185201
GEOSGeom_destroy_r( geosctxt, mGeos );
186202
mGeos = NULL;
@@ -625,44 +641,6 @@ namespace pal
625641
}
626642
}
627643

628-
629-
PointSet* PointSet::createProblemSpecificPointSet( double bbmin[2], double bbmax[2], bool *inside )
630-
{
631-
#ifdef _DEBUG_FULL_
632-
std::cout << "CreateProblemSpecific:" << std::endl;
633-
#endif
634-
PointSet *shape = new PointSet();
635-
shape->x = new double[nbPoints];
636-
shape->y = new double[nbPoints];
637-
shape->nbPoints = nbPoints;
638-
shape->type = type;
639-
640-
shape->xmin = xmin;
641-
shape->xmax = xmax;
642-
shape->ymin = ymin;
643-
shape->ymax = ymax;
644-
645-
*inside = true;
646-
647-
for ( int i = 0; i < nbPoints; i++ )
648-
{
649-
shape->x[i] = this->x[i];
650-
shape->y[i] = this->y[i];
651-
652-
// check whether it's not outside
653-
if ( x[i] < bbmin[0] || x[i] > bbmax[0] || y[i] < bbmin[1] || y[i] > bbmax[1] )
654-
*inside = false;
655-
}
656-
657-
shape->holeOf = NULL;
658-
shape->parent = NULL;
659-
660-
return shape;
661-
}
662-
663-
664-
665-
666644
CHullBox * PointSet::compute_chull_bbox()
667645
{
668646
int i;
@@ -830,74 +808,47 @@ namespace pal
830808
return finalBb;
831809
}
832810

833-
double PointSet::getDist( double px, double py, double *rx, double *ry )
811+
double PointSet::minDistanceToPoint( double px, double py, double *rx, double *ry )
834812
{
835-
if ( nbPoints == 1 || type == GEOS_POINT )
836-
{
837-
if ( rx && ry )
838-
{
839-
*rx = x[0];
840-
*ry = y[0];
841-
}
842-
return dist_euc2d_sq( x[0], y[0], px, py );
843-
}
844-
845-
int a, b;
846-
int nbP = ( type == GEOS_POLYGON ? nbPoints : nbPoints - 1 );
813+
if ( !mGeos )
814+
createGeosGeom();
847815

848-
double best_dist = DBL_MAX;
849-
double d;
816+
if ( !mGeos )
817+
return 0;
850818

851-
for ( a = 0; a < nbP; a++ )
852-
{
853-
b = ( a + 1 ) % nbPoints;
854-
855-
double px2, py2;
856-
px2 = px - y[b] + y[a];
857-
py2 = py + x[b] - x[a];
858-
double ix, iy;
859-
860-
// (px,py)->(px2,py2) is a line perpendicular to a->b
861-
// Check the line p->p2 cross the segment a->b
862-
if ( computeLineSegIntersection( px, py, px2, py2,
863-
x[a], y[a], x[b], y[b],
864-
&ix, &iy ) )
865-
{
866-
d = dist_euc2d_sq( px, py, ix, iy );
867-
}
868-
else
869-
{
870-
double d1 = dist_euc2d_sq( x[a], y[a], px, py );
871-
double d2 = dist_euc2d_sq( x[b], y[b], px, py );
872-
if ( d1 < d2 )
873-
{
874-
d = d1;
875-
ix = x[a];
876-
iy = y[a];
877-
}
878-
else
879-
{
880-
d = d2;
881-
ix = x[b];
882-
iy = y[b];
883-
}
884-
}
819+
GEOSContextHandle_t geosctxt = geosContext();
885820

886-
if ( d < best_dist )
887-
{
888-
best_dist = d;
889-
if ( rx && ry )
890-
{
891-
*rx = ix;
892-
*ry = iy;
893-
}
894-
}
895-
} // end for (a in nbPoints)
821+
GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
822+
GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
823+
GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
824+
GEOSGeometry* geosPt = GEOSGeom_createPoint_r( geosctxt, coord );
896825

897-
return best_dist;
826+
int type = GEOSGeomTypeId_r( geosctxt, mGeos );
827+
const GEOSGeometry* extRing = 0;
828+
if ( type != GEOS_POLYGON )
829+
{
830+
extRing = mGeos;
831+
}
832+
else
833+
{
834+
//for polygons, we want distance to exterior ring (not an interior point)
835+
extRing = GEOSGetExteriorRing_r( geosctxt, mGeos );
836+
}
837+
GEOSCoordSequence *nearestCoord = GEOSNearestPoints_r( geosctxt, extRing, geosPt );
838+
double nx;
839+
double ny;
840+
( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord, 0, &nx );
841+
( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord, 0, &ny );
842+
GEOSGeom_destroy_r( geosctxt, geosPt );
843+
844+
if ( rx )
845+
*rx = nx;
846+
if ( ry )
847+
*ry = ny;
848+
849+
return dist_euc2d_sq( px, py, nx, ny );
898850
}
899851

900-
901852
void PointSet::getCentroid( double &px, double &py, bool forceInside ) const
902853
{
903854
if ( !mGeos )
@@ -932,6 +883,15 @@ namespace pal
932883
GEOSGeom_destroy_r( geosctxt, centroidGeom );
933884
}
934885

886+
const GEOSGeometry *PointSet::geos() const
887+
{
888+
if ( !mGeos )
889+
createGeosGeom();
890+
891+
return mGeos;
892+
}
893+
894+
935895
} // end namespace
936896

937897
#endif

‎src/core/pal/pointset.h

Lines changed: 8 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -44,18 +44,6 @@ namespace pal
4444
class Projection;
4545
class LabelPosition;
4646

47-
typedef struct _cross
48-
{
49-
int pt;
50-
double d;
51-
double x;
52-
double y;
53-
int seg; // seg{0,1,2,3}
54-
int nextCorner; // pt{0,1,2,3}
55-
int way;
56-
57-
} Crossing;
58-
5947
class PointSet;
6048

6149
typedef struct _cHullBox
@@ -70,18 +58,6 @@ namespace pal
7058
} CHullBox;
7159

7260

73-
74-
inline bool ptrCrossingCompare( Crossing * a, Crossing * b )
75-
{
76-
return a == b;
77-
}
78-
79-
inline bool crossingDist( void *a, void *b )
80-
{
81-
return (( Crossing* ) a )->d > (( Crossing* ) b )->d;
82-
}
83-
84-
8561
class CORE_EXPORT PointSet
8662
{
8763
friend class FeaturePart;
@@ -95,10 +71,6 @@ namespace pal
9571
PointSet( int nbPoints, double *x, double *y );
9672
virtual ~PointSet();
9773

98-
/** Returns the point set's GEOS geometry.
99-
*/
100-
const GEOSGeometry* getGeometry() const { return mGeos; }
101-
10274
PointSet* extractShape( int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty );
10375

10476
/** Tests whether point set contains a specified point.
@@ -108,27 +80,23 @@ namespace pal
10880
*/
10981
bool containsPoint( double x, double y ) const;
11082

111-
PointSet* createProblemSpecificPointSet( double bbmin[2], double bbmax[2], bool *inside );
112-
11383
CHullBox * compute_chull_bbox();
11484

115-
11685
/** Split a concave shape into several convex shapes.
11786
*/
11887
static void splitPolygons( QLinkedList<PointSet *> &shapes_toProcess,
11988
QLinkedList<PointSet *> &shapes_final,
12089
double xrm, double yrm, const QString &uid );
12190

122-
123-
124-
/** Return the minimum distance bw this and the point (px,py). Optionally, store the nearest point in (rx,ry).
91+
/** Returns the minimum distance between the point set geometry and the point (px,py)
92+
* Optionally, the nearest point is stored in (rx,ry).
12593
* @param px x coordinate of the point
12694
* @param py y coordinate of the points
12795
* @param rx pointer to x coorinates of the nearest point (can be NULL)
12896
* @param ry pointer to y coorinates of the nearest point (can be NULL)
12997
* @returns minimum distance
13098
*/
131-
double getDist( double px, double py, double *rx, double *ry );
99+
double minDistanceToPoint( double px, double py, double *rx = 0, double *ry = 0 );
132100

133101
void getCentroid( double &px, double &py, bool forceInside = false ) const;
134102

@@ -151,7 +119,7 @@ namespace pal
151119
* @param dl distance to traverse along line
152120
* @param px final x coord on line
153121
* @param py final y coord on line
154-
*/
122+
*/
155123
inline void getPoint( double *d, double *ad, double dl,
156124
double *px, double *py )
157125
{
@@ -192,6 +160,10 @@ namespace pal
192160
}
193161
}
194162

163+
/** Returns the point set's GEOS geometry.
164+
*/
165+
const GEOSGeometry* geos() const;
166+
195167
protected:
196168
mutable GEOSGeometry *mGeos;
197169
mutable bool mOwnsGeom;

0 commit comments

Comments
 (0)
Please sign in to comment.