Fix test failure, report distance from pole

qgis · Nov 14, 2016 · b440939 · b440939
1 parent d6f09c0
commit b440939
Show file tree

Hide file tree

Showing 9 changed files with 41 additions and 10 deletions.
diff --git a/python/core/geometry/qgsgeometry.sip b/python/core/geometry/qgsgeometry.sip
@@ -535,11 +535,12 @@ class QgsGeometry
      * 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.
+     * Optionally, the distance to the polygon boundary from the pole can be stored.
      * @see centroid()
      * @see pointOnSurface()
      * @note added in QGIS 3.0
      */
-    QgsGeometry poleOfInaccessibility( double precision ) const;
+    QgsGeometry poleOfInaccessibility( double precision, double* distanceToBoundary /Out/ = nullptr ) const;
 
     /** Returns the smallest convex polygon that contains all the points in the geometry. */
     QgsGeometry convexHull() const;

diff --git a/python/plugins/processing/algs/help/qgis.yaml b/python/plugins/processing/algs/help/qgis.yaml
@@ -340,6 +340,8 @@ qgis:polarplot: >
 qgis:poleofinaccessibility: >
   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.
 
+  The distance from the calculated pole to the polygon boundary will be stored as a new attribute in the output layer.
+
 qgis:polygoncentroids: >
   This algorithm creates a new point layer, with points representing the centroid of polygons of an input layer.
 

diff --git a/python/plugins/processing/algs/qgis/PoleOfInaccessibility.py b/python/plugins/processing/algs/qgis/PoleOfInaccessibility.py
@@ -27,8 +27,9 @@
 
 import os
 
-from qgis.core import QgsGeometry, QgsWkbTypes
+from qgis.core import QgsGeometry, QgsWkbTypes, QgsField, NULL
 
+from qgis.PyQt.QtCore import QVariant
 from qgis.PyQt.QtGui import QIcon
 
 from processing.core.GeoAlgorithm import GeoAlgorithm
@@ -65,9 +66,12 @@ def processAlgorithm(self, progress):
             self.getParameterValue(self.INPUT_LAYER))
         tolerance = self.getParameterValue(self.TOLERANCE)
 
+        fields = layer.fields()
+        fields.append(QgsField('dist_pole', QVariant.Double))
+
         writer = self.getOutputFromName(
             self.OUTPUT_LAYER).getVectorWriter(
-                layer.fields(),
+                fields,
                 QgsWkbTypes.Point,
                 layer.crs())
 
@@ -78,12 +82,19 @@ def processAlgorithm(self, progress):
             output_feature = input_feature
             input_geometry = input_feature.geometry()
             if input_geometry:
-                output_geometry = input_geometry.poleOfInaccessibility(tolerance)
+                output_geometry, distance = input_geometry.poleOfInaccessibility(tolerance)
                 if not output_geometry:
                     raise GeoAlgorithmExecutionException(
                         self.tr('Error calculating pole of inaccessibility'))
+                attrs = input_feature.attributes()
+                attrs.append(distance)
+                output_feature.setAttributes(attrs)
 
                 output_feature.setGeometry(output_geometry)
+            else:
+                attrs = input_feature.attributes()
+                attrs.append(NULL)
+                output_feature.setAttributes(attrs)
 
             writer.addFeature(output_feature)
             progress.setPercentage(int(current * total))

diff --git a/python/plugins/processing/tests/testdata/expected/pole_inaccessibility_polys.gml b/python/plugins/processing/tests/testdata/expected/pole_inaccessibility_polys.gml
@@ -17,6 +17,7 @@
       <ogr:name>aaaaa</ogr:name>
       <ogr:intval>33</ogr:intval>
       <ogr:floatval>44.123456</ogr:floatval>
+      <ogr:dist_pole>1.5</ogr:dist_pole>
     </ogr:pole_inaccessibility_polys>
   </gml:featureMember>
   <gml:featureMember>
@@ -25,20 +26,23 @@
       <ogr:name>Aaaaa</ogr:name>
       <ogr:intval>-33</ogr:intval>
       <ogr:floatval>0</ogr:floatval>
+      <ogr:dist_pole>0.414211273193359</ogr:dist_pole>
     </ogr:pole_inaccessibility_polys>
   </gml:featureMember>
   <gml:featureMember>
     <ogr:pole_inaccessibility_polys fid="polys.2">
       <ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>2.5,5.5</gml:coordinates></gml:Point></ogr:geometryProperty>
       <ogr:name>bbaaa</ogr:name>
       <ogr:floatval>0.123</ogr:floatval>
+      <ogr:dist_pole>0.5</ogr:dist_pole>
     </ogr:pole_inaccessibility_polys>
   </gml:featureMember>
   <gml:featureMember>
     <ogr:pole_inaccessibility_polys fid="polys.3">
       <ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>6.585784912109375,-2.414215087890625</gml:coordinates></gml:Point></ogr:geometryProperty>
       <ogr:name>ASDF</ogr:name>
       <ogr:intval>0</ogr:intval>
+      <ogr:dist_pole>0.585784912109375</ogr:dist_pole>
     </ogr:pole_inaccessibility_polys>
   </gml:featureMember>
   <gml:featureMember>
@@ -53,6 +57,7 @@
       <ogr:name>elim</ogr:name>
       <ogr:intval>2</ogr:intval>
       <ogr:floatval>3.33</ogr:floatval>
+      <ogr:dist_pole>1.71028189541736</ogr:dist_pole>
     </ogr:pole_inaccessibility_polys>
   </gml:featureMember>
 </ogr:FeatureCollection>
diff --git a/src/core/geometry/qgsgeometry.cpp b/src/core/geometry/qgsgeometry.cpp
@@ -1480,11 +1480,11 @@ QgsGeometry QgsGeometry::pointOnSurface() const
   return QgsGeometry( pt.clone() );
 }
 
-QgsGeometry QgsGeometry::poleOfInaccessibility( double precision ) const
+QgsGeometry QgsGeometry::poleOfInaccessibility( double precision, double* distanceToBoundary ) const
 {
   QgsInternalGeometryEngine engine( *this );
 
-  return engine.poleOfInaccessibility( precision );
+  return engine.poleOfInaccessibility( precision, distanceToBoundary );
 }
 
 QgsGeometry QgsGeometry::convexHull() const

diff --git a/src/core/geometry/qgsgeometry.h b/src/core/geometry/qgsgeometry.h
@@ -588,11 +588,12 @@ class CORE_EXPORT QgsGeometry
      * 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.
+     * Optionally, the distance to the polygon boundary from the pole can be stored.
      * @see centroid()
      * @see pointOnSurface()
      * @note added in QGIS 3.0
      */
-    QgsGeometry poleOfInaccessibility( double precision ) const;
+    QgsGeometry poleOfInaccessibility( double precision, double* distanceToBoundary = nullptr ) const;
 
     //! Returns the smallest convex polygon that contains all the points in the geometry.
     QgsGeometry convexHull() const;

diff --git a/src/core/geometry/qgsinternalgeometryengine.cpp b/src/core/geometry/qgsinternalgeometryengine.cpp
@@ -97,6 +97,7 @@ QgsGeometry QgsInternalGeometryEngine::extrude( double x, double y ) const
 // ported from the original Javascript implementation developed by Vladimir Agafonkin
 // originally licensed under the ISC License
 
+/// @cond PRIVATE
 class Cell
 {
   public:
@@ -153,9 +154,13 @@ Cell* getCentroidCell( const QgsPolygonV2* polygon )
     return new Cell( x / area, y / area, 0.0, polygon );
 }
 
+///@endcond
 
-QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision ) const
+QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision , double* distanceFromBoundary ) const
 {
+  if ( distanceFromBoundary )
+    *distanceFromBoundary = DBL_MAX;
+
   if ( !mGeometry || mGeometry->isEmpty() )
     return QgsGeometry();
 
@@ -232,5 +237,8 @@ QgsGeometry QgsInternalGeometryEngine::poleOfInaccessibility( double precision )
     cellQueue.push( new Cell( currentCell->x + h, currentCell->y + h, h, polygon ) );
   }
 
+  if ( distanceFromBoundary )
+    *distanceFromBoundary = bestCell->d;
+
   return QgsGeometry( new QgsPointV2( bestCell->x, bestCell->y ) );
 }
diff --git a/src/core/geometry/qgsinternalgeometryengine.h b/src/core/geometry/qgsinternalgeometryengine.h
@@ -57,8 +57,9 @@ class QgsInternalGeometryEngine
      * 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.
+     * Optionally, the distance to the polygon boundary from the pole can be stored.
      */
-    QgsGeometry poleOfInaccessibility( double precision ) const;
+    QgsGeometry poleOfInaccessibility( double precision, double* distanceFromBoundary = nullptr ) const;
 
   private:
     const QgsAbstractGeometry* mGeometry;

diff --git a/tests/src/core/testqgsgeometry.cpp b/tests/src/core/testqgsgeometry.cpp
@@ -3960,9 +3960,11 @@ void TestQgsGeometry::poleOfInaccessibility()
   QgsGeometry poly1 = QgsGeometry::fromWkt( poly1Wkt );
   QgsGeometry poly2 = QgsGeometry::fromWkt( poly2Wkt );
 
-  QgsPoint point = poly1.poleOfInaccessibility( 1 ).asPoint();
+  double distance;
+  QgsPoint point = poly1.poleOfInaccessibility( 1 , &distance ).asPoint();
   QGSCOMPARENEAR( point.x(), 3867.37, 0.01 );
   QGSCOMPARENEAR( point.y(), 2126.45, 0.01 );
+  QGSCOMPARENEAR( distance, 289.51, 0.01 );
 
   point = poly1.poleOfInaccessibility( 50 ).asPoint();
   QGSCOMPARENEAR( point.x(), 3855.33, 0.01 );