From c975764c12a182a9be37c1cce125ecc924da4e5f Mon Sep 17 00:00:00 2001 From: Nyall Dawson Date: Wed, 7 Dec 2016 16:41:56 +1000 Subject: [PATCH 1/2] Port processing oriented minimum bounding box alg to QgsGeometry --- python/core/geometry/qgsgeometry.sip | 14 ++- .../algs/qgis/OrientedMinimumBoundingBox.py | 94 +++++-------------- .../tests/testdata/custom/oriented_bbox.gfs | 32 +++++++ .../tests/testdata/custom/oriented_bbox.gml | 59 ++++++++++++ .../testdata/expected/oriented_bounds.gfs | 57 +++++++++++ .../testdata/expected/oriented_bounds.gml | 89 ++++++++++++++++++ .../tests/testdata/qgis_algorithm_tests.yaml | 12 +++ src/core/geometry/qgsgeometry.cpp | 59 ++++++++++++ src/core/geometry/qgsgeometry.h | 14 ++- tests/src/python/test_qgsgeometry.py | 23 +++++ 10 files changed, 381 insertions(+), 72 deletions(-) create mode 100644 python/plugins/processing/tests/testdata/custom/oriented_bbox.gfs create mode 100644 python/plugins/processing/tests/testdata/custom/oriented_bbox.gml create mode 100644 python/plugins/processing/tests/testdata/expected/oriented_bounds.gfs create mode 100644 python/plugins/processing/tests/testdata/expected/oriented_bounds.gml diff --git a/python/core/geometry/qgsgeometry.sip b/python/core/geometry/qgsgeometry.sip index e752c913dad..0193e0f36bf 100644 --- a/python/core/geometry/qgsgeometry.sip +++ b/python/core/geometry/qgsgeometry.sip @@ -405,9 +405,21 @@ class QgsGeometry */ QgsGeometry makeDifference( const QgsGeometry& other ) const; - /** Returns the bounding box of this feature*/ + /** + * Returns the bounding box of the geometry. + * @see orientedMinimumBoundingBox() + */ QgsRectangle boundingBox() const; + /** + * Returns the oriented minimum bounding box for the geometry, which is the smallest (by area) + * rotated rectangle which fully encompasses the geometry. The area, angle (clockwise in degrees from North), + * width and height of the rotated bounding box will also be returned. + * @note added in QGIS 3.0 + * @see boundingBox() + */ + QgsGeometry orientedMinimumBoundingBox( double& area /Out/, double &angle /Out/, double& width /Out/, double& height /Out/ ) const; + /** Test for intersection with a rectangle (uses GEOS) */ bool intersects( const QgsRectangle& r ) const; diff --git a/python/plugins/processing/algs/qgis/OrientedMinimumBoundingBox.py b/python/plugins/processing/algs/qgis/OrientedMinimumBoundingBox.py index febebe96737..c7e00bac339 100644 --- a/python/plugins/processing/algs/qgis/OrientedMinimumBoundingBox.py +++ b/python/plugins/processing/algs/qgis/OrientedMinimumBoundingBox.py @@ -27,7 +27,7 @@ __revision__ = '$Format:%H$' from math import degrees, atan2 from qgis.PyQt.QtCore import QVariant -from qgis.core import Qgis, QgsField, QgsPoint, QgsGeometry, QgsFeature, QgsWkbTypes +from qgis.core import Qgis, QgsField, QgsFields, QgsPoint, QgsGeometry, QgsFeature, QgsWkbTypes, QgsFeatureRequest from processing.core.GeoAlgorithm import GeoAlgorithm from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException from processing.core.parameters import ParameterVector @@ -62,13 +62,15 @@ class OrientedMinimumBoundingBox(GeoAlgorithm): if byFeature and layer.geometryType() == QgsWkbTypes.PointGeometry and layer.featureCount() <= 2: raise GeoAlgorithmExecutionException(self.tr("Can't calculate an OMBB for each point, it's a point. The number of points must be greater than 2")) - fields = [ - QgsField('AREA', QVariant.Double), - QgsField('PERIMETER', QVariant.Double), - QgsField('ANGLE', QVariant.Double), - QgsField('WIDTH', QVariant.Double), - QgsField('HEIGHT', QVariant.Double), - ] + if byFeature: + fields = layer.fields() + else: + fields = QgsFields() + fields.append(QgsField('area', QVariant.Double)) + fields.append(QgsField('perimeter', QVariant.Double)) + fields.append(QgsField('angle', QVariant.Double)) + fields.append(QgsField('width', QVariant.Double)) + fields.append(QgsField('height', QVariant.Double)) writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(fields, QgsWkbTypes.Polygon, layer.crs()) @@ -81,30 +83,27 @@ class OrientedMinimumBoundingBox(GeoAlgorithm): del writer def layerOmmb(self, layer, writer, progress): - current = 0 - - fit = layer.getFeatures() - inFeat = QgsFeature() - total = 100.0 / layer.featureCount() + req = QgsFeatureRequest().setSubsetOfAttributes([]) + features = vector.features(layer, req) + total = 100.0 / len(features) newgeometry = QgsGeometry() first = True - while fit.nextFeature(inFeat): + for current, inFeat in enumerate(features): if first: newgeometry = inFeat.geometry() first = False else: newgeometry = newgeometry.combine(inFeat.geometry()) - current += 1 progress.setPercentage(int(current * total)) - geometry, area, perim, angle, width, height = self.OMBBox(newgeometry) + geometry, area, angle, width, height = newgeometry.orientedMinimumBoundingBox() if geometry: outFeat = QgsFeature() outFeat.setGeometry(geometry) outFeat.setAttributes([area, - perim, + width * 2 + height * 2, angle, width, height]) @@ -115,62 +114,17 @@ class OrientedMinimumBoundingBox(GeoAlgorithm): total = 100.0 / len(features) outFeat = QgsFeature() for current, inFeat in enumerate(features): - geometry, area, perim, angle, width, height = self.OMBBox( - inFeat.geometry()) + geometry, area, angle, width, height = inFeat.geometry().orientedMinimumBoundingBox() if geometry: outFeat.setGeometry(geometry) - outFeat.setAttributes([area, - perim, - angle, - width, - height]) + attrs = inFeat.attributes() + attrs.extend([area, + width * 2 + height * 2, + angle, + width, + height]) + outFeat.setAttributes(attrs) writer.addFeature(outFeat) else: progress.setInfo(self.tr("Can't calculate an OMBB for feature {0}.").format(inFeat.id())) progress.setPercentage(int(current * total)) - - def GetAngleOfLineBetweenTwoPoints(self, p1, p2, angle_unit="degrees"): - xDiff = p2.x() - p1.x() - yDiff = p2.y() - p1.y() - if angle_unit == "radians": - return atan2(yDiff, xDiff) - else: - return degrees(atan2(yDiff, xDiff)) - - def OMBBox(self, geom): - g = geom.convexHull() - - if g.type() != QgsWkbTypes.PolygonGeometry: - return None, None, None, None, None, None - r = g.asPolygon()[0] - - p0 = QgsPoint(r[0][0], r[0][1]) - - i = 0 - l = len(r) - OMBBox = QgsGeometry() - gBBox = g.boundingBox() - OMBBox_area = gBBox.height() * gBBox.width() - OMBBox_angle = 0 - OMBBox_width = 0 - OMBBox_heigth = 0 - OMBBox_perim = 0 - while i < l - 1: - x = QgsGeometry(g) - angle = self.GetAngleOfLineBetweenTwoPoints(r[i], r[i + 1]) - x.rotate(angle, p0) - bbox = x.boundingBox() - bb = QgsGeometry.fromWkt(bbox.asWktPolygon()) - bb.rotate(-angle, p0) - - areabb = bb.area() - if areabb <= OMBBox_area: - OMBBox = QgsGeometry(bb) - OMBBox_area = areabb - OMBBox_angle = angle - OMBBox_width = bbox.width() - OMBBox_heigth = bbox.height() - OMBBox_perim = 2 * OMBBox_width + 2 * OMBBox_heigth - i += 1 - - return OMBBox, OMBBox_area, OMBBox_perim, OMBBox_angle, OMBBox_width, OMBBox_heigth diff --git a/python/plugins/processing/tests/testdata/custom/oriented_bbox.gfs b/python/plugins/processing/tests/testdata/custom/oriented_bbox.gfs new file mode 100644 index 00000000000..f5795f83339 --- /dev/null +++ b/python/plugins/processing/tests/testdata/custom/oriented_bbox.gfs @@ -0,0 +1,32 @@ + + + oriented_bbox + oriented_bbox + + 3 + EPSG:4326 + + 6 + -1.00000 + 10.00000 + -3.00000 + 6.00000 + + + intval + intval + Integer + + + floatval + floatval + Real + + + name + name + String + 5 + + + diff --git a/python/plugins/processing/tests/testdata/custom/oriented_bbox.gml b/python/plugins/processing/tests/testdata/custom/oriented_bbox.gml new file mode 100644 index 00000000000..50bf756fb1c --- /dev/null +++ b/python/plugins/processing/tests/testdata/custom/oriented_bbox.gml @@ -0,0 +1,59 @@ + + + + + -1-3 + 106 + + + + + + 3.8,2.2 6,1 6,-3 2.4,-1.0 3.8,2.2 + 120 + -100291.43213 + + + + + 5.4,5.0 6,4 5.2,3.8 4,4 5.4,5.0 + Aaaaa + -33 + 0 + + + + + -1,-1 -1,3 3,3 3,2 2,2 2,-1 -1,-1 + aaaaa + 33 + 44.12346 + + + + + 6.8,1.8 10,1 9.6,-2.2 6.4,-3.0 7.2,-0.6 6.8,1.88.0,-0.6 7,-2 9,-2 9,0 8.0,-0.6 + ASDF + 0 + + + + + 1.6,4.8 2,6 3,6 1.6,4.8 + bbaaa + 0.123 + + + + + 3.8,2.2 6,1 6,-3 2.4,-1.0 3.8,2.2 + elim + 2 + 3.33 + + + diff --git a/python/plugins/processing/tests/testdata/expected/oriented_bounds.gfs b/python/plugins/processing/tests/testdata/expected/oriented_bounds.gfs new file mode 100644 index 00000000000..2867f28fd42 --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/oriented_bounds.gfs @@ -0,0 +1,57 @@ + + + oriented_bounds + oriented_bounds + + 3 + EPSG:4326 + + 6 + -1.00000 + 10.04414 + -3.27034 + 6.44118 + + + intval + intval + Integer + + + floatval + floatval + Real + + + area + area + Real + + + perimeter + perimeter + Real + + + angle + angle + Real + + + width + width + Real + + + height + height + Real + + + name + name + String + 5 + + + diff --git a/python/plugins/processing/tests/testdata/expected/oriented_bounds.gml b/python/plugins/processing/tests/testdata/expected/oriented_bounds.gml new file mode 100644 index 00000000000..8f543a5b0ff --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/oriented_bounds.gml @@ -0,0 +1,89 @@ + + + + + -1-3.270344827586206 + 10.044137931034486.441176470588236 + + + + + + 6,-3 7.69811320754717,0.056603773584906 3.80943396226415,2.21698113207547 2.11132075471698,-0.839622641509436 6,-3 + 120 + -100291.43213 + 15.5547169811321 + 15.8902367081648 + 119.054604099077 + 3.49662910448615 + 4.44848924959627 + + + + + 4.11764705882353,3.52941176470588 6.0,4.0 5.72941176470588,5.08235294117647 3.84705882352941,4.61176470588235 4.11764705882353,3.52941176470588 + -33 + 0 + Aaaaa + 2.16470588235294 + 6.11189775091559 + 165.963756532073 + 1.94028500029066 + 1.11566387516713 + + + + + -1.0,3.0 -1.0,-1.0 3.0,-1.0 3.0,3.0 -1.0,3.0 + 33 + 44.12346 + aaaaa + 16 + 16 + 90 + 4 + 4 + + + + + 6.4,-3.0 9.64413793103449,-3.27034482758621 10.0441379310345,1.52965517241379 6.8,1.8 6.4,-3.0 + 0 + ASDF + 15.68 + 16.1440412835671 + 4.76364169072617 + 3.25538281026661 + 4.81663783151692 + + + + + 1.36470588235294,4.94117647058824 2.1,4.5 3.0,6.0 2.26470588235294,6.44117647058824 1.36470588235294,4.94117647058824 + 0.123 + bbaaa + 1.5 + 5.21355698833227 + 30.9637565320735 + 0.857492925712544 + 1.74928556845359 + + + + + 6,-3 7.69811320754717,0.056603773584906 3.80943396226415,2.21698113207547 2.11132075471698,-0.839622641509436 6,-3 + 2 + 3.33 + elim + 15.5547169811321 + 15.8902367081648 + 119.054604099077 + 3.49662910448615 + 4.44848924959627 + + + diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml index 57924fa0a51..ad1d74b83d7 100644 --- a/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml +++ b/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml @@ -1854,3 +1854,15 @@ tests: OUTPUT: hash: fe6e018be13c5a3c17f3f4d0f0dc7686c628cb440b74c4642aa0c939 type: rasterhash + + - algorithm: qgis:orientedminimumboundingbox + name: Oriented minimum bounding box polys + params: + BY_FEATURE: true + INPUT_LAYER: + name: custom/oriented_bbox.gml + type: vector + results: + OUTPUT: + name: expected/oriented_bounds.gml + type: vector diff --git a/src/core/geometry/qgsgeometry.cpp b/src/core/geometry/qgsgeometry.cpp index 48ad9e9dd99..a8984e50289 100644 --- a/src/core/geometry/qgsgeometry.cpp +++ b/src/core/geometry/qgsgeometry.cpp @@ -862,6 +862,65 @@ QgsRectangle QgsGeometry::boundingBox() const return QgsRectangle(); } +QgsGeometry QgsGeometry::orientedMinimumBoundingBox( double& area, double &angle, double& width, double& height ) const +{ + QgsRectangle minRect; + area = DBL_MAX; + angle = 0; + width = DBL_MAX; + height = DBL_MAX; + + if ( !d->geometry || d->geometry->nCoordinates() < 2 ) + return QgsGeometry(); + + QgsGeometry hull = convexHull(); + if ( hull.isEmpty() ) + return QgsGeometry(); + + QgsVertexId vertexId; + QgsPointV2 pt0; + QgsPointV2 pt1; + QgsPointV2 pt2; + // get first point + hull.geometry()->nextVertex( vertexId, pt0 ); + pt1 = pt0; + double prevAngle = 0.0; + while ( hull.geometry()->nextVertex( vertexId, pt2 ) ) + { + double currentAngle = QgsGeometryUtils::lineAngle( pt1.x(), pt1.y(), pt2.x(), pt2.y() ); + double rotateAngle = 180.0 / M_PI * ( currentAngle - prevAngle ); + prevAngle = currentAngle; + + QTransform t = QTransform::fromTranslate( pt0.x(), pt0.y() ); + t.rotate( rotateAngle ); + t.translate( -pt0.x(), -pt0.y() ); + + hull.geometry()->transform( t ); + + QgsRectangle bounds = hull.geometry()->boundingBox(); + double currentArea = bounds.width() * bounds.height(); + if ( currentArea < area ) + { + minRect = bounds; + area = currentArea; + angle = 180.0 / M_PI * currentAngle; + width = bounds.width(); + height = bounds.height(); + } + + pt2 = pt1; + } + + QgsGeometry minBounds = QgsGeometry::fromRect( minRect ); + minBounds.rotate( angle, QgsPoint( pt0.x(), pt0.y() ) ); + + // constrain angle to 0 - 180 + if ( angle > 180.0 ) + angle = fmod( angle, 180.0 ); + + return minBounds; +} + bool QgsGeometry::intersects( const QgsRectangle& r ) const { QgsGeometry g = fromRect( r ); diff --git a/src/core/geometry/qgsgeometry.h b/src/core/geometry/qgsgeometry.h index 4cee39704cd..9c0c33ce54f 100644 --- a/src/core/geometry/qgsgeometry.h +++ b/src/core/geometry/qgsgeometry.h @@ -458,9 +458,21 @@ class CORE_EXPORT QgsGeometry */ QgsGeometry makeDifference( const QgsGeometry& other ) const; - //! Returns the bounding box of this feature + /** + * Returns the bounding box of the geometry. + * @see orientedMinimumBoundingBox() + */ QgsRectangle boundingBox() const; + /** + * Returns the oriented minimum bounding box for the geometry, which is the smallest (by area) + * rotated rectangle which fully encompasses the geometry. The area, angle (clockwise in degrees from North), + * width and height of the rotated bounding box will also be returned. + * @note added in QGIS 3.0 + * @see boundingBox() + */ + QgsGeometry orientedMinimumBoundingBox( double& area, double &angle, double& width, double& height ) const; + //! Test for intersection with a rectangle (uses GEOS) bool intersects( const QgsRectangle& r ) const; diff --git a/tests/src/python/test_qgsgeometry.py b/tests/src/python/test_qgsgeometry.py index 6f30dd2af21..6bf534ce97b 100644 --- a/tests/src/python/test_qgsgeometry.py +++ b/tests/src/python/test_qgsgeometry.py @@ -3645,6 +3645,29 @@ class TestQgsGeometry(unittest.TestCase): self.assertTrue(compareWkt(result, exp, 0.00001), "Extend line: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result)) + def testMinimumOrientedBoundingBox(self): + empty = QgsGeometry() + bbox, area, angle, width, height = empty.orientedMinimumBoundingBox() + self.assertFalse(bbox) + + # not a useful geometry + point = QgsGeometry.fromWkt('Point(1 2)') + bbox, area, angle, width, height = point.orientedMinimumBoundingBox() + self.assertFalse(bbox) + + # polygon + polygon = QgsGeometry.fromWkt('Polygon((-0.1 -1.3, 2.1 1, 3 2.8, 6.7 0.2, 3 -1.8, 0.3 -2.7, -0.1 -1.3))') + bbox, area, angle, width, height = polygon.orientedMinimumBoundingBox() + exp = 'Polygon ((-0.94905660 -1.571698, 2.3817055 -4.580453, 6.7000000 0.1999999, 3.36923 3.208754, -0.949056 -1.57169))' + + result = bbox.exportToWkt() + self.assertTrue(compareWkt(result, exp, 0.00001), + "Oriented MBBR: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result)) + self.assertAlmostEqual(area, 28.9152, places=3) + self.assertAlmostEqual(angle, 42.0922, places=3) + self.assertAlmostEqual(width, 4.4884, places=3) + self.assertAlmostEqual(height, 6.4420, places=3) + if __name__ == '__main__': unittest.main() From 1bdb35d63044c24df406d0b7b9b888f84642e705 Mon Sep 17 00:00:00 2001 From: Nyall Dawson Date: Wed, 7 Dec 2016 19:57:07 +1000 Subject: [PATCH 2/2] Avoid key error on fields which should be skipped --- python/testing/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/testing/__init__.py b/python/testing/__init__.py index 7289292f75d..74d6b808000 100644 --- a/python/testing/__init__.py +++ b/python/testing/__init__.py @@ -102,8 +102,6 @@ class TestCase(_TestCase): ) for attr_expected, field_expected in zip(feats[0].attributes(), layer_expected.fields().toList()): - attr_result = feats[1][field_expected.name()] - field_result = [fld for fld in layer_expected.fields().toList() if fld.name() == field_expected.name()][0] try: cmp = compare['fields'][field_expected.name()] except KeyError: @@ -116,6 +114,9 @@ class TestCase(_TestCase): if 'skip' in cmp: continue + attr_result = feats[1][field_expected.name()] + field_result = [fld for fld in layer_expected.fields().toList() if fld.name() == field_expected.name()][0] + # Cast field to a given type if 'cast' in cmp: if cmp['cast'] == 'int':