QGIS API Documentation 3.39.0-Master (67e056379ed)
Loading...
Searching...
No Matches
qgsalgorithmconcavehull.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmconcavehull.cpp
3 ---------------------
4 begin : July 2023
5 copyright : (C) 2023 by Alexander Bruy
6 email : alexander dot bruy at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
19#include "qgsspatialindex.h"
20#include "qgsmultipoint.h"
22#include "qgsvectorlayer.h"
23
25
26QString QgsConcaveHullAlgorithm::name() const
27{
28 return QStringLiteral( "concavehull" );
29}
30
31QString QgsConcaveHullAlgorithm::displayName() const
32{
33 return QObject::tr( "Concave hull" );
34}
35
36QStringList QgsConcaveHullAlgorithm::tags() const
37{
38 return QObject::tr( "concave,hull,bounds,bounding" ).split( ',' );
39}
40
41QString QgsConcaveHullAlgorithm::group() const
42{
43 return QObject::tr( "Vector geometry" );
44}
45
46QString QgsConcaveHullAlgorithm::groupId() const
47{
48 return QStringLiteral( "vectorgeometry" );
49}
50
51QString QgsConcaveHullAlgorithm::shortHelpString() const
52{
53 return QObject::tr( "This algorithm computes the concave hull of the features from an input layer." );
54}
55
56QgsConcaveHullAlgorithm *QgsConcaveHullAlgorithm::createInstance() const
57{
58 return new QgsConcaveHullAlgorithm();
59}
60
61void QgsConcaveHullAlgorithm::initAlgorithm( const QVariantMap & )
62{
63 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList< int >() << static_cast< int >( Qgis::ProcessingSourceType::VectorPoint ) ) );
64 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "ALPHA" ), QObject::tr( "Threshold (0-1, where 1 is equivalent with Convex Hull)" ), Qgis::ProcessingNumberParameterType::Double, 0.3, false, 0, 1 ) );
65 addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "HOLES" ), QObject::tr( "Allow holes" ), true ) );
66 addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "NO_MULTIGEOMETRY" ), QObject::tr( "Split multipart geometry into singleparts" ), false ) );
67 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Concave hull" ), Qgis::ProcessingSourceType::VectorPolygon ) );
68}
69
70bool QgsConcaveHullAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
71{
72 mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
73 if ( !mSource )
74 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
75
76 if ( mSource->featureCount() < 3 )
77 throw QgsProcessingException( QObject::tr( "Input layer should contain at least 3 points." ) );
78
79 mStep = mSource->featureCount() > 0 ? 50.0 / mSource->featureCount() : 1;
80
81 mPercentage = parameterAsDouble( parameters, QStringLiteral( "ALPHA" ), context );
82 mAllowHoles = parameterAsBool( parameters, QStringLiteral( "HOLES" ), context );
83 mSplitMultipart = parameterAsBool( parameters, QStringLiteral( "NO_MULTIGEOMETRY" ), context );
84
85 return true;
86}
87
88QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
89{
90 QString dest;
91 std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, QgsFields(), Qgis::WkbType::Polygon, mSource->sourceCrs() ) );
92 if ( !sink )
93 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
94
95#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
96 concaveHullQgis( sink, parameters, context, feedback );
97#else
98 concaveHullGeos( sink, parameters, feedback );
99#endif
100
101 QVariantMap outputs;
102 outputs.insert( QStringLiteral( "OUTPUT" ), dest );
103 return outputs;
104}
105
106void QgsConcaveHullAlgorithm::concaveHullGeos( std::unique_ptr< QgsFeatureSink > &sink, const QVariantMap &parameters, QgsProcessingFeedback *feedback )
107{
108 long long i = 0;
109
111 QgsFeature f;
112 QgsGeometry allPoints;
113 while ( it.nextFeature( f ) )
114 {
115 i++;
116 if ( feedback->isCanceled() )
117 return;
118
119 feedback->setProgress( i * mStep );
120
121 if ( !f.hasGeometry() )
122 continue;
123
124 const QgsAbstractGeometry *geom = f.geometry().constGet();
125 if ( QgsWkbTypes::isMultiType( geom->wkbType() ) )
126 {
127 const QgsMultiPoint mp( *qgsgeometry_cast< const QgsMultiPoint * >( geom ) );
128 for ( auto pit = mp.const_parts_begin(); pit != mp.const_parts_end(); ++pit )
129 {
130 allPoints.addPartV2( qgsgeometry_cast< QgsPoint * >( *pit )->clone(), Qgis::WkbType::Point );
131 }
132 }
133 else
134 {
135 allPoints.addPartV2( qgsgeometry_cast< QgsPoint * >( geom )->clone(), Qgis::WkbType::Point );
136 }
137 }
138 const QgsGeometry concaveHull = allPoints.concaveHull( mPercentage, mAllowHoles );
139
140 if ( mSplitMultipart && concaveHull.isMultipart() )
141 {
142 QVector< QgsGeometry > collection = concaveHull.asGeometryCollection();
143 mStep = collection.length() > 0 ? 50.0 / collection.length() : 1;
144 for ( int i = 0; i < collection.length(); i++ )
145 {
146 if ( feedback->isCanceled() )
147 {
148 break;
149 }
150
151 QgsGeometry geom = collection[i];
152 if ( !mAllowHoles )
153 {
154 geom = collection[i].removeInteriorRings();
155 }
156 QgsFeature f;
157 f.setGeometry( geom );
158 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
159 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
160
161 feedback->setProgress( 50 + i * mStep );
162 }
163 }
164 else
165 {
166 QgsGeometry geom( concaveHull );
167 if ( !mAllowHoles )
168 {
169 geom = concaveHull.removeInteriorRings();
170 }
171 QgsFeature f;
172 f.setGeometry( geom );
173 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
174 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
175 feedback->setProgress( 100 );
176 }
177}
178
179void QgsConcaveHullAlgorithm::concaveHullQgis( std::unique_ptr< QgsFeatureSink > &sink, const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
180{
181 QgsProcessingMultiStepFeedback multiStepFeedback( 5, feedback );
182
183 feedback->setProgressText( QObject::tr( "Creating Delaunay triangles…" ) );
184 multiStepFeedback.setCurrentStep( 1 );
185
186 QVariantMap params;
187 params["INPUT"] = parameters["INPUT"];
188 params["TOLERANCE"] = 0.0;
189 params["ADD_ATTRIBUTES"] = false;
190 params["OUTPUT"] = QgsProcessing::TEMPORARY_OUTPUT;
191 const QgsProcessingAlgorithm *delaunayAlg = QgsApplication::processingRegistry()->algorithmById( QStringLiteral( "native:delaunaytriangulation" ) );
192 if ( !delaunayAlg )
193 {
194 feedback->reportError( QObject::tr( "Failed to compute concave hull: Delaunay triangulation algorithm not found!" ), true );
195 }
196 std::unique_ptr< QgsProcessingAlgorithm > algorithm;
197 algorithm.reset( delaunayAlg->create() );
198 QVariantMap results = algorithm->run( params, context, &multiStepFeedback );
199 QgsVectorLayer *layer = qobject_cast<QgsVectorLayer *>( QgsProcessingUtils::mapLayerFromString( results["OUTPUT"].toString(), context ) );
200
201 if ( !layer )
202 {
203 throw QgsProcessingException( QObject::tr( "Failed to compute Delaunay triangulation." ) );
204 }
205
206 if ( layer->featureCount() == 0 )
207 {
208 throw QgsProcessingException( QObject::tr( "No Delaunay triangles created." ) );
209 }
210
211 feedback->setProgressText( QObject::tr( "Computing edges max length…" ) );
212 multiStepFeedback.setCurrentStep( 2 );
213
214 QVector<double> length;
215 QMap<long, double> edges;
216 long i = 0;
217 double step = layer->featureCount() > 0 ? 100.0 / layer->featureCount() : 1;
218 QgsFeatureIterator it = layer->getFeatures( QgsFeatureRequest().setNoAttributes() );
219 QgsFeature f;
220 while ( it.nextFeature( f ) )
221 {
222 i++;
223 if ( feedback->isCanceled() )
224 return;
225
226 multiStepFeedback.setProgress( i * step );
227
228 if ( !f.hasGeometry() )
229 continue;
230
231 QVector<QgsPointXY> points = f.geometry().asPolygon().at( 0 );
232 for ( int j = 0; j < points.size() - 1; j++ )
233 {
234 length << std::sqrt( points.at( j ).sqrDist( points.at( j + 1 ) ) );
235 }
236 QVector<double> vec = length.mid( length.size() - 3, -1 );
237 edges[f.id()] = *std::max_element( vec.constBegin(), vec.constEnd() );
238 }
239 const double maxLength = *std::max_element( length.constBegin(), length.constEnd() );
240
241 feedback->setProgressText( QObject::tr( "Removing features…" ) );
242 multiStepFeedback.setCurrentStep( 3 );
243 i = 0;
244 step = edges.size() > 0 ? 100.0 / edges.size() : 1;
245 QgsFeatureIds toDelete;
246 QMap<long, double>::iterator edgesIt = edges.begin();
247 while ( edgesIt != edges.end() )
248 {
249 if ( feedback->isCanceled() )
250 return;
251
252 if ( edgesIt.value() > mPercentage * maxLength )
253 {
254 toDelete << edgesIt.key();
255 }
256
257 ++edgesIt;
258 i++;
259 multiStepFeedback.setProgress( i * step );
260 }
261 layer->dataProvider()->deleteFeatures( toDelete );
262
263 feedback->setProgressText( QObject::tr( "Dissolving Delaunay triangles…" ) );
264 multiStepFeedback.setCurrentStep( 4 );
265 params.clear();
266 params["INPUT"] = layer->source();
267 params["OUTPUT"] = QgsProcessing::TEMPORARY_OUTPUT;
268 const QgsProcessingAlgorithm *dissolveAlg = QgsApplication::processingRegistry()->algorithmById( QStringLiteral( "native:dissolve" ) );
269 if ( !dissolveAlg )
270 {
271 feedback->reportError( QObject::tr( "Failed to compute concave hull: Dissolve algorithm not found!" ), true );
272 }
273 algorithm.reset( dissolveAlg->create() );
274 results = algorithm->run( params, context, &multiStepFeedback );
275 layer = qobject_cast<QgsVectorLayer *>( QgsProcessingUtils::mapLayerFromString( results["OUTPUT"].toString(), context ) );
276 if ( !layer )
277 {
278 throw QgsProcessingException( QObject::tr( "Failed to dissolve Delaunay triangles." ) );
279 }
280 if ( layer->featureCount() == 0 )
281 {
282 throw QgsProcessingException( QObject::tr( "There are no features in the dissolved layer." ) );
283 }
284
285 layer->getFeatures().nextFeature( f );
286 QgsGeometry concaveHull = f.geometry();
287
288 // save result
289 multiStepFeedback.setCurrentStep( 5 );
290
291 if ( mSplitMultipart && concaveHull.isMultipart() )
292 {
293 const QVector< QgsGeometry > collection = concaveHull.asGeometryCollection();
294 step = collection.length() > 0 ? 50.0 / collection.length() : 1;
295 for ( int i = 0; i < collection.length(); i++ )
296 {
297 if ( feedback->isCanceled() )
298 {
299 break;
300 }
301
302 QgsGeometry geom = collection[i];
303 if ( !mAllowHoles )
304 {
305 geom = collection[i].removeInteriorRings();
306 }
307 QgsFeature f;
308 f.setGeometry( geom );
309 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
310 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
311
312 multiStepFeedback.setProgress( i * step );
313 }
314 }
315 else
316 {
317 QgsGeometry geom( concaveHull );
318 if ( !mAllowHoles )
319 {
320 geom = concaveHull.removeInteriorRings();
321 }
322 QgsFeature f;
323 f.setGeometry( geom );
324 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
325 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
326 multiStepFeedback.setProgress( 100 );
327 }
328}
329
@ VectorPoint
Vector point layers.
@ VectorPolygon
Vector polygon layers.
@ SkipGeometryValidityChecks
Invalid geometry checks should always be skipped. This flag can be useful for algorithms which always...
@ Polygon
Polygon.
Abstract base class for all geometries.
Qgis::WkbType wkbType() const
Returns the WKB type of the geometry.
static QgsProcessingRegistry * processingRegistry()
Returns the application's processing registry, used for managing processing providers,...
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsFeatureId id
Definition qgsfeature.h:66
QgsGeometry geometry
Definition qgsfeature.h:69
bool hasGeometry() const
Returns true if the feature has an associated geometry.
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Container of fields for a vector layer.
Definition qgsfields.h:46
A geometry is the spatial representation of a feature.
QgsPolygonXY asPolygon() const
Returns the contents of the geometry as a polygon.
QVector< QgsGeometry > asGeometryCollection() const
Returns contents of the geometry as a list of geometries.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
QgsGeometry concaveHull(double targetPercent, bool allowHoles=false) const
Returns a possibly concave polygon that contains all the points in the geometry.
Qgis::GeometryOperationResult addPartV2(const QVector< QgsPointXY > &points, Qgis::WkbType wkbType=Qgis::WkbType::Unknown)
Adds a new part to a the geometry.
QgsGeometry removeInteriorRings(double minimumAllowedArea=-1) const
Removes the interior rings from a (multi)polygon geometry.
QString source() const
Returns the source for the layer.
Multi point geometry collection.
Abstract base class for processing algorithms.
QVariantMap run(const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback, bool *ok=nullptr, const QVariantMap &configuration=QVariantMap(), bool catchExceptions=true) const
Executes the algorithm using the specified parameters.
QgsProcessingAlgorithm * create(const QVariantMap &configuration=QVariantMap()) const
Creates a copy of the algorithm, ready for execution.
Contains information about the context in which a processing algorithm is executed.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
virtual void reportError(const QString &error, bool fatalError=false)
Reports that the algorithm encountered an error while executing.
virtual void setProgressText(const QString &text)
Sets a progress report text string.
Processing feedback object for multi-step operations.
A boolean parameter for processing algorithms.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A numeric parameter for processing algorithms.
const QgsProcessingAlgorithm * algorithmById(const QString &id) const
Finds an algorithm by its ID.
static QgsMapLayer * mapLayerFromString(const QString &string, QgsProcessingContext &context, bool allowLoadingNewLayers=true, QgsProcessingUtils::LayerHint typeHint=QgsProcessingUtils::LayerHint::UnknownType, QgsProcessing::LayerOptionsFlags flags=QgsProcessing::LayerOptionsFlags())
Interprets a string as a map layer within the supplied context.
static const QString TEMPORARY_OUTPUT
Constant used to indicate that a Processing algorithm output should be a temporary layer/file.
virtual bool deleteFeatures(const QgsFeatureIds &id)
Deletes one or more features from the provider.
Represents a vector layer which manages a vector based data sets.
long long featureCount(const QString &legendKey) const
Number of features rendered with specified legend key.
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const FINAL
Queries the layer for features specified in request.
QgsVectorDataProvider * dataProvider() FINAL
Returns the layer's data provider, it may be nullptr.
static bool isMultiType(Qgis::WkbType type)
Returns true if the WKB type is a multi type.
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into allowing algorithms to be written in pure substantial changes are required in order to port existing x Processing algorithms for QGIS x The most significant changes are outlined not GeoAlgorithm For algorithms which operate on features one by consider subclassing the QgsProcessingFeatureBasedAlgorithm class This class allows much of the boilerplate code for looping over features from a vector layer to be bypassed and instead requires implementation of a processFeature method Ensure that your algorithm(or algorithm 's parent class) implements the new pure virtual createInstance(self) call
QSet< QgsFeatureId > QgsFeatureIds