Skip to content

Commit

Permalink
Merge pull request #3843 from nyalldawson/oriented
Browse files Browse the repository at this point in the history 
Port minimum oriented bounding box to QgsGeometry
  • Loading branch information
@nyalldawson
nyalldawson committed Dec 7, 2016
2 parents ed5a2bc + 1bdb35d commit 37edb69
Show file tree
Hide file tree
Showing 11 changed files with 384 additions and 74 deletions.
14 changes: 13 additions & 1 deletion python/core/geometry/qgsgeometry.sip
View file
Expand Up @@ -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;

Expand Down
94 changes: 24 additions & 70 deletions python/plugins/processing/algs/qgis/OrientedMinimumBoundingBox.py
View file
Expand Up @@ -27,7 +27,7 @@

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
Expand Down Expand Up @@ -62,13 +62,15 @@ def processAlgorithm(self, progress):
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())
Expand All @@ -81,30 +83,27 @@ def processAlgorithm(self, progress):
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])
Expand All @@ -115,62 +114,17 @@ def featureOmbb(self, layer, writer, progress):
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
32 changes: 32 additions & 0 deletions python/plugins/processing/tests/testdata/custom/oriented_bbox.gfs
View file
@@ -0,0 +1,32 @@
<GMLFeatureClassList>
<GMLFeatureClass>
<Name>oriented_bbox</Name>
<ElementPath>oriented_bbox</ElementPath>
<!--POLYGON-->
<GeometryType>3</GeometryType>
<SRSName>EPSG:4326</SRSName>
<DatasetSpecificInfo>
<FeatureCount>6</FeatureCount>
<ExtentXMin>-1.00000</ExtentXMin>
<ExtentXMax>10.00000</ExtentXMax>
<ExtentYMin>-3.00000</ExtentYMin>
<ExtentYMax>6.00000</ExtentYMax>
</DatasetSpecificInfo>
<PropertyDefn>
<Name>intval</Name>
<ElementPath>intval</ElementPath>
<Type>Integer</Type>
</PropertyDefn>
<PropertyDefn>
<Name>floatval</Name>
<ElementPath>floatval</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>name</Name>
<ElementPath>name</ElementPath>
<Type>String</Type>
<Width>5</Width>
</PropertyDefn>
</GMLFeatureClass>
</GMLFeatureClassList>
59 changes: 59 additions & 0 deletions python/plugins/processing/tests/testdata/custom/oriented_bbox.gml
View file
@@ -0,0 +1,59 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation=""
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>-1</gml:X><gml:Y>-3</gml:Y></gml:coord>
<gml:coord><gml:X>10</gml:X><gml:Y>6</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>

<gml:featureMember>
<ogr:oriented_bbox fid="polys.4">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>3.8,2.2 6,1 6,-3 2.4,-1.0 3.8,2.2</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:intval>120</ogr:intval>
<ogr:floatval>-100291.43213</ogr:floatval>
</ogr:oriented_bbox>
</gml:featureMember>
<gml:featureMember>
<ogr:oriented_bbox fid="polys.1">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>5.4,5.0 6,4 5.2,3.8 4,4 5.4,5.0</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>Aaaaa</ogr:name>
<ogr:intval>-33</ogr:intval>
<ogr:floatval>0</ogr:floatval>
</ogr:oriented_bbox>
</gml:featureMember>
<gml:featureMember>
<ogr:oriented_bbox fid="polys.0">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>-1,-1 -1,3 3,3 3,2 2,2 2,-1 -1,-1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>aaaaa</ogr:name>
<ogr:intval>33</ogr:intval>
<ogr:floatval>44.12346</ogr:floatval>
</ogr:oriented_bbox>
</gml:featureMember>
<gml:featureMember>
<ogr:oriented_bbox fid="polys.3">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>6.8,1.8 10,1 9.6,-2.2 6.4,-3.0 7.2,-0.6 6.8,1.8</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs><gml:innerBoundaryIs><gml:LinearRing><gml:coordinates>8.0,-0.6 7,-2 9,-2 9,0 8.0,-0.6</gml:coordinates></gml:LinearRing></gml:innerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>ASDF</ogr:name>
<ogr:intval>0</ogr:intval>
</ogr:oriented_bbox>
</gml:featureMember>
<gml:featureMember>
<ogr:oriented_bbox fid="polys.2">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>1.6,4.8 2,6 3,6 1.6,4.8</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>bbaaa</ogr:name>
<ogr:floatval>0.123</ogr:floatval>
</ogr:oriented_bbox>
</gml:featureMember>
<gml:featureMember>
<ogr:oriented_bbox fid="polys.5">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>3.8,2.2 6,1 6,-3 2.4,-1.0 3.8,2.2</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>elim</ogr:name>
<ogr:intval>2</ogr:intval>
<ogr:floatval>3.33</ogr:floatval>
</ogr:oriented_bbox>
</gml:featureMember>
</ogr:FeatureCollection>
57 changes: 57 additions & 0 deletions python/plugins/processing/tests/testdata/expected/oriented_bounds.gfs
View file
@@ -0,0 +1,57 @@
<GMLFeatureClassList>
<GMLFeatureClass>
<Name>oriented_bounds</Name>
<ElementPath>oriented_bounds</ElementPath>
<!--POLYGON-->
<GeometryType>3</GeometryType>
<SRSName>EPSG:4326</SRSName>
<DatasetSpecificInfo>
<FeatureCount>6</FeatureCount>
<ExtentXMin>-1.00000</ExtentXMin>
<ExtentXMax>10.04414</ExtentXMax>
<ExtentYMin>-3.27034</ExtentYMin>
<ExtentYMax>6.44118</ExtentYMax>
</DatasetSpecificInfo>
<PropertyDefn>
<Name>intval</Name>
<ElementPath>intval</ElementPath>
<Type>Integer</Type>
</PropertyDefn>
<PropertyDefn>
<Name>floatval</Name>
<ElementPath>floatval</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>area</Name>
<ElementPath>area</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>perimeter</Name>
<ElementPath>perimeter</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>angle</Name>
<ElementPath>angle</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>width</Name>
<ElementPath>width</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>height</Name>
<ElementPath>height</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>name</Name>
<ElementPath>name</ElementPath>
<Type>String</Type>
<Width>5</Width>
</PropertyDefn>
</GMLFeatureClass>
</GMLFeatureClassList>

0 comments on commit 37edb69

Please sign in to comment.