Make pole of inaccessibility calculation handle multipolygons

This commit is contained in:
Nyall Dawson 2017-07-17 07:20:03 +10:00
parent 6487fbb2d1
commit 38a13ff5af
2 changed files with 69 additions and 22 deletions

View File

@ -155,25 +155,10 @@ Cell *getCentroidCell( const QgsPolygonV2 *polygon )
return new Cell( x / area, y / area, 0.0, polygon );
}
///@endcond
QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision, double *distanceFromBoundary ) const
QgsPoint surfacePoleOfInaccessibility( const QgsSurface *surface, double precision, double &distanceFromBoundary )
{
if ( distanceFromBoundary )
*distanceFromBoundary = DBL_MAX;
if ( !mGeometry || mGeometry->isEmpty() )
return QgsGeometry();
if ( precision <= 0 )
return QgsGeometry();
const QgsSurface *surface = dynamic_cast< const QgsSurface * >( mGeometry );
if ( !surface )
return QgsGeometry();
std::unique_ptr< QgsPolygonV2 > segmentizedPoly;
const QgsPolygonV2 *polygon = dynamic_cast< const QgsPolygonV2 * >( mGeometry );
const QgsPolygonV2 *polygon = dynamic_cast< const QgsPolygonV2 * >( surface );
if ( !polygon )
{
segmentizedPoly.reset( static_cast< QgsPolygonV2 *>( surface->segmentize() ) );
@ -187,7 +172,7 @@ QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision,
double cellSize = qMin( bounds.width(), bounds.height() );
if ( qgsDoubleNear( cellSize, 0.0 ) )
return QgsGeometry( new QgsPoint( bounds.xMinimum(), bounds.yMinimum() ) );
return QgsPoint( bounds.xMinimum(), bounds.yMinimum() );
double h = cellSize / 2.0;
std::priority_queue< Cell *, std::vector<Cell *>, GreaterThanByMax > cellQueue;
@ -238,15 +223,70 @@ QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision,
cellQueue.push( new Cell( currentCell->x + h, currentCell->y + h, h, polygon ) );
}
if ( distanceFromBoundary )
*distanceFromBoundary = bestCell->d;
distanceFromBoundary = bestCell->d;
return QgsGeometry( new QgsPoint( bestCell->x, bestCell->y ) );
return QgsPoint( bestCell->x, bestCell->y );
}
///@endcond
QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision, double *distanceFromBoundary ) const
{
if ( distanceFromBoundary )
*distanceFromBoundary = DBL_MAX;
if ( !mGeometry || mGeometry->isEmpty() )
return QgsGeometry();
if ( precision <= 0 )
return QgsGeometry();
if ( const QgsGeometryCollection *gc = dynamic_cast< const QgsGeometryCollection *>( mGeometry ) )
{
int numGeom = gc->numGeometries();
double maxDist = 0;
QgsPoint bestPoint;
bool found = false;
for ( int i = 0; i < numGeom; ++i )
{
const QgsSurface *surface = dynamic_cast< const QgsSurface * >( gc->geometryN( i ) );
if ( !surface )
continue;
found = true;
double dist = DBL_MAX;
QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
if ( dist > maxDist )
{
maxDist = dist;
bestPoint = p;
}
}
if ( !found )
return QgsGeometry();
if ( distanceFromBoundary )
*distanceFromBoundary = maxDist;
return QgsGeometry( new QgsPoint( bestPoint ) );
}
else
{
const QgsSurface *surface = dynamic_cast< const QgsSurface * >( mGeometry );
if ( !surface )
return QgsGeometry();
double dist = DBL_MAX;
QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
if ( distanceFromBoundary )
*distanceFromBoundary = dist;
return QgsGeometry( new QgsPoint( p ) );
}
}
// helpers for orthogonalize
// adapted from original code in potlach/id osm editor
// adapted from original code in potlatch/id osm editor
bool dotProductWithinAngleTolerance( double dotProduct, double lowerThreshold, double upperThreshold )
{

View File

@ -5451,6 +5451,13 @@ void TestQgsGeometry::poleOfInaccessibility()
point = curved.poleOfInaccessibility( 0.01 ).asPoint();
QGSCOMPARENEAR( point.x(), -0.4324, 0.0001 );
QGSCOMPARENEAR( point.y(), -0.2434, 0.0001 );
// multipolygon
QgsGeometry multiPoly = QgsGeometry::fromWkt( QStringLiteral( "MultiPolygon (((0 0, 10 0, 10 10, 0 10, 0 0)),((30 30, 50 30, 50 60, 30 60, 30 30)))" ) );
point = multiPoly.poleOfInaccessibility( 0.01, &distance ).asPoint();
QGSCOMPARENEAR( point.x(), 40, 0.0001 );
QGSCOMPARENEAR( point.y(), 45, 0.0001 );
QGSCOMPARENEAR( distance, 10.0, 0.00001 );
}
void TestQgsGeometry::makeValid()