Skip to content

Commit 38a13ff

Browse files
committedJul 16, 2017
Make pole of inaccessibility calculation handle multipolygons
1 parent 6487fbb commit 38a13ff

File tree

2 files changed

+68
-21
lines changed

2 files changed

+68
-21
lines changed
 

‎src/core/geometry/qgsinternalgeometryengine.cpp

Lines changed: 61 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -155,25 +155,10 @@ Cell *getCentroidCell( const QgsPolygonV2 *polygon )
155155
return new Cell( x / area, y / area, 0.0, polygon );
156156
}
157157

158-
///@endcond
159-
160-
QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision, double *distanceFromBoundary ) const
158+
QgsPoint surfacePoleOfInaccessibility( const QgsSurface *surface, double precision, double &distanceFromBoundary )
161159
{
162-
if ( distanceFromBoundary )
163-
*distanceFromBoundary = DBL_MAX;
164-
165-
if ( !mGeometry || mGeometry->isEmpty() )
166-
return QgsGeometry();
167-
168-
if ( precision <= 0 )
169-
return QgsGeometry();
170-
171-
const QgsSurface *surface = dynamic_cast< const QgsSurface * >( mGeometry );
172-
if ( !surface )
173-
return QgsGeometry();
174-
175160
std::unique_ptr< QgsPolygonV2 > segmentizedPoly;
176-
const QgsPolygonV2 *polygon = dynamic_cast< const QgsPolygonV2 * >( mGeometry );
161+
const QgsPolygonV2 *polygon = dynamic_cast< const QgsPolygonV2 * >( surface );
177162
if ( !polygon )
178163
{
179164
segmentizedPoly.reset( static_cast< QgsPolygonV2 *>( surface->segmentize() ) );
@@ -187,7 +172,7 @@ QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision,
187172
double cellSize = qMin( bounds.width(), bounds.height() );
188173

189174
if ( qgsDoubleNear( cellSize, 0.0 ) )
190-
return QgsGeometry( new QgsPoint( bounds.xMinimum(), bounds.yMinimum() ) );
175+
return QgsPoint( bounds.xMinimum(), bounds.yMinimum() );
191176

192177
double h = cellSize / 2.0;
193178
std::priority_queue< Cell *, std::vector<Cell *>, GreaterThanByMax > cellQueue;
@@ -238,15 +223,70 @@ QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision,
238223
cellQueue.push( new Cell( currentCell->x + h, currentCell->y + h, h, polygon ) );
239224
}
240225

226+
distanceFromBoundary = bestCell->d;
227+
228+
return QgsPoint( bestCell->x, bestCell->y );
229+
}
230+
231+
///@endcond
232+
233+
QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision, double *distanceFromBoundary ) const
234+
{
241235
if ( distanceFromBoundary )
242-
*distanceFromBoundary = bestCell->d;
236+
*distanceFromBoundary = DBL_MAX;
243237

244-
return QgsGeometry( new QgsPoint( bestCell->x, bestCell->y ) );
238+
if ( !mGeometry || mGeometry->isEmpty() )
239+
return QgsGeometry();
240+
241+
if ( precision <= 0 )
242+
return QgsGeometry();
243+
244+
if ( const QgsGeometryCollection *gc = dynamic_cast< const QgsGeometryCollection *>( mGeometry ) )
245+
{
246+
int numGeom = gc->numGeometries();
247+
double maxDist = 0;
248+
QgsPoint bestPoint;
249+
bool found = false;
250+
for ( int i = 0; i < numGeom; ++i )
251+
{
252+
const QgsSurface *surface = dynamic_cast< const QgsSurface * >( gc->geometryN( i ) );
253+
if ( !surface )
254+
continue;
255+
256+
found = true;
257+
double dist = DBL_MAX;
258+
QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
259+
if ( dist > maxDist )
260+
{
261+
maxDist = dist;
262+
bestPoint = p;
263+
}
264+
}
265+
266+
if ( !found )
267+
return QgsGeometry();
268+
269+
if ( distanceFromBoundary )
270+
*distanceFromBoundary = maxDist;
271+
return QgsGeometry( new QgsPoint( bestPoint ) );
272+
}
273+
else
274+
{
275+
const QgsSurface *surface = dynamic_cast< const QgsSurface * >( mGeometry );
276+
if ( !surface )
277+
return QgsGeometry();
278+
279+
double dist = DBL_MAX;
280+
QgsPoint p = surfacePoleOfInaccessibility( surface, precision, dist );
281+
if ( distanceFromBoundary )
282+
*distanceFromBoundary = dist;
283+
return QgsGeometry( new QgsPoint( p ) );
284+
}
245285
}
246286

247287

248288
// helpers for orthogonalize
249-
// adapted from original code in potlach/id osm editor
289+
// adapted from original code in potlatch/id osm editor
250290

251291
bool dotProductWithinAngleTolerance( double dotProduct, double lowerThreshold, double upperThreshold )
252292
{

‎tests/src/core/testqgsgeometry.cpp

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5451,6 +5451,13 @@ void TestQgsGeometry::poleOfInaccessibility()
54515451
point = curved.poleOfInaccessibility( 0.01 ).asPoint();
54525452
QGSCOMPARENEAR( point.x(), -0.4324, 0.0001 );
54535453
QGSCOMPARENEAR( point.y(), -0.2434, 0.0001 );
5454+
5455+
// multipolygon
5456+
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)))" ) );
5457+
point = multiPoly.poleOfInaccessibility( 0.01, &distance ).asPoint();
5458+
QGSCOMPARENEAR( point.x(), 40, 0.0001 );
5459+
QGSCOMPARENEAR( point.y(), 45, 0.0001 );
5460+
QGSCOMPARENEAR( distance, 10.0, 0.00001 );
54545461
}
54555462

54565463
void TestQgsGeometry::makeValid()

0 commit comments

Comments
 (0)
Please sign in to comment.