Skip to content

Commit

Permalink
[FEATURE] Subdivide algorithm for QgsGeometry
Browse files Browse the repository at this point in the history
Subdivides the geometry. The returned geometry will be a collection
containing subdivided parts from the original geometry, where no
part has more then the specified maximum number of nodes.
  • Loading branch information
nyalldawson committed Jun 14, 2017
1 parent bde0c72 commit e74395d
Show file tree
Hide file tree
Showing 6 changed files with 210 additions and 0 deletions.
17 changes: 17 additions & 0 deletions python/core/geometry/qgsgeometry.sip
Original file line number Diff line number Diff line change
Expand Up @@ -782,6 +782,23 @@ Returns the smallest convex polygon that contains all the points in the geometry
If ``edgesOnly`` is true than line string boundary geometries will be returned
instead of polygons.
An empty geometry will be returned if the diagram could not be calculated.
.. versionadded:: 3.0
:rtype: QgsGeometry
%End

QgsGeometry subdivide( int maxNodes = 256 ) const;
%Docstring
Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
from the original geometry, where no part has more then the specified maximum number of nodes (``maxNodes``).

This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
indexed and faster to perform further operations such as intersects on. The returned geometry parts may
not be valid and may contain self-intersections.

The minimum allowed value for ``maxNodes`` is 8.

Curved geometries will be segmentized before subdivision.

.. versionadded:: 3.0
:rtype: QgsGeometry
%End
Expand Down
24 changes: 24 additions & 0 deletions src/core/geometry/qgsgeometry.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1655,6 +1655,30 @@ QgsGeometry QgsGeometry::delaunayTriangulation( double tolerance, bool edgesOnly
return geos.delaunayTriangulation( tolerance, edgesOnly );
}

QgsGeometry QgsGeometry::subdivide( int maxNodes ) const
{
if ( !d->geometry )
{
return QgsGeometry();
}

const QgsAbstractGeometry *geom = d->geometry;
std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
if ( QgsWkbTypes::isCurvedType( d->geometry->wkbType() ) )
{
segmentizedCopy.reset( d->geometry->segmentize() );
geom = segmentizedCopy.get();
}

QgsGeos geos( geom );
QgsAbstractGeometry *result = geos.subdivide( maxNodes );
if ( !result )
{
return QgsGeometry();
}
return QgsGeometry( result );
}

QgsGeometry QgsGeometry::interpolate( double distance ) const
{
if ( !d->geometry )
Expand Down
16 changes: 16 additions & 0 deletions src/core/geometry/qgsgeometry.h
Original file line number Diff line number Diff line change
Expand Up @@ -733,6 +733,22 @@ class CORE_EXPORT QgsGeometry
*/
QgsGeometry delaunayTriangulation( double tolerance = 0.0, bool edgesOnly = false ) const;

/**
* Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
* from the original geometry, where no part has more then the specified maximum number of nodes (\a maxNodes).
*
* This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
* indexed and faster to perform further operations such as intersects on. The returned geometry parts may
* not be valid and may contain self-intersections.
*
* The minimum allowed value for \a maxNodes is 8.
*
* Curved geometries will be segmentized before subdivision.
*
* \since QGIS 3.0
*/
QgsGeometry subdivide( int maxNodes = 256 ) const;

/**
* Return interpolated point on line at distance
* \since QGIS 1.9
Expand Down
116 changes: 116 additions & 0 deletions src/core/geometry/qgsgeos.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,122 @@ QgsAbstractGeometry *QgsGeos::clip( const QgsRectangle &rect, QString *errorMsg
}
}




void QgsGeos::subdivideRecursive( const GEOSGeometry *currentPart, int maxNodes, int depth, QgsGeometryCollection *parts, const QgsRectangle &clipRect ) const
{
int partType = GEOSGeomTypeId_r( geosinit.ctxt, currentPart );
if ( qgsDoubleNear( clipRect.width(), 0.0 ) && qgsDoubleNear( clipRect.height(), 0.0 ) )
{
if ( partType == GEOS_POINT )
{
parts->addGeometry( fromGeos( currentPart ) );
return;
}
else
{
return;
}
}

if ( partType == GEOS_MULTILINESTRING || partType == GEOS_MULTIPOLYGON || partType == GEOS_GEOMETRYCOLLECTION )
{
int partCount = GEOSGetNumGeometries_r( geosinit.ctxt, currentPart );
for ( int i = 0; i < partCount; ++i )
{
subdivideRecursive( GEOSGetGeometryN_r( geosinit.ctxt, currentPart, i ), maxNodes, depth, parts, clipRect );
}
return;
}

if ( depth > 50 )
{
parts->addGeometry( fromGeos( currentPart ) );
return;
}

int vertexCount = GEOSGetNumCoordinates_r( geosinit.ctxt, currentPart );
if ( vertexCount == 0 )
{
return;
}
else if ( vertexCount < maxNodes )
{
parts->addGeometry( fromGeos( currentPart ) );
return;
}

// chop clipping rect in half by longest side
double width = clipRect.width();
double height = clipRect.height();
QgsRectangle halfClipRect1 = clipRect;
QgsRectangle halfClipRect2 = clipRect;
if ( width > height )
{
halfClipRect1.setXMaximum( clipRect.xMinimum() + width / 2.0 );
halfClipRect2.setXMinimum( halfClipRect1.xMaximum() );
}
else
{
halfClipRect1.setYMaximum( clipRect.yMinimum() + height / 2.0 );
halfClipRect2.setYMinimum( halfClipRect1.yMaximum() );
}

if ( height <= 0 )
{
halfClipRect1.setYMinimum( halfClipRect1.yMinimum() - DBL_EPSILON );
halfClipRect2.setYMinimum( halfClipRect2.yMinimum() - DBL_EPSILON );
halfClipRect1.setYMaximum( halfClipRect1.yMaximum() + DBL_EPSILON );
halfClipRect2.setYMaximum( halfClipRect2.yMaximum() + DBL_EPSILON );
}
if ( width <= 0 )
{
halfClipRect1.setXMinimum( halfClipRect1.xMinimum() - DBL_EPSILON );
halfClipRect2.setXMinimum( halfClipRect2.xMinimum() - DBL_EPSILON );
halfClipRect1.setXMaximum( halfClipRect1.xMaximum() + DBL_EPSILON );
halfClipRect2.setXMaximum( halfClipRect2.xMaximum() + DBL_EPSILON );
}

GEOSGeomScopedPtr clipPart1;
clipPart1.reset( GEOSClipByRect_r( geosinit.ctxt, currentPart, halfClipRect1.xMinimum(), halfClipRect1.yMinimum(), halfClipRect1.xMaximum(), halfClipRect1.yMaximum() ) );
GEOSGeomScopedPtr clipPart2;
clipPart2.reset( GEOSClipByRect_r( geosinit.ctxt, currentPart, halfClipRect2.xMinimum(), halfClipRect2.yMinimum(), halfClipRect2.xMaximum(), halfClipRect2.yMaximum() ) );

++depth;

if ( clipPart1 )
{
subdivideRecursive( clipPart1.get(), maxNodes, depth, parts, halfClipRect1 );
}
if ( clipPart2 )
{
subdivideRecursive( clipPart2.get(), maxNodes, depth, parts, halfClipRect2 );
}

return;
}

QgsAbstractGeometry *QgsGeos::subdivide( int maxNodes, QString *errorMsg ) const
{
if ( !mGeos )
{
return nullptr;
}

// minimum allowed max is 8
maxNodes = qMax( maxNodes, 8 );

std::unique_ptr< QgsGeometryCollection > parts = QgsGeometryFactory::createCollectionOfType( mGeometry->wkbType() );
try
{
subdivideRecursive( mGeos, maxNodes, 0, parts.get(), mGeometry->boundingBox() );
}
CATCH_GEOS_WITH_ERRMSG( nullptr )

return parts.release();
}

QgsAbstractGeometry *QgsGeos::combine( const QgsAbstractGeometry &geom, QString *errorMsg ) const
{
return overlay( geom, UNION, errorMsg );
Expand Down
18 changes: 18 additions & 0 deletions src/core/geometry/qgsgeos.h
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ email : marco.hugentobler at sourcepole dot com
class QgsLineString;
class QgsPolygonV2;
class QgsGeometry;
class QgsGeometryCollection;

/** \ingroup core
* Does vector analysis using the geos library and handles import, export, exception handling*
Expand Down Expand Up @@ -52,6 +53,22 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
*/
QgsAbstractGeometry *clip( const QgsRectangle &rectangle, QString *errorMsg = nullptr ) const;

/**
* Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
* from the original geometry, where no part has more then the specified maximum number of nodes (\a maxNodes).
*
* This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
* indexed and faster to perform further operations such as intersects on. The returned geometry parts may
* not be valid and may contain self-intersections.
*
* The minimum allowed value for \a maxNodes is 8.
*
* Curved geometries are not supported.
*
* \since QGIS 3.0
*/
QgsAbstractGeometry *subdivide( int maxNodes, QString *errorMsg = nullptr ) const;

QgsAbstractGeometry *combine( const QgsAbstractGeometry &geom, QString *errorMsg = nullptr ) const override;
QgsAbstractGeometry *combine( const QList< QgsAbstractGeometry *> &, QString *errorMsg = nullptr ) const override;
QgsAbstractGeometry *symDifference( const QgsAbstractGeometry &geom, QString *errorMsg = nullptr ) const override;
Expand Down Expand Up @@ -253,6 +270,7 @@ class CORE_EXPORT QgsGeos: public QgsGeometryEngine
static int lineContainedInLine( const GEOSGeometry *line1, const GEOSGeometry *line2 );
static int pointContainedInLine( const GEOSGeometry *point, const GEOSGeometry *line );
static int geomDigits( const GEOSGeometry *geom );
void subdivideRecursive( const GEOSGeometry *currentPart, int maxNodes, int depth, QgsGeometryCollection *parts, const QgsRectangle &clipRect ) const;
};

/// @cond PRIVATE
Expand Down
19 changes: 19 additions & 0 deletions tests/src/python/test_qgsgeometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -4165,6 +4165,25 @@ def testPoint(self):
self.assertEqual(point_zm.z(), 3)
self.assertEqual(point_zm.m(), 4)

def testSubdivide(self):
tests = [["LINESTRING (1 1,1 9,9 9,9 1)", 8, "MULTILINESTRING ((1 1,1 9,9 9,9 1))"],
["Point (1 1)", 8, "MultiPoint ((1 1))"],
["GeometryCollection ()", 8, "GeometryCollection ()"],
["LINESTRING (1 1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1 9)", 8, "MultiLineString ((1 1, 1 2, 1 3, 1 4, 1 5),(1 5, 1 6, 1 7, 1 8, 1 9))"],
["LINESTRING (1 1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1 9)", 1,
"MultiLineString ((1 1, 1 2, 1 3, 1 4, 1 5),(1 5, 1 6, 1 7, 1 8, 1 9))"],
["LINESTRING (1 1,1 2,1 3,1 4,1 5,1 6,1 7,1 8,1 9)", 16,
"MultiLineString ((1 1, 1 2, 1 3, 1 4, 1 5, 1 6, 1 7, 1 8, 1 9))"],
["POLYGON ((0 0, 0 200, 200 200, 200 0, 0 0),(60 180, 20 180, 20 140, 60 140, 60 180),(180 60, 140 60, 140 20, 180 20, 180 60))", 8, "MultiPolygon (((0 0, 0 100, 100 100, 100 0, 0 0)),((100 0, 100 50, 140 50, 140 20, 150 20, 150 0, 100 0)),((150 0, 150 20, 180 20, 180 50, 200 50, 200 0, 150 0)),((100 50, 100 100, 150 100, 150 60, 140 60, 140 50, 100 50)),((150 60, 150 100, 200 100, 200 50, 180 50, 180 60, 150 60)),((0 100, 0 150, 20 150, 20 140, 50 140, 50 100, 0 100)),((50 100, 50 140, 60 140, 60 150, 100 150, 100 100, 50 100)),((0 150, 0 200, 50 200, 50 180, 20 180, 20 150, 0 150)),((50 180, 50 200, 100 200, 100 150, 60 150, 60 180, 50 180)),((100 100, 100 200, 200 200, 200 100, 100 100)))"]
]
for t in tests:
input = QgsGeometry.fromWkt(t[0])
o = input.subdivide(t[1])
exp = t[2]
result = o.exportToWkt()
self.assertTrue(compareWkt(result, exp, 0.00001),
"clipped: mismatch Expected:\n{}\nGot:\n{}\n".format(exp, result))

def testClipped(self):
tests = [["LINESTRING (1 1,1 9,9 9,9 1)", QgsRectangle(0, 0, 10, 10), "LINESTRING (1 1,1 9,9 9,9 1)"],
["LINESTRING (-1 -9,-1 11,9 11)", QgsRectangle(0, 0, 10, 10), "GEOMETRYCOLLECTION ()"],
Expand Down

0 comments on commit e74395d

Please sign in to comment.