processing: Add approximate medial axis processing

SFCGAL backend is needed.

src/analysis/CMakeLists.txt

@@ -44,6 +44,7 @@ set(QGIS_ANALYSIS_SRCS
   processing/qgsalgorithmangletonearest.cpp
   processing/qgsalgorithmannotations.cpp
   processing/qgsalgorithmapplylayerstyle.cpp
+  processing/qgsalgorithmapproximatemedialaxis.cpp
   processing/qgsalgorithmarraytranslatedfeatures.cpp
   processing/qgsalgorithmaspect.cpp
   processing/qgsalgorithmassignprojection.cpp
@@ -632,6 +633,10 @@ if (WITH_DRACO)
   target_link_libraries(qgis_analysis ${DRACO_LIBRARY})
 endif()
 
+if (WITH_SFCGAL)
+  target_link_libraries(qgis_analysis SFCGAL::SFCGAL)
+endif()
+
 target_compile_definitions(qgis_analysis PRIVATE "-DQT_NO_FOREACH")
 
 # clang-tidy

src/analysis/processing/qgsalgorithmapproximatemedialaxis.cpp

@@ -0,0 +1,167 @@
+/***************************************************************************
+                         qgsalgorithmapproximatemedialaxis.cpp
+                         ---------------------
+    begin                : September 2025
+    copyright            : (C) 2025 by Jean Felder
+    email                : jean dot felder at oslandia dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+
+#include "qgsalgorithmapproximatemedialaxis.h"
+#include "qgsexception.h"
+#ifdef WITH_SFCGAL
+#include "qgssfcgalgeometry.h"
+#endif
+
+///@cond PRIVATE
+
+QString QgsApproximateMedialAxisAlgorithm::name() const
+{
+  return QStringLiteral( "approximatemedialaxis" );
+}
+
+QString QgsApproximateMedialAxisAlgorithm::displayName() const
+{
+  return QObject::tr( "Approximate medial axis" );
+}
+
+QStringList QgsApproximateMedialAxisAlgorithm::tags() const
+{
+  return QObject::tr( "medial,axis,create,lines,straight,skeleton" ).split( ',' );
+}
+
+QString QgsApproximateMedialAxisAlgorithm::group() const
+{
+  return QObject::tr( "Vector geometry" );
+}
+
+QString QgsApproximateMedialAxisAlgorithm::groupId() const
+{
+  return QStringLiteral( "vectorgeometry" );
+}
+
+QString QgsApproximateMedialAxisAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "The Approximate Medial Axis algorithm generates a simplified skeleton of a shape by approximating its medial axis. \n\n"
+                      "The output is a collection of lines that follow the central structure of the shape. The result is a thin, stable set "
+                      "of curves that capture the main topology while ignoring noise.\n\n"
+                      "This algorithm ignores the Z dimensions. If the geometry is 3D, the approximate medial axis will be calculated from "
+                      "its 2D projection." );
+}
+
+QString QgsApproximateMedialAxisAlgorithm::shortDescription() const
+{
+  return QObject::tr( "Returns an approximate medial axis for a polygon layer input based on its straight skeleton." );
+}
+
+QgsApproximateMedialAxisAlgorithm *QgsApproximateMedialAxisAlgorithm::createInstance() const
+{
+  return new QgsApproximateMedialAxisAlgorithm();
+}
+
+QgsFields QgsApproximateMedialAxisAlgorithm::outputFields( const QgsFields &inputFields ) const
+{
+  QgsFields newFields;
+  newFields.append( QgsField( QStringLiteral( "length" ), QMetaType::Type::Double, QString(), 20, 6 ) );
+  return QgsProcessingUtils::combineFields( inputFields, newFields );
+}
+
+QList<int> QgsApproximateMedialAxisAlgorithm::inputLayerTypes() const
+{
+  return QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorPolygon );
+}
+
+QString QgsApproximateMedialAxisAlgorithm::outputName() const
+{
+  return QObject::tr( "Medial axis" );
+}
+
+Qgis::ProcessingSourceType QgsApproximateMedialAxisAlgorithm::outputLayerType() const
+{
+  return Qgis::ProcessingSourceType::VectorLine;
+}
+
+Qgis::WkbType QgsApproximateMedialAxisAlgorithm::outputWkbType( Qgis::WkbType inputWkbType ) const
+{
+  Q_UNUSED( inputWkbType )
+  return Qgis::WkbType::MultiLineString;
+}
+
+bool QgsApproximateMedialAxisAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  Q_UNUSED( parameters )
+  Q_UNUSED( context )
+  Q_UNUSED( feedback )
+
+#ifdef WITH_SFCGAL
+  return true;
+#else
+  throw QgsProcessingException( QObject::tr( "This processing requires a QGIS installation with SFCGAL support enabled. Please use a version of QGIS that includes SFCGAL." ) );
+#endif
+}
+
+QgsFeatureList QgsApproximateMedialAxisAlgorithm::processFeature( const QgsFeature &feature, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  Q_UNUSED( context )
+
+#ifdef WITH_SFCGAL
+  QgsFeature modifiedFeature = feature;
+  if ( modifiedFeature.hasGeometry() )
+  {
+    QgsGeometry outputGeometry;
+    QgsSfcgalGeometry inputSfcgalGeometry;
+    try
+    {
+      inputSfcgalGeometry = QgsSfcgalGeometry( modifiedFeature.geometry() );
+    }
+    catch ( const QgsSfcgalException &exception )
+    {
+      feedback->reportError( QObject::tr( "Cannot calculate approximate medial axis for feature %1: %2" ).arg( feature.id() ).arg( exception.what() ) );
+      modifiedFeature.clearGeometry();
+    }
+
+    if ( !inputSfcgalGeometry.isEmpty() )
+    {
+      try
+      {
+        std::unique_ptr<QgsSfcgalGeometry> outputSfcgalGeometry = inputSfcgalGeometry.approximateMedialAxis();
+        outputGeometry = QgsGeometry( outputSfcgalGeometry->asQgisGeometry() );
+        modifiedFeature.setGeometry( outputGeometry );
+      }
+      catch ( const QgsSfcgalException &medialAxisException )
+      {
+        feedback->reportError( QObject::tr( "Cannot calculate approximate medial axis for feature %1: %2" ).arg( feature.id() ).arg( medialAxisException.what() ) );
+      }
+    }
+
+    if ( !outputGeometry.isNull() )
+    {
+      QgsAttributes attrs = modifiedFeature.attributes();
+      attrs << outputGeometry.constGet()->length();
+      modifiedFeature.setAttributes( attrs );
+    }
+    else
+    {
+      QgsAttributes attrs = modifiedFeature.attributes();
+      attrs << QVariant();
+      modifiedFeature.setAttributes( attrs );
+    }
+  }
+  return QgsFeatureList() << modifiedFeature;
+#else
+  Q_UNUSED( feature )
+  Q_UNUSED( feedback )
+  throw QgsProcessingException( QObject::tr( "This processing requires a QGIS installation with SFCGAL support enabled. Please use a version of QGIS that includes SFCGAL." ) );
+#endif
+}
+
+///@endcond PRIVATE

src/analysis/processing/qgsalgorithmapproximatemedialaxis.h

@@ -0,0 +1,56 @@
+/***************************************************************************
+                         qgsalgorithmapproximatemedialaxis.h
+                         ---------------------
+    begin                : September 2025
+    copyright            : (C) 2025 by Jean Felder
+    email                : jean dot felder at oslandia dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSALGORITHMAPPROXIMATEMEDIALAXIS_H
+#define QGSALGORITHMAPPROXIMATEMEDIALAXIS_H
+
+#define SIP_NO_FILE
+
+#include "qgis_sip.h"
+#include "qgsprocessingalgorithm.h"
+
+///@cond PRIVATE
+
+
+/**
+ * Approximate medial axis algorithm with SFCGAL backend.
+ */
+class QgsApproximateMedialAxisAlgorithm : public QgsProcessingFeatureBasedAlgorithm
+{
+  public:
+    QgsApproximateMedialAxisAlgorithm() = default;
+    QString name() const override;
+    QString displayName() const override;
+    QStringList tags() const override;
+    QString group() const override;
+    QString groupId() const override;
+    QString shortHelpString() const override;
+    QString shortDescription() const override;
+    QgsApproximateMedialAxisAlgorithm *createInstance() const override SIP_FACTORY;
+    QList<int> inputLayerTypes() const override;
+
+  protected:
+    QString outputName() const override;
+    Qgis::ProcessingSourceType outputLayerType() const override;
+    Qgis::WkbType outputWkbType( Qgis::WkbType inputWkbType ) const override;
+    QgsFields outputFields( const QgsFields &inputFields ) const override;
+    bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+    QgsFeatureList processFeature( const QgsFeature &feature, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+};
+
+///@endcond PRIVATE
+#endif // QGSALGORITHMAPPROXIMATEMEDIALAXIS_H

src/analysis/processing/qgsnativealgorithms.cpp

@@ -28,6 +28,7 @@
 #include "qgsalgorithmangletonearest.h"
 #include "qgsalgorithmannotations.h"
 #include "qgsalgorithmapplylayerstyle.h"
+#include "qgsalgorithmapproximatemedialaxis.h"
 #include "qgsalgorithmarraytranslatedfeatures.h"
 #include "qgsalgorithmaspect.h"
 #include "qgsalgorithmassignprojection.h"
@@ -351,6 +352,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
   addAlgorithm( new QgsAngleToNearestAlgorithm() );
   addAlgorithm( new QgsApplyLayerMetadataAlgorithm() );
   addAlgorithm( new QgsApplyLayerStyleAlgorithm() );
+  addAlgorithm( new QgsApproximateMedialAxisAlgorithm() );
   addAlgorithm( new QgsArrayTranslatedFeaturesAlgorithm() );
   addAlgorithm( new QgsAspectAlgorithm() );
   addAlgorithm( new QgsAssignProjectionAlgorithm() );

tests/src/analysis/CMakeLists.txt

@@ -43,6 +43,10 @@ if (WITH_PDAL AND PDAL_2_5_OR_HIGHER)
   set(TESTS ${TESTS} testqgsprocessingpdalalgs.cpp)
 endif()
 
+if (WITH_SFCGAL)
+  set(TESTS ${TESTS} testqgsprocessingsfcgalalgs.cpp)
+endif()
+
 foreach(TESTSRC ${TESTS})
     add_qgis_test(${TESTSRC} MODULE analysis LINKEDLIBRARIES qgis_analysis)
 endforeach(TESTSRC)

tests/src/analysis/testqgsprocessingsfcgalalgs.cpp

@@ -0,0 +1,157 @@
+/***************************************************************************
+                         testqgsprocessingsfcgalalgs.cpp
+                         ---------------------
+    begin                : September 2025
+    copyright            : (C) 2025 by Jean Felder
+    email                : jean dot felder at oslandia dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "qgsabstractgeometry.h"
+#include "qgsgeometry.h"
+#include "qgslinestring.h"
+#include "qgsmultilinestring.h"
+#include "qgstest.h"
+#include "qgsnativealgorithms.h"
+#include "qgsprocessingregistry.h"
+#include "qgsprocessingcontext.h"
+#include "qgsvectorlayer.h"
+#include "qgssfcgalgeometry.h"
+#include <qtestcase.h>
+#include <memory>
+
+// #include <QThread>
+
+class TestQgsProcessingSfcgalAlgs : public QgsTest
+{
+    Q_OBJECT
+
+  public:
+    TestQgsProcessingSfcgalAlgs()
+      : QgsTest( QStringLiteral( "Processing SFCGAL Algorithms Test" ) )
+    {}
+
+  private slots:
+    void initTestCase();    // will be called before the first testfunction is executed.
+    void cleanupTestCase(); // will be called after the last testfunction was executed.
+    void init();            // will be called before each testfunction is executed.
+    void cleanup();         // will be called after every testfunction.
+
+    void medialAxis();
+
+  private:
+    QgsGeometry openWktFile( const QString &wktFile );
+    QgsVectorLayer *mPolygonLayer = nullptr;
+};
+
+void TestQgsProcessingSfcgalAlgs::initTestCase()
+{
+  QgsApplication::init();
+  QgsApplication::initQgis();
+
+  // Set up the QgsSettings environment
+  QCoreApplication::setOrganizationName( QStringLiteral( "QGIS" ) );
+  QCoreApplication::setOrganizationDomain( QStringLiteral( "qgis.org" ) );
+  QCoreApplication::setApplicationName( QStringLiteral( "QGIS-TEST" ) );
+
+  QgsApplication::processingRegistry()->addProvider( new QgsNativeAlgorithms( QgsApplication::processingRegistry() ) );
+
+  const QString dataDir( TEST_DATA_DIR ); //defined in CmakeLists.txt
+
+  const QString polysFileName = dataDir + "/polys.shp";
+  const QFileInfo polyFileInfo( polysFileName );
+  mPolygonLayer = new QgsVectorLayer( polyFileInfo.filePath(), QStringLiteral( "polygons" ), QStringLiteral( "ogr" ) );
+  QVERIFY( mPolygonLayer->isValid() );
+
+  // Register the layer with the registry
+  QgsProject::instance()->addMapLayer( mPolygonLayer );
+}
+
+void TestQgsProcessingSfcgalAlgs::cleanupTestCase()
+{
+  QgsProject::instance()->removeMapLayer( mPolygonLayer );
+  QgsApplication::exitQgis();
+}
+
+void TestQgsProcessingSfcgalAlgs::init()
+{
+}
+
+void TestQgsProcessingSfcgalAlgs::cleanup()
+{
+}
+
+QgsGeometry TestQgsProcessingSfcgalAlgs::openWktFile( const QString &wktFile )
+{
+  QString expectedPath = testDataPath( QString( "control_files/expected_sfcgal/expected_%1" ).arg( wktFile ) );
+  QFile expectedFile( expectedPath );
+  if ( !expectedFile.open( QFile::ReadOnly | QIODevice::Text ) )
+  {
+    qWarning() << "Unable to open expected data file" << expectedPath;
+    return QgsGeometry();
+  }
+
+  // remove '\n' from dumped file
+  QByteArray expectedBA = expectedFile.readAll();
+  QString expectedStr;
+  for ( int i = 0; i < expectedBA.length(); i++ )
+  {
+    if ( expectedBA.at( i ) != '\n' )
+      expectedStr += expectedBA.at( i );
+  }
+
+  // load geom from corrected wkt
+  return QgsGeometry::fromWkt( expectedStr );
+}
+
+void TestQgsProcessingSfcgalAlgs::medialAxis()
+{
+  std::unique_ptr<QgsProcessingAlgorithm> alg( QgsApplication::processingRegistry()->createAlgorithmById( QStringLiteral( "native:approximatemedialaxis" ) ) );
+  QVERIFY( alg != nullptr );
+
+  auto context = std::make_unique<QgsProcessingContext>();
+  context->setProject( QgsProject::instance() );
+
+  QVariantMap parameters;
+  parameters.insert( QStringLiteral( "INPUT" ), QVariant::fromValue( mPolygonLayer ) );
+  parameters.insert( QStringLiteral( "OUTPUT" ), QgsProcessing::TEMPORARY_OUTPUT );
+
+  bool ok = false;
+  QgsProcessingFeedback feedback;
+  QVariantMap results = alg->run( parameters, *context, &feedback, &ok );
+  QVERIFY( ok );
+
+  QgsVectorLayer *outputLayer = qobject_cast<QgsVectorLayer *>( context->getMapLayer( results.value( QStringLiteral( "OUTPUT" ) ).toString() ) );
+  QVERIFY( outputLayer );
+  QVERIFY( outputLayer->isValid() );
+  QCOMPARE( outputLayer->featureCount(), 10 );
+
+  // retrieve first feature
+  QgsFeatureIterator outputFeatureIterator = outputLayer->dataProvider()->getFeatures();
+  QgsFeature outputFeature;
+  QVERIFY( outputFeatureIterator.nextFeature( outputFeature ) );
+
+  // check attribute
+  QGSCOMPARENEAR( outputFeature.attribute( "length" ).toDouble(), 14.01481, 0.00001 );
+
+  // check geometry
+  const QgsGeometry outputGeom( outputFeature.geometry() );
+  QVERIFY( outputGeom.wkbType() == Qgis::WkbType::MultiLineString );
+  const QgsGeometry expectedGeom = openWktFile( "processing_medial_axis.wkt" );
+  QVERIFY( expectedGeom.wkbType() == Qgis::WkbType::MultiLineString );
+
+  QCOMPARE( outputGeom.asWkt( 4 ), expectedGeom.asWkt( 4 ) );
+
+  QgsProject::instance()->removeMapLayer( outputLayer );
+}
+
+QGSTEST_MAIN( TestQgsProcessingSfcgalAlgs )
+#include "testqgsprocessingsfcgalalgs.moc"