Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
[pal] Use GEOS for point in polygon test, also use prepared GEOS
geometries where possible
  • Loading branch information
nyalldawson committed Jul 19, 2015
1 parent cb249c3 commit 1b126d3
Show file tree
Hide file tree
Showing 8 changed files with 60 additions and 43 deletions.
2 changes: 1 addition & 1 deletion src/core/pal/costcalculator.cpp
Expand Up @@ -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
Expand Down
10 changes: 3 additions & 7 deletions src/core/pal/feature.cpp
Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 )
Expand Down
19 changes: 0 additions & 19 deletions src/core/pal/geomfunction.cpp
Expand Up @@ -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 )
Expand Down
2 changes: 0 additions & 2 deletions src/core/pal/geomfunction.h
Expand Up @@ -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 );
Expand Down
8 changes: 4 additions & 4 deletions src/core/pal/labelposition.cpp
Expand Up @@ -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;
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/core/pal/labelposition.h
Expand Up @@ -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 );
Expand Down
40 changes: 35 additions & 5 deletions src/core/pal/pointset.cpp
Expand Up @@ -43,7 +43,6 @@
namespace pal
{


PointSet::PointSet()
: mGeos( 0 )
, mOwnsGeom( false )
Expand All @@ -53,6 +52,7 @@ namespace pal
, xmax( -DBL_MAX )
, ymin( DBL_MAX )
, ymax( -DBL_MAX )
, mPreparedGeom( 0 )
{
nbPoints = cHullSize = 0;
x = NULL;
Expand All @@ -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];
Expand All @@ -93,6 +94,7 @@ namespace pal
, xmax( aY )
, ymin( aX )
, ymax( aY )
, mPreparedGeom( 0 )
{
nbPoints = cHullSize = 1;
x = new double[1];
Expand All @@ -115,6 +117,7 @@ namespace pal
, xmax( -DBL_MAX )
, ymin( DBL_MAX )
, ymax( -DBL_MAX )
, mPreparedGeom( 0 )
{
int i;

Expand Down Expand Up @@ -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 );
Expand All @@ -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();

Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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();
Expand All @@ -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 );

Expand Down
20 changes: 16 additions & 4 deletions src/core/pal/pointset.h
Expand Up @@ -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();
Expand All @@ -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; }

Expand Down Expand Up @@ -189,8 +196,8 @@ namespace pal
}

protected:
GEOSGeometry *mGeos;
bool mOwnsGeom;
mutable GEOSGeometry *mGeos;
mutable bool mOwnsGeom;

int nbPoints;
double *x;
Expand All @@ -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
Expand Down

0 comments on commit 1b126d3

Please sign in to comment.