use geos to calculate concave hull on systems with geos >=3.11.0

2025-12-07 00:03:52 -05:00 · 2023-08-09 13:56:46 +03:00 · 2023-08-09 13:56:46 +03:00 · 514e8b5b24
commit 514e8b5b24
parent 496d8267cc
5 changed files with 197 additions and 42 deletions
--- a/python/plugins/processing/tests/testdata/expected/concave_hull_points_09.gml
+++ b/python/plugins/processing/tests/testdata/expected/concave_hull_points_09.gml
@ -0,0 +1,16 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<ogr:FeatureCollection
+     gml:id="aFeatureCollection"
+     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+     xsi:schemaLocation="http://ogr.maptools.org/ concave_hull_points_09.xsd"
+     xmlns:ogr="http://ogr.maptools.org/"
+     xmlns:gml="http://www.opengis.net/gml/3.2">
+  <gml:boundedBy><gml:Envelope srsName="urn:ogc:def:crs:EPSG::4326"><gml:lowerCorner>-5 0</gml:lowerCorner><gml:upperCorner>3 8</gml:upperCorner></gml:Envelope></gml:boundedBy>
+                                                                                                                                                                             
+  <ogr:featureMember>
+    <ogr:concave_hull_points_09 gml:id="concave_hull_points_09.0">
+      <gml:boundedBy><gml:Envelope srsName="urn:ogc:def:crs:EPSG::4326"><gml:lowerCorner>-5 0</gml:lowerCorner><gml:upperCorner>3 8</gml:upperCorner></gml:Envelope></gml:boundedBy>
+      <ogr:geometryProperty><gml:Polygon srsName="urn:ogc:def:crs:EPSG::4326" gml:id="concave_hull_points_09.geom.0"><gml:exterior><gml:LinearRing><gml:posList>-1 0 1 1 2 2 3 3 2 5 -1 8 -1 7 -5 0 -1 0</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon></ogr:geometryProperty>
+    </ogr:concave_hull_points_09>
+  </ogr:featureMember>
+</ogr:FeatureCollection>
--- a/python/plugins/processing/tests/testdata/expected/concave_hull_points_09.xsd
+++ b/python/plugins/processing/tests/testdata/expected/concave_hull_points_09.xsd
@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<xs:schema 
+    targetNamespace="http://ogr.maptools.org/"
+    xmlns:ogr="http://ogr.maptools.org/"
+    xmlns:xs="http://www.w3.org/2001/XMLSchema"
+    xmlns:gml="http://www.opengis.net/gml/3.2"
+    xmlns:gmlsf="http://www.opengis.net/gmlsf/2.0"
+    elementFormDefault="qualified"
+    version="1.0">
+<xs:annotation>
+  <xs:appinfo source="http://schemas.opengis.net/gmlsfProfile/2.0/gmlsfLevels.xsd">
+    <gmlsf:ComplianceLevel>0</gmlsf:ComplianceLevel>
+  </xs:appinfo>
+</xs:annotation>
+<xs:import namespace="http://www.opengis.net/gml/3.2" schemaLocation="http://schemas.opengis.net/gml/3.2.1/gml.xsd"/>
+<xs:import namespace="http://www.opengis.net/gmlsf/2.0" schemaLocation="http://schemas.opengis.net/gmlsfProfile/2.0/gmlsfLevels.xsd"/>
+<xs:element name="FeatureCollection" type="ogr:FeatureCollectionType" substitutionGroup="gml:AbstractFeature"/>
+<xs:complexType name="FeatureCollectionType">
+  <xs:complexContent>
+    <xs:extension base="gml:AbstractFeatureType">
+      <xs:sequence minOccurs="0" maxOccurs="unbounded">
+        <xs:element name="featureMember">
+          <xs:complexType>
+            <xs:complexContent>
+              <xs:extension base="gml:AbstractFeatureMemberType">
+                <xs:sequence>
+                  <xs:element ref="gml:AbstractFeature"/>
+                </xs:sequence>
+              </xs:extension>
+            </xs:complexContent>
+          </xs:complexType>
+        </xs:element>
+      </xs:sequence>
+    </xs:extension>
+  </xs:complexContent>
+</xs:complexType>
+<xs:element name="concave_hull_points_09" type="ogr:concave_hull_points_09_Type" substitutionGroup="gml:AbstractFeature"/>
+<xs:complexType name="concave_hull_points_09_Type">
+  <xs:complexContent>
+    <xs:extension base="gml:AbstractFeatureType">
+      <xs:sequence>
+        <xs:element name="geometryProperty" type="gml:SurfacePropertyType" nillable="true" minOccurs="0" maxOccurs="1"/> <!-- restricted to Polygon --><!-- srsName="urn:ogc:def:crs:EPSG::4326" -->
+      </xs:sequence>
+    </xs:extension>
+  </xs:complexContent>
+</xs:complexType>
+</xs:schema>
--- a/python/plugins/processing/tests/testdata/qgis_algorithm_tests4.yaml
+++ b/python/plugins/processing/tests/testdata/qgis_algorithm_tests4.yaml
@ -268,23 +268,6 @@ tests:
              cast: str
            fid: skip

-  - algorithm: native:concavehull
-    name: Concave Hull - Points (0.3)
-    params:
-      ALPHA: 0.3
-      HOLES: true
-      INPUT:
-        name: points.gml
-        type: vector
-      NO_MULTIGEOMETRY: false
-    results:
-      OUTPUT:
-        name: expected/concave_hull_points_03.gml
-        type: vector
-        compare:
-          geometry:
-            normalize: True
-
  - algorithm: native:concavehull
    name: Concave Hull - Points (0.7)
    params:
@ -302,6 +285,23 @@ tests:
          geometry:
            normalize: True

+  - algorithm: native:concavehull
+    name: Concave Hull - Points (0.9)
+    params:
+      ALPHA: 1
+      HOLES: true
+      INPUT:
+        name: points.gml
+        type: vector
+      NO_MULTIGEOMETRY: false
+    results:
+      OUTPUT:
+        name: expected/concave_hull_points_09.gml
+        type: vector
+        compare:
+          geometry:
+            normalize: True
+
  - algorithm: qgis:knearestconcavehull
    name: K-nearest Neighbor Concave Hull - Points (k=7)
    params:
--- a/src/analysis/processing/qgsalgorithmconcavehull.cpp
+++ b/src/analysis/processing/qgsalgorithmconcavehull.cpp
@ -67,27 +67,117 @@ void QgsConcaveHullAlgorithm::initAlgorithm( const QVariantMap & )
  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Concave hull" ), QgsProcessing::TypeVectorPolygon ) );
 }

-QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+bool QgsConcaveHullAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
 {
-  std::unique_ptr< QgsProcessingFeatureSource > source( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
-  if ( !source )
+  mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
+  if ( !mSource )
    throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );

-  if ( source->featureCount() < 3 )
+  if ( mSource->featureCount() < 3 )
    throw QgsProcessingException( QObject::tr( "Input layer should contain at least 3 points." ) );

-  const double percentage = parameterAsDouble( parameters, QStringLiteral( "ALPHA" ), context );
-  const bool allowHoles = parameterAsBool( parameters, QStringLiteral( "HOLES" ), context );
-  const bool splitMultipart = parameterAsBool( parameters, QStringLiteral( "NO_MULTIGEOMETRY" ), context );
+  mPercentage = parameterAsDouble( parameters, QStringLiteral( "ALPHA" ), context );
+  mAllowHoles = parameterAsBool( parameters, QStringLiteral( "HOLES" ), context );
+  mSplitMultipart = parameterAsBool( parameters, QStringLiteral( "NO_MULTIGEOMETRY" ), context );

-  QgsFields fields;
-  fields.append( QgsField( QStringLiteral( "id" ), QVariant::LongLong ) );
+  return true;
+}

+QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
  QString dest;
-  std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, Qgis::WkbType::Polygon, source->sourceCrs() ) );
+  std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, QgsFields(), Qgis::WkbType::Polygon, mSource->sourceCrs() ) );
  if ( !sink )
    throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );

+  QgsGeometry concaveHull;
+#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
+  concaveHullQgis( sink, parameters, context, feedback );
+#else
+  concaveHullGeos( sink, parameters, feedback );
+#endif
+
+  QVariantMap outputs;
+  outputs.insert( QStringLiteral( "OUTPUT" ), dest );
+  return outputs;
+}
+
+void QgsConcaveHullAlgorithm::concaveHullGeos( std::unique_ptr< QgsFeatureSink > &sink, const QVariantMap &parameters, QgsProcessingFeedback *feedback )
+{
+  long i = 0;
+  double step = mSource->featureCount() > 0 ? 50.0 / mSource->featureCount() : 1;
+
+  QgsFeatureIterator it = mSource->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QList< int >() ), QgsProcessingFeatureSource::FlagSkipGeometryValidityChecks );
+  QgsFeature f;
+  QgsGeometry allPoints;
+  while ( it.nextFeature( f ) )
+  {
+    i++;
+    if ( feedback->isCanceled() )
+      return;
+
+    feedback->setProgress( i * step );
+
+    if ( !f.hasGeometry() )
+      continue;
+
+    const QgsAbstractGeometry *geom = f.geometry().constGet();
+    if ( QgsWkbTypes::isMultiType( geom->wkbType() ) )
+    {
+      const QgsMultiPoint mp( *qgsgeometry_cast< const QgsMultiPoint * >( geom ) );
+      for ( auto pit = mp.const_parts_begin(); pit != mp.const_parts_end(); ++pit )
+      {
+        allPoints.addPart( qgsgeometry_cast< QgsPoint * >( *pit )->clone(), Qgis::GeometryType::Point );
+      }
+    }
+    else
+    {
+      allPoints.addPart( qgsgeometry_cast< QgsPoint * >( geom )->clone(), Qgis::GeometryType::Point );
+    }
+  }
+  QgsGeometry concaveHull = allPoints.concaveHull( mPercentage, mAllowHoles );
+
+  if ( mSplitMultipart && concaveHull.isMultipart() )
+  {
+    QVector< QgsGeometry > collection = concaveHull.asGeometryCollection();
+    step = collection.length() > 0 ? 50.0 / collection.length() : 1;
+    for ( int i = 0; i < collection.length(); i++ )
+    {
+      if ( feedback->isCanceled() )
+      {
+        break;
+      }
+
+      QgsGeometry geom = collection[i];
+      if ( !mAllowHoles )
+      {
+        geom = collection[i].removeInteriorRings();
+      }
+      QgsFeature f;
+      f.setGeometry( geom );
+      if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
+        throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
+
+      feedback->setProgress( 50 + i * step );
+    }
+  }
+  else
+  {
+    QgsGeometry geom( concaveHull );
+    if ( !mAllowHoles )
+    {
+      geom = concaveHull.removeInteriorRings();
+    }
+    QgsFeature f;
+    f.setGeometry( geom );
+    if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
+      throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
+    feedback->setProgress( 100 );
+  }
+}
+
+void QgsConcaveHullAlgorithm::concaveHullQgis( std::unique_ptr< QgsFeatureSink > &sink, const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
  QgsProcessingMultiStepFeedback multiStepFeedback( 5, feedback );

  feedback->setProgressText( QObject::tr( "Creating Delaunay triangles…" ) );
@ -132,7 +222,7 @@ QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parame
  {
    i++;
    if ( feedback->isCanceled() )
-      break;
+      return;

    multiStepFeedback.setProgress( i * step );

@ -158,14 +248,14 @@ QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parame
  while ( edgesIt != edges.end() )
  {
    if ( feedback->isCanceled() )
-      break;
+      return;

-    if ( edgesIt.value() > percentage * maxLength )
+    if ( edgesIt.value() > mPercentage * maxLength )
    {
      toDelete << edgesIt.key();
    }

-    ++edgesIt;
+    edgesIt++;
    i++;
    multiStepFeedback.setProgress( i * step );
  }
@ -200,7 +290,7 @@ QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parame
  // save result
  multiStepFeedback.setCurrentStep( 5 );

-  if ( splitMultipart && concaveHull.isMultipart() )
+  if ( mSplitMultipart && concaveHull.isMultipart() )
  {
    QVector< QgsGeometry > collection = concaveHull.asGeometryCollection();
    step = collection.length() > 0 ? 50.0 / collection.length() : 1;
@ -212,13 +302,11 @@ QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parame
      }

      QgsGeometry geom = collection[i];
-      if ( !allowHoles )
+      if ( !mAllowHoles )
      {
        geom = collection[i].removeInteriorRings();
      }
      QgsFeature f;
-      f.setFields( fields );
-      f.setAttributes( QgsAttributes() << i );
      f.setGeometry( geom );
      if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
        throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
@ -229,22 +317,16 @@ QVariantMap QgsConcaveHullAlgorithm::processAlgorithm( const QVariantMap &parame
  else
  {
    QgsGeometry geom( concaveHull );
-    if ( !allowHoles )
+    if ( !mAllowHoles )
    {
      geom = concaveHull.removeInteriorRings();
    }
    QgsFeature f;
-    f.setFields( fields );
-    f.setAttributes( QgsAttributes() << 1 );
    f.setGeometry( geom );
    if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
      throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
    multiStepFeedback.setProgress( 100 );
  }
-
-  QVariantMap outputs;
-  outputs.insert( QStringLiteral( "OUTPUT" ), dest );
-  return outputs;
 }

 ///@endcond
--- a/src/analysis/processing/qgsalgorithmconcavehull.h
+++ b/src/analysis/processing/qgsalgorithmconcavehull.h
@ -47,8 +47,18 @@ class QgsConcaveHullAlgorithm : public QgsProcessingAlgorithm

  protected:

+    bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
    QVariantMap processAlgorithm( const QVariantMap &parameters,
                                  QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+
+  private:
+    void concaveHullGeos( std::unique_ptr< QgsFeatureSink > &sink, const QVariantMap &parameters, QgsProcessingFeedback *feedback );
+    void concaveHullQgis( std::unique_ptr< QgsFeatureSink > &sink, const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback );
+
+    std::unique_ptr< QgsProcessingFeatureSource > mSource;
+    double mPercentage;
+    bool mAllowHoles;
+    bool mSplitMultipart;
 };

 ///@endcond PRIVATE