qgsvectorlayerdirector.cpp

Go to the documentation of this file.

/***************************************************************************
  qgslinevectorlayerdirector.cpp
  --------------------------------------
  Date                 : 2010-10-20
  Copyright            : (C) 2010 by Yakushev Sergey
  Email                : [email protected]
****************************************************************************
*                                                                          *
*   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 "qgsvectorlayerdirector.h"
 
#include <spatialindex/SpatialIndex.h>
 
#include "qgsdistancearea.h"
#include "qgsfeatureiterator.h"
#include "qgsfeaturesource.h"
#include "qgsgeometry.h"
#include "qgsgraphbuilderinterface.h"
#include "qgslogger.h"
#include "qgsvectordataprovider.h"
#include "qgswkbtypes.h"
 
#include <QString>
#include <QtAlgorithms>
 
#include "moc_qgsvectorlayerdirector.cpp"
 
using namespace Qt::StringLiterals;
 
using namespace SpatialIndex;
 
struct TiePointInfo
{
    TiePointInfo() = default;
    TiePointInfo( int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end )
      : additionalPointId( additionalPointId )
      , mNetworkFeatureId( featureId )
      , mFirstPoint( start )
      , mLastPoint( end )
    {}
 
    int additionalPointId = -1;
    QgsPointXY mTiedPoint;
    double mLength = std::numeric_limits<double>::max();
    QgsFeatureId mNetworkFeatureId = -1;
    QgsPointXY mFirstPoint;
    QgsPointXY mLastPoint;
};
 
QgsVectorLayerDirector::QgsVectorLayerDirector( QgsFeatureSource *source, int directionFieldId, const QString &directDirectionValue, const QString &reverseDirectionValue, const QString &bothDirectionValue, const Direction defaultDirection )
  : mSource( source )
  , mDirectionFieldId( directionFieldId )
  , mDirectDirectionValue( directDirectionValue )
  , mReverseDirectionValue( reverseDirectionValue )
  , mBothDirectionValue( bothDirectionValue )
  , mDefaultDirection( defaultDirection )
{
}
 
QString QgsVectorLayerDirector::name() const
{
  return u"Vector line"_s;
}
 
QgsAttributeList QgsVectorLayerDirector::requiredAttributes() const
{
  QSet<int> attrs;
 
  if ( mDirectionFieldId != -1 )
    attrs.insert( mDirectionFieldId );
 
  for ( const QgsNetworkStrategy *strategy : mStrategies )
  {
    attrs.unite( strategy->requiredAttributes() );
  }
  return qgis::setToList( attrs );
}
 
QgsVectorLayerDirector::Direction QgsVectorLayerDirector::directionForFeature( const QgsFeature &feature ) const
{
  if ( mDirectionFieldId < 0 )
    return mDefaultDirection;
 
  QString str = feature.attribute( mDirectionFieldId ).toString();
  if ( str == mBothDirectionValue )
  {
    return Direction::DirectionBoth;
  }
  else if ( str == mDirectDirectionValue )
  {
    return Direction::DirectionForward;
  }
  else if ( str == mReverseDirectionValue )
  {
    return Direction::DirectionBackward;
  }
  else
  {
    return mDefaultDirection;
  }
}

class QgsNetworkVisitor : public SpatialIndex::IVisitor
{
  public:
    explicit QgsNetworkVisitor( QVector<int> &pointIndexes )
      : mPoints( pointIndexes ) {}
 
    void visitNode( const INode &n ) override
    {
      Q_UNUSED( n )
    }
 
    void visitData( const IData &d ) override
    {
      mPoints.append( static_cast<int>( d.getIdentifier() ) );
    }
 
    void visitData( std::vector<const IData *> &v ) override
    {
      Q_UNUSED( v )
    }
 
  private:
    QVector<int> &mPoints;
};

 
std::unique_ptr<SpatialIndex::ISpatialIndex> createVertexSpatialIndex( SpatialIndex::IStorageManager &storageManager )
{
  // R-Tree parameters
  double fillFactor = 0.7;
  unsigned long indexCapacity = 10;
  unsigned long leafCapacity = 10;
  unsigned long dimension = 2;
  RTree::RTreeVariant variant = RTree::RV_RSTAR;
 
  // create R-tree
  SpatialIndex::id_type indexId;
  std::unique_ptr<SpatialIndex::ISpatialIndex> iRTree( RTree::createNewRTree( storageManager, fillFactor, indexCapacity, leafCapacity, dimension, variant, indexId ) );
  return iRTree;
}
 
int findClosestVertex( const QgsPointXY &point, SpatialIndex::ISpatialIndex *index, double tolerance )
{
  QVector<int> matching;
  QgsNetworkVisitor visitor( matching );
 
  double pt1[2] = { point.x() - tolerance, point.y() - tolerance },
         pt2[2] = { point.x() + tolerance, point.y() + tolerance };
  SpatialIndex::Region searchRegion( pt1, pt2, 2 );
 
  index->intersectsWithQuery( searchRegion, visitor );
 
  return matching.empty() ? -1 : matching.at( 0 );
}
 
void QgsVectorLayerDirector::makeGraph( QgsGraphBuilderInterface *builder, const QVector<QgsPointXY> &additionalPoints, QVector<QgsPointXY> &snappedPoints, QgsFeedback *feedback ) const
{
  long featureCount = mSource->featureCount() * 2;
  int step = 0;
 
  QgsCoordinateTransform ct;
  ct.setSourceCrs( mSource->sourceCrs() );
  if ( builder->coordinateTransformationEnabled() )
  {
    ct.setDestinationCrs( builder->destinationCrs() );
  }
 
  // clear existing snapped points list, and resize to length of provided additional points
  snappedPoints = QVector<QgsPointXY>( additionalPoints.size(), QgsPointXY( 0.0, 0.0 ) );
  // tie points = snapped location of specified additional points to network lines
  QVector<TiePointInfo> additionalTiePoints( additionalPoints.size() );
 
  // graph's vertices = all vertices in graph, with vertices within builder's tolerance collapsed together
  QVector<QgsPointXY> graphVertices;
 
  // spatial index for graph vertices
  std::unique_ptr<SpatialIndex::IStorageManager> iStorage( StorageManager::createNewMemoryStorageManager() );
  std::unique_ptr<SpatialIndex::ISpatialIndex> iRTree = createVertexSpatialIndex( *iStorage );
 
  double tolerance = std::max( builder->topologyTolerance(), 1e-10 );
  auto findPointWithinTolerance = [&iRTree, tolerance]( const QgsPointXY &point ) -> int {
    return findClosestVertex( point, iRTree.get(), tolerance );
  };
  auto addPointToIndex = [&iRTree]( const QgsPointXY &point, int index ) {
    double coords[] = { point.x(), point.y() };
    iRTree->insertData( 0, nullptr, SpatialIndex::Point( coords, 2 ), index );
  };
 
  // first iteration - get all nodes from network, and snap additional points to network
  QgsFeatureIterator fit = mSource->getFeatures( QgsFeatureRequest().setNoAttributes() );
  QgsFeature feature;
  while ( fit.nextFeature( feature ) )
  {
    if ( feedback && feedback->isCanceled() )
      return;
 
    QgsMultiPolylineXY mpl;
    const Qgis::WkbType wkbType = QgsWkbTypes::flatType( feature.geometry().wkbType() );
    if ( QgsWkbTypes::geometryType( wkbType ) == Qgis::GeometryType::Line && QgsWkbTypes::isMultiType( wkbType ) )
    {
      mpl = feature.geometry().asMultiPolyline();
    }
    else if ( QgsWkbTypes::geometryType( wkbType ) == Qgis::GeometryType::Line && !QgsWkbTypes::isMultiType( wkbType ) )
    {
      mpl.push_back( feature.geometry().asPolyline() );
    }
 
    for ( const QgsPolylineXY &line : std::as_const( mpl ) )
    {
      QgsPointXY pt1, pt2;
      bool isFirstPoint = true;
      for ( const QgsPointXY &point : line )
      {
        pt2 = ct.transform( point );
 
        int pt2Idx = findPointWithinTolerance( pt2 );
        if ( pt2Idx == -1 )
        {
          // no vertex already exists within tolerance - add to points, and index
          addPointToIndex( pt2, static_cast<int>( graphVertices.count() ) );
          graphVertices.push_back( pt2 );
        }
        else
        {
          // vertex already exists within tolerance - use that
          pt2 = graphVertices.at( pt2Idx );
        }
 
        if ( !isFirstPoint )
        {
          // check if this line segment is a candidate for being closest to each additional point
          int i = 0;
          for ( const QgsPointXY &additionalPoint : additionalPoints )
          {
            QgsPointXY snappedPoint;
            double thisSegmentClosestDist = std::numeric_limits<double>::max();
            if ( pt1 == pt2 )
            {
              thisSegmentClosestDist = additionalPoint.sqrDist( pt1 );
              snappedPoint = pt1;
            }
            else
            {
              thisSegmentClosestDist = additionalPoint.sqrDistToSegment( pt1.x(), pt1.y(), pt2.x(), pt2.y(), snappedPoint, 0 );
            }
 
            if ( thisSegmentClosestDist < additionalTiePoints[i].mLength )
            {
              // found a closer segment for this additional point
              TiePointInfo info( i, feature.id(), pt1, pt2 );
              info.mLength = thisSegmentClosestDist;
              info.mTiedPoint = snappedPoint;
 
              additionalTiePoints[i] = info;
              snappedPoints[i] = info.mTiedPoint;
            }
            i++;
          }
        }
        pt1 = pt2;
        isFirstPoint = false;
      }
    }
    if ( feedback )
      feedback->setProgress( 100.0 * static_cast<double>( ++step ) / static_cast<double>( featureCount ) );
  }
 
  // build a hash of feature ids to tie points which depend on this feature
  QHash<QgsFeatureId, QList<int>> tiePointNetworkFeatures;
  int i = 0;
  for ( TiePointInfo &info : additionalTiePoints )
  {
    tiePointNetworkFeatures[info.mNetworkFeatureId] << i;
    i++;
  }
 
  // add tied point to graph
  for ( int i = 0; i < snappedPoints.size(); ++i )
  {
    // check index to see if vertex exists within tolerance of tie point
    const QgsPointXY point = snappedPoints.at( i );
    int ptIdx = findPointWithinTolerance( point );
    if ( ptIdx == -1 )
    {
      // no vertex already within tolerance, add to index and network vertices
      addPointToIndex( point, static_cast<int>( graphVertices.count() ) );
      graphVertices.push_back( point );
    }
    else
    {
      // otherwise snap tie point to vertex
      snappedPoints[i] = graphVertices.at( ptIdx );
    }
  }
  // also need to update tie points - they need to be matched for snapped points
  for ( int i = 0; i < additionalTiePoints.count(); ++i )
  {
    additionalTiePoints[i].mTiedPoint = snappedPoints.at( additionalTiePoints.at( i ).additionalPointId );
  }
 
  // begin graph construction
 
  // add vertices to graph
  {
    int i = 0;
    for ( const QgsPointXY &point : graphVertices )
    {
      builder->addVertex( i, point );
      i++;
    }
  }
 
  fit = mSource->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes() ) );
  while ( fit.nextFeature( feature ) )
  {
    if ( feedback && feedback->isCanceled() )
      return;
 
    Direction direction = directionForFeature( feature );
 
    // begin features segments and add arc to the Graph;
    QgsMultiPolylineXY mpl;
    const Qgis::WkbType wkbType = QgsWkbTypes::flatType( feature.geometry().wkbType() );
    if ( QgsWkbTypes::geometryType( wkbType ) == Qgis::GeometryType::Line && QgsWkbTypes::isMultiType( wkbType ) )
    {
      mpl = feature.geometry().asMultiPolyline();
    }
    else if ( QgsWkbTypes::geometryType( wkbType ) == Qgis::GeometryType::Line && !QgsWkbTypes::isMultiType( wkbType ) )
    {
      mpl.push_back( feature.geometry().asPolyline() );
    }
 
    for ( const QgsPolylineXY &line : std::as_const( mpl ) )
    {
      QgsPointXY pt1, pt2;
 
      bool isFirstPoint = true;
      for ( const QgsPointXY &point : line )
      {
        pt2 = ct.transform( point );
        int pPt2idx = findPointWithinTolerance( pt2 );
        Q_ASSERT_X( pPt2idx >= 0, "QgsVectorLayerDirector::makeGraph", "encountered a vertex which was not present in graph" );
        pt2 = graphVertices.at( pPt2idx );
 
        if ( !isFirstPoint )
        {
          QMap<double, QgsPointXY> pointsOnArc;
          pointsOnArc[0.0] = pt1;
          pointsOnArc[pt1.sqrDist( pt2 )] = pt2;
 
          const QList<int> tiePointsForCurrentFeature = tiePointNetworkFeatures.value( feature.id() );
          for ( int tiePointIdx : tiePointsForCurrentFeature )
          {
            const TiePointInfo &t = additionalTiePoints.at( tiePointIdx );
            if ( t.mFirstPoint == pt1 && t.mLastPoint == pt2 )
            {
              pointsOnArc[pt1.sqrDist( t.mTiedPoint )] = t.mTiedPoint;
            }
          }
 
          QgsPointXY arcPt1;
          QgsPointXY arcPt2;
          int pt1idx = -1;
          int pt2idx = -1;
          bool isFirstPoint = true;
          for ( auto arcPointIt = pointsOnArc.constBegin(); arcPointIt != pointsOnArc.constEnd(); ++arcPointIt )
          {
            arcPt2 = arcPointIt.value();
 
            pt2idx = findPointWithinTolerance( arcPt2 );
            Q_ASSERT_X( pt2idx >= 0, "QgsVectorLayerDirector::makeGraph", "encountered a vertex which was not present in graph" );
            arcPt2 = graphVertices.at( pt2idx );
 
            if ( !isFirstPoint && arcPt1 != arcPt2 )
            {
              double distance = 0;
              try
              {
                distance = builder->distanceArea()->measureLine( arcPt1, arcPt2 );
              }
              catch ( QgsCsException & )
              {
                // TODO report errors to user
                QgsDebugError( u"An error occurred while calculating length"_s );
              }
 
              QVector<QVariant> prop;
              prop.reserve( mStrategies.size() );
              for ( QgsNetworkStrategy *strategy : mStrategies )
              {
                prop.push_back( strategy->cost( distance, feature ) );
              }
 
              if ( direction == Direction::DirectionForward || direction == Direction::DirectionBoth )
              {
                builder->addEdge( pt1idx, arcPt1, pt2idx, arcPt2, prop );
              }
              if ( direction == Direction::DirectionBackward || direction == Direction::DirectionBoth )
              {
                builder->addEdge( pt2idx, arcPt2, pt1idx, arcPt1, prop );
              }
            }
            pt1idx = pt2idx;
            arcPt1 = arcPt2;
            isFirstPoint = false;
          }
        }
        pt1 = pt2;
        isFirstPoint = false;
      }
    }
    if ( feedback )
    {
      feedback->setProgress( 100.0 * static_cast<double>( ++step ) / static_cast<double>( featureCount ) );
    }
  }
}