qgsalgorithmdbscanclustering.cpp

Go to the documentation of this file.

/***************************************************************************
                         qgsalgorithmdbscanclustering.cpp
                         ---------------------
    begin                : July 2018
    copyright            : (C) 2018 by Nyall Dawson
    email                : nyall dot dawson at gmail 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 "qgsalgorithmdbscanclustering.h"
 
#include <unordered_set>
 
#include "qgsspatialindexkdbush.h"
 
#include <QString>
 
using namespace Qt::StringLiterals;

 
QString QgsDbscanClusteringAlgorithm::name() const
{
  return u"dbscanclustering"_s;
}
 
QString QgsDbscanClusteringAlgorithm::displayName() const
{
  return QObject::tr( "DBSCAN clustering" );
}
 
QString QgsDbscanClusteringAlgorithm::shortDescription() const
{
  return QObject::tr( "Clusters point features using a density based scan algorithm." );
}
 
QStringList QgsDbscanClusteringAlgorithm::tags() const
{
  return QObject::tr( "clustering,clusters,density,based,points,distance" ).split( ',' );
}
 
QString QgsDbscanClusteringAlgorithm::group() const
{
  return QObject::tr( "Vector analysis" );
}
 
QString QgsDbscanClusteringAlgorithm::groupId() const
{
  return u"vectoranalysis"_s;
}
 
void QgsDbscanClusteringAlgorithm::initAlgorithm( const QVariantMap & )
{
  addParameter( new QgsProcessingParameterFeatureSource( u"INPUT"_s, QObject::tr( "Input layer" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorPoint ) ) );
  addParameter( new QgsProcessingParameterNumber( u"MIN_SIZE"_s, QObject::tr( "Minimum cluster size" ), Qgis::ProcessingNumberParameterType::Integer, 5, false, 1 ) );
  addParameter( new QgsProcessingParameterDistance( u"EPS"_s, QObject::tr( "Maximum distance between clustered points" ), 1, u"INPUT"_s, false, 0 ) );
 
  auto dbscanStarParam = std::make_unique<QgsProcessingParameterBoolean>( u"DBSCAN*"_s, QObject::tr( "Treat border points as noise (DBSCAN*)" ), false );
  dbscanStarParam->setFlags( dbscanStarParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
  addParameter( dbscanStarParam.release() );
 
  auto fieldNameParam = std::make_unique<QgsProcessingParameterString>( u"FIELD_NAME"_s, QObject::tr( "Cluster field name" ), u"CLUSTER_ID"_s );
  fieldNameParam->setFlags( fieldNameParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
  addParameter( fieldNameParam.release() );
  auto sizeFieldNameParam = std::make_unique<QgsProcessingParameterString>( u"SIZE_FIELD_NAME"_s, QObject::tr( "Cluster size field name" ), u"CLUSTER_SIZE"_s );
  sizeFieldNameParam->setFlags( sizeFieldNameParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
  addParameter( sizeFieldNameParam.release() );
 
  addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT"_s, QObject::tr( "Clusters" ), Qgis::ProcessingSourceType::VectorPoint ) );
 
  addOutput( new QgsProcessingOutputNumber( u"NUM_CLUSTERS"_s, QObject::tr( "Number of clusters" ) ) );
}
 
QString QgsDbscanClusteringAlgorithm::shortHelpString() const
{
  return QObject::tr(
    "This algorithm clusters point features based on a 2D implementation of Density-based spatial clustering of applications with noise (DBSCAN) algorithm.\n\n"
    "The algorithm requires two parameters, a minimum cluster size (“minPts”), and the maximum distance allowed between clustered points (“eps”)."
  );
}
 
QgsDbscanClusteringAlgorithm *QgsDbscanClusteringAlgorithm::createInstance() const
{
  return new QgsDbscanClusteringAlgorithm();
}
 
struct KDBushDataEqualById
{
    bool operator()( const QgsSpatialIndexKDBushData &a, const QgsSpatialIndexKDBushData &b ) const { return a.id == b.id; }
};
 
struct KDBushDataHashById
{
    std::size_t operator()( const QgsSpatialIndexKDBushData &a ) const { return std::hash<QgsFeatureId> {}( a.id ); }
};
 
QVariantMap QgsDbscanClusteringAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
  std::unique_ptr<QgsProcessingFeatureSource> source( parameterAsSource( parameters, u"INPUT"_s, context ) );
  if ( !source )
    throw QgsProcessingException( invalidSourceError( parameters, u"INPUT"_s ) );
 
  const std::size_t minSize = static_cast<std::size_t>( parameterAsInt( parameters, u"MIN_SIZE"_s, context ) );
  const double eps1 = parameterAsDouble( parameters, u"EPS"_s, context );
  const double eps2 = parameterAsDouble( parameters, u"EPS2"_s, context );
  const bool borderPointsAreNoise = parameterAsBoolean( parameters, u"DBSCAN*"_s, context );
 
  QgsFields outputFields = source->fields();
  QgsFields newFields;
  const QString clusterFieldName = parameterAsString( parameters, u"FIELD_NAME"_s, context );
  newFields.append( QgsField( clusterFieldName, QMetaType::Type::Int ) );
  const QString clusterSizeFieldName = parameterAsString( parameters, u"SIZE_FIELD_NAME"_s, context );
  newFields.append( QgsField( clusterSizeFieldName, QMetaType::Type::Int ) );
  outputFields = QgsProcessingUtils::combineFields( outputFields, newFields );
 
  QString dest;
  std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, u"OUTPUT"_s, context, dest, outputFields, source->wkbType(), source->sourceCrs() ) );
  if ( !sink )
    throw QgsProcessingException( invalidSinkError( parameters, u"OUTPUT"_s ) );
 
  QgsFeatureRequest indexRequest;
 
  std::unordered_map<QgsFeatureId, QDateTime> idToDateTime;
  const QString dateTimeFieldName = parameterAsString( parameters, u"DATETIME_FIELD"_s, context );
  int dateTimefieldIndex = -1;
  if ( !dateTimeFieldName.isEmpty() )
  {
    dateTimefieldIndex = source->fields().lookupField( dateTimeFieldName );
    if ( dateTimefieldIndex == -1 )
      throw QgsProcessingException( QObject::tr( "Datetime field missing" ) );
 
    indexRequest.setSubsetOfAttributes( QgsAttributeList() << dateTimefieldIndex );
  }
  else
  {
    indexRequest.setNoAttributes();
  }
 
  // build spatial index, also collecting feature datetimes if required
  feedback->pushInfo( QObject::tr( "Building spatial index" ) );
  QgsFeatureIterator indexIterator = source->getFeatures( indexRequest );
  QgsSpatialIndexKDBush index(
    indexIterator,
    [&idToDateTime, dateTimefieldIndex]( const QgsFeature &feature ) -> bool {
      if ( dateTimefieldIndex >= 0 )
        idToDateTime[feature.id()] = feature.attributes().at( dateTimefieldIndex ).toDateTime();
      return true;
    },
    feedback
  );
 
  if ( feedback->isCanceled() )
    return QVariantMap();
 
  // stdbscan!
  feedback->pushInfo( QObject::tr( "Analysing clusters" ) );
  std::unordered_map<QgsFeatureId, int> idToCluster;
  idToCluster.reserve( index.size() );
  const long featureCount = source->featureCount();
  QgsFeatureIterator features = source->getFeatures( QgsFeatureRequest().setNoAttributes() );
  stdbscan( minSize, eps1, eps2, borderPointsAreNoise, featureCount, features, index, idToCluster, idToDateTime, feedback );
 
  // cluster size
  std::unordered_map<int, int> clusterSize;
  std::for_each( idToCluster.begin(), idToCluster.end(), [&clusterSize]( std::pair<QgsFeatureId, int> idCluster ) { clusterSize[idCluster.second]++; } );
 
  // write clusters
  const double writeStep = featureCount > 0 ? 10.0 / featureCount : 1;
  features = source->getFeatures();
  int i = 0;
  QgsFeature feat;
  while ( features.nextFeature( feat ) )
  {
    i++;
    if ( feedback->isCanceled() )
    {
      break;
    }
 
    feedback->setProgress( 90 + i * writeStep );
    QgsAttributes attr = feat.attributes();
    const auto cluster = idToCluster.find( feat.id() );
    if ( cluster != idToCluster.end() )
    {
      attr << cluster->second << clusterSize[cluster->second];
    }
    else
    {
      attr << QVariant() << QVariant();
    }
    feat.setAttributes( attr );
    if ( !sink->addFeature( feat, QgsFeatureSink::FastInsert ) )
      throw QgsProcessingException( writeFeatureError( sink.get(), parameters, u"OUTPUT"_s ) );
  }
 
  sink->finalize();
 
  QVariantMap outputs;
  outputs.insert( u"OUTPUT"_s, dest );
  outputs.insert( u"NUM_CLUSTERS"_s, static_cast<unsigned int>( clusterSize.size() ) );
  return outputs;
}
 
void QgsDbscanClusteringAlgorithm::stdbscan(
  const std::size_t minSize,
  const double eps1,
  const double eps2,
  const bool borderPointsAreNoise,
  const long featureCount,
  QgsFeatureIterator features,
  QgsSpatialIndexKDBush &index,
  std::unordered_map<QgsFeatureId, int> &idToCluster,
  std::unordered_map<QgsFeatureId, QDateTime> &idToDateTime,
  QgsProcessingFeedback *feedback
)
{
  const double step = featureCount > 0 ? 90.0 / featureCount : 1;
 
  std::unordered_set<QgsFeatureId> visited;
  visited.reserve( index.size() );
 
  QgsFeature feat;
  int i = 0;
  int clusterCount = 0;
 
  while ( features.nextFeature( feat ) )
  {
    if ( feedback->isCanceled() )
    {
      break;
    }
 
    if ( !feat.hasGeometry() )
    {
      feedback->setProgress( ++i * step );
      continue;
    }
 
    if ( visited.find( feat.id() ) != visited.end() )
    {
      // already visited!
      continue;
    }
 
    QgsPointXY point;
    if ( QgsWkbTypes::flatType( feat.geometry().wkbType() ) == Qgis::WkbType::Point )
      point = QgsPointXY( *qgsgeometry_cast<const QgsPoint *>( feat.geometry().constGet() ) );
    else
    {
      // not a point geometry
      feedback->reportError( QObject::tr( "Feature %1 is a %2 feature, not a point." ).arg( feat.id() ).arg( QgsWkbTypes::displayString( feat.geometry().wkbType() ) ) );
      feedback->setProgress( ++i * step );
      continue;
    }
 
    if ( !idToDateTime.empty() && !idToDateTime[feat.id()].isValid() )
    {
      // missing datetime value
      feedback->reportError( QObject::tr( "Feature %1 is missing a valid datetime value." ).arg( feat.id() ).arg( QgsWkbTypes::displayString( feat.geometry().wkbType() ) ) );
      feedback->setProgress( ++i * step );
      continue;
    }
 
    std::unordered_set<QgsSpatialIndexKDBushData, KDBushDataHashById, KDBushDataEqualById> within;
 
    if ( minSize > 1 )
    {
      index.within( point, eps1, [&within, pointId = feat.id(), &idToDateTime, &eps2]( const QgsSpatialIndexKDBushData &data ) {
        if ( idToDateTime.empty() || ( idToDateTime[data.id].isValid() && std::abs( idToDateTime[pointId].msecsTo( idToDateTime[data.id] ) ) <= eps2 ) )
          within.insert( data );
      } );
      if ( within.size() < minSize )
        continue;
 
      visited.insert( feat.id() );
    }
    else
    {
      // optimised case for minSize == 1, we can skip the initial check
      within.insert( QgsSpatialIndexKDBushData( feat.id(), point.x(), point.y() ) );
    }
 
    // start new cluster
    clusterCount++;
    idToCluster[feat.id()] = clusterCount;
    feedback->setProgress( ++i * step );
 
    while ( !within.empty() )
    {
      if ( feedback->isCanceled() )
      {
        break;
      }
 
      const QgsSpatialIndexKDBushData j = *within.begin();
      within.erase( within.begin() );
 
      if ( visited.find( j.id ) != visited.end() )
      {
        // already visited!
        continue;
      }
 
      visited.insert( j.id );
      feedback->setProgress( ++i * step );
 
      // check from this point
      const QgsPointXY point2 = j.point();
 
      std::unordered_set<QgsSpatialIndexKDBushData, KDBushDataHashById, KDBushDataEqualById> within2;
      index.within( point2, eps1, [&within2, point2Id = j.id, &idToDateTime, &eps2]( const QgsSpatialIndexKDBushData &data ) {
        if ( idToDateTime.empty() || ( idToDateTime[data.id].isValid() && std::abs( idToDateTime[point2Id].msecsTo( idToDateTime[data.id] ) ) <= eps2 ) )
          within2.insert( data );
      } );
 
      if ( within2.size() >= minSize )
      {
        // expand neighbourhood
        std::copy_if( within2.begin(), within2.end(), std::inserter( within, within.end() ), [&visited]( const QgsSpatialIndexKDBushData &needle ) {
          return visited.find( needle.id ) == visited.end();
        } );
      }
      if ( !borderPointsAreNoise || within2.size() >= minSize )
      {
        idToCluster[j.id] = clusterCount;
      }
    }
  }
}