Skip to content

Commit d6f09c0

Browse files
nyalldawsonnyalldawson
committedNov 14, 2016
[FEATURE] Add method to calculate pole of inaccessibility for polygons
Implements a method in QgsGeometry and a processing algorithm to calculate the pole of inaccessibility for a surface, which is the most distant internal point from the boundary of the surface. This function uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative approach guaranteed to find the true pole of inaccessibility within a specified tolerance. More precise tolerances require more iterations and will take longer to calculate.
·
ltr-3_40final-3_0_0
1 parent d5c307e commit d6f09c0

File tree

13 files changed

+538
-7
lines changed

13 files changed

+538
-7
lines changed
 

‎python/core/geometry/qgsgeometry.sip

Lines changed: 24 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -510,15 +510,37 @@ class QgsGeometry
510510
/** Returns a simplified version of this geometry using a specified tolerance value */
511511
QgsGeometry simplify( double tolerance ) const;
512512

513-
/** Returns the center of mass of a geometry
513+
/**
514+
* Returns the center of mass of a geometry.
514515
* @note for line based geometries, the center point of the line is returned,
515516
* and for point based geometries, the point itself is returned
517+
* @see pointOnSurface()
518+
* @see poleOfInaccessibility()
516519
*/
517520
QgsGeometry centroid() const;
518521

519-
/** Returns a point within a geometry */
522+
/**
523+
* Returns a point guaranteed to lie on the surface of a geometry. While the centroid()
524+
* of a geometry may be located outside of the geometry itself (eg for concave shapes),
525+
* the point on surface will always be inside the geometry.
526+
* @see centroid()
527+
* @see poleOfInaccessibility()
528+
*/
520529
QgsGeometry pointOnSurface() const;
521530

531+
/**
532+
* Calculates the approximate pole of inaccessibility for a surface, which is the
533+
* most distant internal point from the boundary of the surface. This function
534+
* uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative
535+
* approach guaranteed to find the true pole of inaccessibility within a specified
536+
* tolerance. More precise tolerances require more iterations and will take longer
537+
* to calculate.
538+
* @see centroid()
539+
* @see pointOnSurface()
540+
* @note added in QGIS 3.0
541+
*/
542+
QgsGeometry poleOfInaccessibility( double precision ) const;
543+
522544
/** Returns the smallest convex polygon that contains all the points in the geometry. */
523545
QgsGeometry convexHull() const;
524546

‎python/plugins/processing/algs/help/qgis.yaml

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -337,6 +337,9 @@ qgis:polarplot: >
337337

338338
Two fields must be entered as parameters: one that define the category each feature two (to group features) and another one with the variable to plot (this has to be a numeric one)
339339

340+
qgis:poleofinaccessibility: >
341+
This algorithm calculates the pole of inaccessibility for a polygon layer, which is the most distant internal point from the boundary of the surface. This algorithm uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative approach guaranteed to find the true pole of inaccessibility within a specified tolerance (in layer units). More precise tolerances require more iterations and will take longer to calculate.
342+
340343
qgis:polygoncentroids: >
341344
This algorithm creates a new point layer, with points representing the centroid of polygons of an input layer.
342345

‎python/plugins/processing/algs/qgis/PointOnSurface.py

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,9 @@ class PointOnSurface(GeoAlgorithm):
4545
INPUT_LAYER = 'INPUT_LAYER'
4646
OUTPUT_LAYER = 'OUTPUT_LAYER'
4747

48+
def getIcon(self):
49+
return QIcon(os.path.join(pluginPath, 'images', 'ftools', 'centroids.png'))
50+
4851
def defineCharacteristics(self):
4952
self.name, self.i18n_name = self.trAlgorithm('Point on surface')
5053
self.group, self.i18n_group = self.trAlgorithm('Vector geometry tools')

‎python/plugins/processing/algs/qgis/PoleOfInaccessibility.py

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
# -*- coding: utf-8 -*-
2+
3+
"""
4+
***************************************************************************
5+
PoleOfInaccessibility.py
6+
------------------------
7+
Date : November 2016
8+
Copyright : (C) 2016 by Nyall Dawson
9+
Email : nyall dot dawson at gmail dot com
10+
***************************************************************************
11+
* *
12+
* This program is free software; you can redistribute it and/or modify *
13+
* it under the terms of the GNU General Public License as published by *
14+
* the Free Software Foundation; either version 2 of the License, or *
15+
* (at your option) any later version. *
16+
* *
17+
***************************************************************************
18+
"""
19+
20+
__author__ = 'Nyall Dawson'
21+
__date__ = 'November 2016'
22+
__copyright__ = '(C) 2016, Nyall Dawson'
23+
24+
# This will get replaced with a git SHA1 when you do a git archive323
25+
26+
__revision__ = '$Format:%H$'
27+
28+
import os
29+
30+
from qgis.core import QgsGeometry, QgsWkbTypes
31+
32+
from qgis.PyQt.QtGui import QIcon
33+
34+
from processing.core.GeoAlgorithm import GeoAlgorithm
35+
from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException
36+
from processing.core.parameters import ParameterVector, ParameterNumber
37+
from processing.core.outputs import OutputVector
38+
from processing.tools import dataobjects, vector
39+
40+
pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
41+
42+
43+
class PoleOfInaccessibility(GeoAlgorithm):
44+
45+
INPUT_LAYER = 'INPUT_LAYER'
46+
TOLERANCE = 'TOLERANCE'
47+
OUTPUT_LAYER = 'OUTPUT_LAYER'
48+
49+
def getIcon(self):
50+
return QIcon(os.path.join(pluginPath, 'images', 'ftools', 'centroids.png'))
51+
52+
def defineCharacteristics(self):
53+
self.name, self.i18n_name = self.trAlgorithm('Pole of Inaccessibility')
54+
self.group, self.i18n_group = self.trAlgorithm('Vector geometry tools')
55+
56+
self.addParameter(ParameterVector(self.INPUT_LAYER,
57+
self.tr('Input layer'),
58+
[dataobjects.TYPE_VECTOR_POLYGON]))
59+
self.addParameter(ParameterNumber(self.TOLERANCE,
60+
self.tr('Tolerance (layer units)'), default=1.0, minValue=0.0))
61+
self.addOutput(OutputVector(self.OUTPUT_LAYER, self.tr('Point'), datatype=[dataobjects.TYPE_VECTOR_POINT]))
62+
63+
def processAlgorithm(self, progress):
64+
layer = dataobjects.getObjectFromUri(
65+
self.getParameterValue(self.INPUT_LAYER))
66+
tolerance = self.getParameterValue(self.TOLERANCE)
67+
68+
writer = self.getOutputFromName(
69+
self.OUTPUT_LAYER).getVectorWriter(
70+
layer.fields(),
71+
QgsWkbTypes.Point,
72+
layer.crs())
73+
74+
features = vector.features(layer)
75+
total = 100.0 / len(features)
76+
77+
for current, input_feature in enumerate(features):
78+
output_feature = input_feature
79+
input_geometry = input_feature.geometry()
80+
if input_geometry:
81+
output_geometry = input_geometry.poleOfInaccessibility(tolerance)
82+
if not output_geometry:
83+
raise GeoAlgorithmExecutionException(
84+
self.tr('Error calculating pole of inaccessibility'))
85+
86+
output_feature.setGeometry(output_geometry)
87+
88+
writer.addFeature(output_feature)
89+
progress.setPercentage(int(current * total))
90+
91+
del writer

‎python/plugins/processing/algs/qgis/QGISAlgorithmProvider.py

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -176,6 +176,7 @@
176176
from .ExtractSpecificNodes import ExtractSpecificNodes
177177
from .GeometryByExpression import GeometryByExpression
178178
from .SnapGeometries import SnapGeometriesToLayer
179+
from .PoleOfInaccessibility import PoleOfInaccessibility
179180

180181
pluginPath = os.path.normpath(os.path.join(
181182
os.path.split(os.path.dirname(__file__))[0], os.pardir))
@@ -238,7 +239,8 @@ def __init__(self):
238239
IdwInterpolationZValue(), IdwInterpolationAttribute(),
239240
TinInterpolationZValue(), TinInterpolationAttribute(),
240241
RemoveNullGeometry(), ExtractByExpression(), ExtendLines(),
241-
ExtractSpecificNodes(), GeometryByExpression(), SnapGeometriesToLayer()
242+
ExtractSpecificNodes(), GeometryByExpression(), SnapGeometriesToLayer(),
243+
PoleOfInaccessibility()
242244
]
243245

244246
if hasMatplotlib:

‎python/plugins/processing/tests/testdata/expected/pole_inaccessibility_polys.gfs

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
<GMLFeatureClassList>
2+
<GMLFeatureClass>
3+
<Name>pole_inaccessibility_polys</Name>
4+
<ElementPath>pole_inaccessibility_polys</ElementPath>
5+
<!--POINT-->
6+
<GeometryType>1</GeometryType>
7+
<SRSName>EPSG:4326</SRSName>
8+
<DatasetSpecificInfo>
9+
<FeatureCount>6</FeatureCount>
10+
<ExtentXMin>0.50000</ExtentXMin>
11+
<ExtentXMax>6.58578</ExtentXMax>
12+
<ExtentYMin>-2.41422</ExtentYMin>
13+
<ExtentYMax>5.50000</ExtentYMax>
14+
</DatasetSpecificInfo>
15+
<PropertyDefn>
16+
<Name>name</Name>
17+
<ElementPath>name</ElementPath>
18+
<Type>String</Type>
19+
<Width>5</Width>
20+
</PropertyDefn>
21+
<PropertyDefn>
22+
<Name>intval</Name>
23+
<ElementPath>intval</ElementPath>
24+
<Type>Integer</Type>
25+
</PropertyDefn>
26+
<PropertyDefn>
27+
<Name>floatval</Name>
28+
<ElementPath>floatval</ElementPath>
29+
<Type>Real</Type>
30+
</PropertyDefn>
31+
</GMLFeatureClass>
32+
</GMLFeatureClassList>

‎python/plugins/processing/tests/testdata/expected/pole_inaccessibility_polys.gml

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
<?xml version="1.0" encoding="utf-8" ?>
2+
<ogr:FeatureCollection
3+
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
4+
xsi:schemaLocation=""
5+
xmlns:ogr="http://ogr.maptools.org/"
6+
xmlns:gml="http://www.opengis.net/gml">
7+
<gml:boundedBy>
8+
<gml:Box>
9+
<gml:coord><gml:X>0.5</gml:X><gml:Y>-2.414215087890625</gml:Y></gml:coord>
10+
<gml:coord><gml:X>6.585784912109375</gml:X><gml:Y>5.5</gml:Y></gml:coord>
11+
</gml:Box>
12+
</gml:boundedBy>
13+
14+
<gml:featureMember>
15+
<ogr:pole_inaccessibility_polys fid="polys.0">
16+
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>0.5,0.5</gml:coordinates></gml:Point></ogr:geometryProperty>
17+
<ogr:name>aaaaa</ogr:name>
18+
<ogr:intval>33</ogr:intval>
19+
<ogr:floatval>44.123456</ogr:floatval>
20+
</ogr:pole_inaccessibility_polys>
21+
</gml:featureMember>
22+
<gml:featureMember>
23+
<ogr:pole_inaccessibility_polys fid="polys.1">
24+
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>4.999996185302734,4.414211273193359</gml:coordinates></gml:Point></ogr:geometryProperty>
25+
<ogr:name>Aaaaa</ogr:name>
26+
<ogr:intval>-33</ogr:intval>
27+
<ogr:floatval>0</ogr:floatval>
28+
</ogr:pole_inaccessibility_polys>
29+
</gml:featureMember>
30+
<gml:featureMember>
31+
<ogr:pole_inaccessibility_polys fid="polys.2">
32+
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>2.5,5.5</gml:coordinates></gml:Point></ogr:geometryProperty>
33+
<ogr:name>bbaaa</ogr:name>
34+
<ogr:floatval>0.123</ogr:floatval>
35+
</ogr:pole_inaccessibility_polys>
36+
</gml:featureMember>
37+
<gml:featureMember>
38+
<ogr:pole_inaccessibility_polys fid="polys.3">
39+
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>6.585784912109375,-2.414215087890625</gml:coordinates></gml:Point></ogr:geometryProperty>
40+
<ogr:name>ASDF</ogr:name>
41+
<ogr:intval>0</ogr:intval>
42+
</ogr:pole_inaccessibility_polys>
43+
</gml:featureMember>
44+
<gml:featureMember>
45+
<ogr:pole_inaccessibility_polys fid="polys.4">
46+
<ogr:intval>120</ogr:intval>
47+
<ogr:floatval>-100291.43213</ogr:floatval>
48+
</ogr:pole_inaccessibility_polys>
49+
</gml:featureMember>
50+
<gml:featureMember>
51+
<ogr:pole_inaccessibility_polys fid="polys.5">
52+
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>4.289703369140625,-0.232696533203125</gml:coordinates></gml:Point></ogr:geometryProperty>
53+
<ogr:name>elim</ogr:name>
54+
<ogr:intval>2</ogr:intval>
55+
<ogr:floatval>3.33</ogr:floatval>
56+
</ogr:pole_inaccessibility_polys>
57+
</gml:featureMember>
58+
</ogr:FeatureCollection>

‎python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1464,3 +1464,15 @@ tests:
14641464
OUTPUT:
14651465
name: expected/snap_point_to_lines_prefer_closest.gml
14661466
type: vector
1467+
1468+
- algorithm: qgis:poleofinaccessibility
1469+
name: Pole of inaccessibility (polygons)
1470+
params:
1471+
INPUT_LAYER:
1472+
name: polys.gml
1473+
type: vector
1474+
TOLERANCE: 1.0e-05
1475+
results:
1476+
OUTPUT_LAYER:
1477+
name: expected/pole_inaccessibility_polys.gml
1478+
type: vector

‎src/core/geometry/qgsgeometry.cpp

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1480,6 +1480,13 @@ QgsGeometry QgsGeometry::pointOnSurface() const
14801480
return QgsGeometry( pt.clone() );
14811481
}
14821482

1483+
QgsGeometry QgsGeometry::poleOfInaccessibility( double precision ) const
1484+
{
1485+
QgsInternalGeometryEngine engine( *this );
1486+
1487+
return engine.poleOfInaccessibility( precision );
1488+
}
1489+
14831490
QgsGeometry QgsGeometry::convexHull() const
14841491
{
14851492
if ( !d->geometry )

‎src/core/geometry/qgsgeometry.h

Lines changed: 24 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -563,15 +563,37 @@ class CORE_EXPORT QgsGeometry
563563
//! Returns a simplified version of this geometry using a specified tolerance value
564564
QgsGeometry simplify( double tolerance ) const;
565565

566-
/** Returns the center of mass of a geometry
566+
/**
567+
* Returns the center of mass of a geometry.
567568
* @note for line based geometries, the center point of the line is returned,
568569
* and for point based geometries, the point itself is returned
570+
* @see pointOnSurface()
571+
* @see poleOfInaccessibility()
569572
*/
570573
QgsGeometry centroid() const;
571574

572-
//! Returns a point within a geometry
575+
/**
576+
* Returns a point guaranteed to lie on the surface of a geometry. While the centroid()
577+
* of a geometry may be located outside of the geometry itself (eg for concave shapes),
578+
* the point on surface will always be inside the geometry.
579+
* @see centroid()
580+
* @see poleOfInaccessibility()
581+
*/
573582
QgsGeometry pointOnSurface() const;
574583

584+
/**
585+
* Calculates the approximate pole of inaccessibility for a surface, which is the
586+
* most distant internal point from the boundary of the surface. This function
587+
* uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative
588+
* approach guaranteed to find the true pole of inaccessibility within a specified
589+
* tolerance. More precise tolerances require more iterations and will take longer
590+
* to calculate.
591+
* @see centroid()
592+
* @see pointOnSurface()
593+
* @note added in QGIS 3.0
594+
*/
595+
QgsGeometry poleOfInaccessibility( double precision ) const;
596+
575597
//! Returns the smallest convex polygon that contains all the points in the geometry.
576598
QgsGeometry convexHull() const;
577599

‎src/core/geometry/qgsinternalgeometryengine.cpp

Lines changed: 148 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,8 +21,11 @@
2121
#include "qgspolygon.h"
2222
#include "qgsmulticurve.h"
2323
#include "qgsgeometry.h"
24+
#include "qgsgeometryutils.h"
25+
2426

2527
#include <QTransform>
28+
#include <queue>
2629

2730
QgsInternalGeometryEngine::QgsInternalGeometryEngine( const QgsGeometry& geometry )
2831
: mGeometry( geometry.geometry() )
@@ -36,7 +39,7 @@ QgsInternalGeometryEngine::QgsInternalGeometryEngine( const QgsGeometry& geometr
3639
* See details in QEP #17
3740
****************************************************************************/
3841

39-
QgsGeometry QgsInternalGeometryEngine::extrude( double x, double y )
42+
QgsGeometry QgsInternalGeometryEngine::extrude( double x, double y ) const
4043
{
4144
QList<QgsLineString*> linesToProcess;
4245

@@ -87,3 +90,147 @@ QgsGeometry QgsInternalGeometryEngine::extrude( double x, double y )
8790

8891
return QgsGeometry();
8992
}
93+
94+
95+
96+
// polylabel implementation
97+
// ported from the original Javascript implementation developed by Vladimir Agafonkin
98+
// originally licensed under the ISC License
99+
100+
class Cell
101+
{
102+
public:
103+
Cell( double x, double y, double h, const QgsPolygonV2* polygon )
104+
: x( x )
105+
, y( y )
106+
, h( h )
107+
, d( polygon->pointDistanceToBoundary( x, y ) )
108+
, max( d + h * M_SQRT2 )
109+
{}
110+
111+
//! Cell center x
112+
double x;
113+
//! Cell center y
114+
double y;
115+
//! Half the cell size
116+
double h;
117+
//! Distance from cell center to polygon
118+
double d;
119+
//! Maximum distance to polygon within a cell
120+
double max;
121+
};
122+
123+
struct GreaterThanByMax
124+
{
125+
bool operator()( const Cell* lhs, const Cell* rhs )
126+
{
127+
return rhs->max > lhs->max;
128+
}
129+
};
130+
131+
Cell* getCentroidCell( const QgsPolygonV2* polygon )
132+
{
133+
double area = 0;
134+
double x = 0;
135+
double y = 0;
136+
137+
const QgsLineString* exterior = static_cast< const QgsLineString*>( polygon->exteriorRing() );
138+
int len = exterior->numPoints() - 1; //assume closed
139+
for ( int i = 0, j = len - 1; i < len; j = i++ )
140+
{
141+
double aX = exterior->xAt( i );
142+
double aY = exterior->yAt( i );
143+
double bX = exterior->xAt( j );
144+
double bY = exterior->yAt( j );
145+
double f = aX * bY - bX * aY;
146+
x += ( aX + bX ) * f;
147+
y += ( aY + bY ) * f;
148+
area += f * 3;
149+
}
150+
if ( area == 0 )
151+
return new Cell( exterior->xAt( 0 ), exterior->yAt( 0 ), 0, polygon );
152+
else
153+
return new Cell( x / area, y / area, 0.0, polygon );
154+
}
155+
156+
157+
QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision ) const
158+
{
159+
if ( !mGeometry || mGeometry->isEmpty() )
160+
return QgsGeometry();
161+
162+
if ( precision <= 0 )
163+
return QgsGeometry();
164+
165+
const QgsSurface* surface = dynamic_cast< const QgsSurface* >( mGeometry );
166+
if ( !surface )
167+
return QgsGeometry();
168+
169+
QScopedPointer< QgsPolygonV2 > segmentizedPoly;
170+
const QgsPolygonV2* polygon = dynamic_cast< const QgsPolygonV2* >( mGeometry );
171+
if ( !polygon )
172+
{
173+
segmentizedPoly.reset( static_cast< QgsPolygonV2*>( surface->segmentize() ) );
174+
polygon = segmentizedPoly.data();
175+
}
176+
177+
// start with the bounding box
178+
QgsRectangle bounds = polygon->boundingBox();
179+
180+
// initial parameters
181+
double cellSize = qMin( bounds.width(), bounds.height() );
182+
183+
if ( qgsDoubleNear( cellSize, 0.0 ) )
184+
return QgsGeometry( new QgsPointV2( bounds.xMinimum(), bounds.yMinimum() ) );
185+
186+
double h = cellSize / 2.0;
187+
std::priority_queue< Cell*, std::vector<Cell*>, GreaterThanByMax > cellQueue;
188+
189+
// cover polygon with initial cells
190+
for ( double x = bounds.xMinimum(); x < bounds.xMaximum(); x += cellSize )
191+
{
192+
for ( double y = bounds.yMinimum(); y < bounds.yMaximum(); y += cellSize )
193+
{
194+
cellQueue.push( new Cell( x + h, y + h, h, polygon ) );
195+
}
196+
}
197+
198+
// take centroid as the first best guess
199+
QScopedPointer< Cell > bestCell( getCentroidCell( polygon ) );
200+
201+
// special case for rectangular polygons
202+
QScopedPointer< Cell > bboxCell( new Cell( bounds.xMinimum() + bounds.width() / 2.0,
203+
bounds.yMinimum() + bounds.height() / 2.0,
204+
0, polygon ) );
205+
if ( bboxCell->d > bestCell->d )
206+
{
207+
bestCell.reset( bboxCell.take() );
208+
}
209+
210+
while ( cellQueue.size() > 0 )
211+
{
212+
// pick the most promising cell from the queue
213+
QScopedPointer< Cell > cell( cellQueue.top() );
214+
cellQueue.pop();
215+
Cell* currentCell = cell.data();
216+
217+
// update the best cell if we found a better one
218+
if ( currentCell->d > bestCell->d )
219+
{
220+
bestCell.reset( cell.take() );
221+
}
222+
223+
// do not drill down further if there's no chance of a better solution
224+
if ( currentCell->max - bestCell->d <= precision )
225+
continue;
226+
227+
// split the cell into four cells
228+
h = currentCell->h / 2.0;
229+
cellQueue.push( new Cell( currentCell->x - h, currentCell->y - h, h, polygon ) );
230+
cellQueue.push( new Cell( currentCell->x + h, currentCell->y - h, h, polygon ) );
231+
cellQueue.push( new Cell( currentCell->x - h, currentCell->y + h, h, polygon ) );
232+
cellQueue.push( new Cell( currentCell->x + h, currentCell->y + h, h, polygon ) );
233+
}
234+
235+
return QgsGeometry( new QgsPointV2( bestCell->x, bestCell->y ) );
236+
}

‎src/core/geometry/qgsinternalgeometryengine.h

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,17 @@ class QgsInternalGeometryEngine
4848
* @param y offset in y direction
4949
* @return an extruded polygon
5050
*/
51-
QgsGeometry extrude( double x, double y );
51+
QgsGeometry extrude( double x, double y ) const;
52+
53+
/**
54+
* Calculates the approximate pole of inaccessibility for a surface, which is the
55+
* most distant internal point from the boundary of the surface. This function
56+
* uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative
57+
* approach guaranteed to find the true pole of inaccessibility within a specified
58+
* tolerance. More precise tolerances require more iterations and will take longer
59+
* to calculate.
60+
*/
61+
QgsGeometry poleOfInaccessibility( double precision ) const;
5262

5363
private:
5464
const QgsAbstractGeometry* mGeometry;

‎tests/src/core/testqgsgeometry.cpp

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,7 @@ class TestQgsGeometry : public QObject
109109
void wkbInOut();
110110

111111
void segmentizeCircularString();
112+
void poleOfInaccessibility();
112113

113114
private:
114115
//! A helper method to do a render check to see if the geometry op is as expected
@@ -3880,5 +3881,126 @@ void TestQgsGeometry::segmentizeCircularString()
38803881
QVERIFY( points.contains( QgsPointV2( 0.5, 0.5 ) ) );
38813882
}
38823883

3884+
void TestQgsGeometry::poleOfInaccessibility()
3885+
{
3886+
QString poly1Wkt( "Polygon ((3116 3071, 3394 3431, 3563 3362, 3611 3205, 3599 3181, 3477 3281, 3449 3160, 3570 3127, 3354 3116,"
3887+
" 3436 3008, 3158 2907, 2831 2438, 3269 2916, 3438 2885, 3295 2799, 3407 2772, 3278 2629, 3411 2689, 3329 2611,"
3888+
" 3360 2531, 3603 2800, 3598 2501, 3317 2429, 3329 2401, 3170 2340, 3142 2291, 3524 2403, 3598 2233, 3460 2117,"
3889+
" 3590 1931, 3364 1753, 3597 1875, 3639 1835, 3660 1733, 3600 1771, 3538 1694, 3661 1732, 3359 1554, 3334 1554,"
3890+
" 3341 1588, 3317 1588, 3305 1644, 3286 1656, 3115 1255, 3072 1252, 3078 1335, 3046 1355, 2984 1234, 2983 1409,"
3891+
" 2876 1222, 2525 1161, 2488 787, 2162 913, 2079 661, 2270 380, 2188 823, 2510 592, 2659 992, 2911 1118, 2943 938,"
3892+
" 2957 1097, 3092 1002, 3006 1097, 3233 1282, 3325 1291, 3296 1116, 3333 1115, 3391 1333, 3434 1274, 3413 1326,"
3893+
" 3449 1327, 3473 1408, 3490 1430, 3526 1434, 3198 -112, 3263 -87, 3289 -128, 4224 -128, 4224 -128, 3474 -128,"
3894+
" 3475 -116, 3486 -120, 3458 -49, 3523 -78, 3513 -128, 3509 -119, 3509 -128, 3717 -128, 3705 -60, 3735 -16,"
3895+
" 3714 38, 3758 88, 3825 47, 3812 -11, 3859 -51, 3871 49, 3760 149, 3636 -74, 3510 126, 3501 245, 3504 270,"
3896+
" 3511 284, 3582 16, 3631 19, 3569 125, 3570 193, 3610 212, 3583 119, 3655 29, 3738 170, 3561 466, 3826 549,"
3897+
" 3527 604, 3609 833, 3681 798, 3956 1127, 3917 964, 4043 850, 4049 1096, 4193 1052, 4191 1078, 4208 1106,"
3898+
" 4222 1110, 4224 1109, 4224 1144, 4202 1158, 4177 1161, 4182 1181, 4075 1201, 4141 1275, 4215 1215, 4221 1223,"
3899+
" 4219 1231, 4224 1243, 4224 1257, 4224 1262, 4224 1345, 4224 1339, 4224 1328, 4215 1335, 4203 1355, 4215 1369,"
3900+
" 4224 1363, 4215 1377, 4208 1387, 4217 1401, 4224 1403, 4224 1520, 4219 1535, 4221 1544, 4199 1593, 4223 1595,"
3901+
" 4206 1626, 4214 1648, 4224 1645, 4224 1640, 4224 2108, 4220 2125, 4224 2125, 4224 2143, 4205 2141, 4180 2159,"
3902+
" 4201 2182, 4163 2189, 4176 2229, 4199 2211, 4210 2218, 4212 2210, 4223 2214, 4224 2207, 4224 2216, 4217 2225,"
3903+
" 4221 2233, 4203 2227, 4209 2248, 4185 2240, 4198 2276, 4144 2218, 4091 2343, 4119 2332, 4121 2347, 4155 2337,"
3904+
" 4180 2355, 4200 2342, 4201 2354, 4213 2352, 4224 2348, 4224 2356, 4207 2361, 4184 2358, 4176 2367, 4106 2355,"
3905+
" 3983 2765, 4050 3151, 4139 2720, 4209 2589, 4211 2600, 4219 2599, 4224 2592, 4224 2574, 4223 2566, 4224 2562,"
3906+
" 4224 2553, 4224 2552, 4224 -128, 4224 4224, 4205 4224, 4015 3513, 3993 3494, 3873 3533, 3887 3539, 3923 3524,"
3907+
" 3950 3529, 4018 3572, 3987 3633, 3983 3571, 3955 3583, 3936 3547, 3882 3539, 3913 3557, 3920 3598, 3901 3596,"
3908+
" 3923 3631, 3914 3628, 3919 3647, 3922 3656, 3917 3666, 3523 3438, 3631 3564, 3527 3597, 3718 3655, 3578 3672,"
3909+
" 3660 3867, 3543 3628, 3416 3725, 3487 3503, 3274 3583, 3271 3644, 3197 3671, 3210 3775, 3184 3788, 3181 3672,"
3910+
" 3306 3521, 3292 3508, 3229 3565, 3219 3564, 3216 3574, 3192 3578, 3297 3444, 3089 3395, 3029 3028, 2973 3133,"
3911+
" 2529 2945, 2538 2811, 2461 2901, 2170 2839, 2121 2797, 2156 2733, 2105 2709, 2096 2695, 2114 2621, 2102 2693,"
3912+
" 2168 2738, 2167 2778, 2447 2765, 2441 2866, 2527 2793, 2670 2938, 2626 2651, 2688 2623, 2740 2922, 3084 2960,"
3913+
" 3116 3071),(4016 1878, 4029 1859, 4008 1839, 4006 1863, 4016 1878),(3315 1339, 3331 1324, 3290 1293, 3305 1315,"
3914+
" 3315 1339),(4136 3071, 4136 3080, 4143 3020, 4137 3036, 4136 3071),(4218 3073, 4183 3114, 4117 3157, 4159 3147,"
3915+
" 4218 3073),(3912 3542, 3934 3536, 3955 3536, 3937 3527, 3900 3540, 3912 3542),(4050 3172, 4043 3210, 4085 3209,"
3916+
" 4090 3179, 4050 3172),(4151 2998, 4158 2977, 4159 2946, 4147 2988, 4151 2998),(2920 3005, 2935 2994, 2864 2973,"
3917+
" 2905 3016, 2920 3005),(3571 2424, 3545 2469, 3608 2480, 3596 2434, 3571 2424),(4095 1229, 4073 1234, 4076 1293,"
3918+
" 4121 1285, 4095 1229),(4173 1536, 4153 1576, 4166 1585, 4198 1571, 4188 1532, 4213 1535, 4224 1512, 4224 1511,"
3919+
" 4209 1511, 4198 1506, 4190 1517, 4194 1509, 4192 1499, 4200 1496, 4202 1504, 4224 1510, 4224 1488, 4215 1486,"
3920+
" 4216 1478, 4224 1472, 4224 1464, 4207 1458, 4193 1464, 4173 1536),(3934 1537, 3968 1630, 3960 1666, 3968 1673,"
3921+
" 3975 1562, 3934 1537),(4182 1653, 4196 1624, 4166 1614, 4157 1674, 4216 1671, 4182 1653),(4200 1619, 4196 1620,"
3922+
" 4200 1632, 4195 1642, 4207 1648, 4200 1619),(4026 1835, 4025 1830, 4016 1808, 4007 1836, 4026 1835),(4199 1384,"
3923+
" 4182 1389, 4206 1412, 4216 1401, 4199 1384),(3926 1251, 3969 1206, 3913 1149, 3878 1173, 3876 1229, 3926 1251),"
3924+
" (3926 1354, 3958 1389, 3997 1384, 3991 1352, 3960 1322, 3955 1299, 3926 1354),(3964 1319, 3974 1329, 3984 1285,"
3925+
" 3963 1301, 3964 1319),(3687 959, 3696 903, 3678 885, 3665 930, 3687 959),(3452 79, 3437 124, 3456 149, 3476 141,"
3926+
" 3452 79),(3751 927, 3738 906, 3719 942, 3739 929, 3751 927))" );
3927+
3928+
QString poly2Wkt( "Polygon ((-128 4224, -128 2734, -11 2643, -101 2919, 120 2654, 19 2846, 217 2897, -13 2873, -79 2998, -106 2989,"
3929+
" -128 3002, -128 4224, -128 4224, 1249 4224, 1247 4199, 1231 4132, 909 3902, 775 3278, 509 2903, 470 2603, 663 2311,"
3930+
" 889 2134, 989 2146, 534 2585, 575 2880, 833 3112, 833 3037, 1025 2977, 834 3096, 946 3217, 923 3153, 971 3094,"
3931+
" 997 3103, 1166 3423, 1198 3329, 1233 3367, 1270 3327, 1274 3354, 1290 3360, 1321 3322, 1331 3318, 1343 3325,"
3932+
" 1310 3073, 1363 3093, 1413 3186, 1325 3398, 1458 3364, 1493 3426, 1526 3394, 1588 3465, 1417 3824, 1458 3825,"
3933+
" 1522 3849, 1729 3488, 2018 3223, 1908 3291, 1924 3238, 1899 3187, 1913 3150, 2070 3041, 2102 3063, 2112 3053,"
3934+
" 2168 3085, 2101 2994, 2265 2863, 1890 2593, 2106 2713, 2130 2706, 2108 2678, 2117 2647, 2105 2636, 2099 2604,"
3935+
" 2094 2591, 2088 2575, 2088 2564, 2249 2751, 2441 2618, 2689 1728, 2681 2323, 2776 1685, 2711 1708, 2673 1680,"
3936+
" 2675 1639, 2810 1522, 2765 1535, 2734 1510, 2850 1332, 2863 1186, 2847 1025, 2832 985, 2739 1025, 2850 1090,"
3937+
" 2859 1174, 2853 1235, 2827 1235, 2839 1279, 2811 1286, 2751 1272, 2851 1160, 2788 1159, 2724 1204, 2713 1198,"
3938+
" 2708 1213, 2699 1221, 2724 1179, 2797 1145, 2785 1123, 2848 1143, 2821 1092, 2284 980, 2288 963, 2264 962,"
3939+
" 2261 948, 2247 958, 2194 947, 2120 900, 2047 809, 2434 923, 2450 817, 2528 737, 2527 667, 2641 544, 2641 460,"
3940+
" 2710 452, 2687 447, 2681 435, 2693 330, 2772 327, 2766 291, 2779 245, 2769 219, 2771 198, 2816 219, 2781 342,"
3941+
" 2612 647, 2903 188, 2803 539, 2950 206, 2927 342, 3123 339, 3092 300, 3073 175, 3082 160, 3412 -103, 4224 -128,"
3942+
" 4032 167, 3572 63, 3554 109, 3606 211, 3604 232, 3454 124, 3332 345, 3237 212, 3181 415, 2953 364, 2761 692,"
3943+
" 2819 788, 2997 769, 2997 626, 3199 659, 3155 730, 2929 805, 2908 841, 3218 717, 3276 744, 3246 723, 3270 658,"
3944+
" 4223 473, 4224 589, 3906 681, 3818 677, 3915 686, 4224 592, 4224 4224, -128 4224, -128 4224),(2049 3181,"
3945+
" 2224 3084, 2204 3040, 2169 3036, 2174 3084, 2155 3102, 2126 3123, 2109 3127, 2103 3102, 2049 3181),(1578 3811,"
3946+
" 1567 3883, 1675 3768, 1658 3712, 1659 3751, 1578 3811),(2304 2930, 2335 2918, 2287 2880, 2279 2926, 2304 2930),"
3947+
" (2316 2895, 2317 2895, 2316 2901, 2331 2901, 2332 2889, 2316 2895),(2304 2850, 2335 2861, 2344 2910, 2357 2828,"
3948+
" 2304 2850),(1714 3583, 1725 3638, 1682 3717, 1797 3564, 1714 3583),(1537 3873, 1456 3827, 1405 3876, 1461 3898,"
3949+
" 1537 3873),(2547 2560, 2375 2815, 2384 2825, 2470 2722, 2685 2336, 2648 2310, 2547 2560),(2107 3073, 2107 3098,"
3950+
" 2144 3086, 2122 3073, 2113 3059, 2107 3073),(2117 3120, 2156 3099, 2164 3083, 2106 3103, 2117 3120),(2303 2981,"
3951+
" 2226 3010, 2232 3064, 2286 3012, 2344 2928, 2303 2981),(2304 2858, 2291 2864, 2303 2885, 2324 2870, 2304 2858),"
3952+
" (2175 3002, 2179 2983, 2152 2991, 2175 3002, 2175 3002),(2175 3022, 2150 3017, 2144 3048, 2159 3031, 2184 3030,"
3953+
" 2175 3022),(2265 2953, 2270 2950, 2262 2950, 2254 2963, 2253 2979, 2265 2953),(2229 2928, 2250 2928, 2241 2922,"
3954+
" 2239 2914, 2224 2900, 2237 2914, 2229 2928),(2531 2473, 2520 2466, 2521 2474, 2511 2503, 2531 2473),(2547 2526,"
3955+
" 2529 2519, 2528 2541, 2544 2538, 2547 2526),(2559 646, 2513 810, 2463 835, 2462 930, 2746 731, 2559 646),(3840 653,"
3956+
" 3809 641, 3718 689, 3547 733, 3553 712, 3440 780, 3840 653),(3327 741, 3242 750, 3195 745, 3180 758, 3326 764,"
3957+
" 3440 742, 3327 741),(3282 702, 3265 699, 3273 721, 3290 714, 3282 702),(3762 662, 3783 654, 3786 638, 3742 653,"
3958+
" 3762 662),(3071 703, 3082 741, 3079 655, 3064 657, 3071 703),(3881 637, 3904 647, 3914 626, 3861 638, 3881 637))" );
3959+
3960+
QgsGeometry poly1 = QgsGeometry::fromWkt( poly1Wkt );
3961+
QgsGeometry poly2 = QgsGeometry::fromWkt( poly2Wkt );
3962+
3963+
QgsPoint point = poly1.poleOfInaccessibility( 1 ).asPoint();
3964+
QGSCOMPARENEAR( point.x(), 3867.37, 0.01 );
3965+
QGSCOMPARENEAR( point.y(), 2126.45, 0.01 );
3966+
3967+
point = poly1.poleOfInaccessibility( 50 ).asPoint();
3968+
QGSCOMPARENEAR( point.x(), 3855.33, 0.01 );
3969+
QGSCOMPARENEAR( point.y(), 2117.55, 0.01 );
3970+
3971+
point = poly2.poleOfInaccessibility( 1 ).asPoint();
3972+
QGSCOMPARENEAR( point.x(), 3263.50, 0.01 );
3973+
QGSCOMPARENEAR( point.y(), 3263.50, 0.01 );
3974+
3975+
//test degenerate polygons
3976+
QgsGeometry degen1 = QgsGeometry::fromWkt( "Polygon(( 0 0, 1 0, 2 0, 0 0 ))" );
3977+
point = degen1.poleOfInaccessibility( 1 ).asPoint();
3978+
QGSCOMPARENEAR( point.x(), 0, 0.01 );
3979+
QGSCOMPARENEAR( point.y(), 0, 0.01 );
3980+
3981+
QgsGeometry degen2 = QgsGeometry::fromWkt( "Polygon(( 0 0, 1 0, 1 1 , 1 0, 0 0 ))" );
3982+
point = degen2.poleOfInaccessibility( 1 ).asPoint();
3983+
QGSCOMPARENEAR( point.x(), 0, 0.01 );
3984+
QGSCOMPARENEAR( point.y(), 0, 0.01 );
3985+
3986+
//empty geometry
3987+
QVERIFY( QgsGeometry().poleOfInaccessibility( 1 ).isEmpty() );
3988+
3989+
// not a polygon
3990+
QgsGeometry lineString = QgsGeometry::fromWkt( "LineString(1 0, 2 2 )" );
3991+
QVERIFY( lineString.poleOfInaccessibility( 1 ).isEmpty() );
3992+
3993+
// invalid threshold
3994+
QVERIFY( poly1.poleOfInaccessibility( -1 ).isEmpty() );
3995+
QVERIFY( poly1.poleOfInaccessibility( 0 ).isEmpty() );
3996+
3997+
// curved geometry
3998+
QgsGeometry curved = QgsGeometry::fromWkt( "CurvePolygon( CompoundCurve( CircularString(-0.44 0.35, 0.51 0.34, 0.56 0.21, 0.11 -0.33, 0.15 -0.35,"
3999+
"-0.93 -0.30, -1.02 -0.22, -0.49 0.01, -0.23 -0.04),(-0.23 -0.04, -0.44, 0.35)))" );
4000+
point = curved.poleOfInaccessibility( 0.01 ).asPoint();
4001+
QGSCOMPARENEAR( point.x(), -0.4324, 0.0001 );
4002+
QGSCOMPARENEAR( point.y(), -0.2434, 0.0001 );
4003+
}
4004+
38834005
QTEST_MAIN( TestQgsGeometry )
38844006
#include "testqgsgeometry.moc"

0 commit comments

Comments
 (0)
Please sign in to comment.