Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgslinevectorlayerdirector.cpp
3 : : --------------------------------------
4 : : Date : 2010-10-20
5 : : Copyright : (C) 2010 by Yakushev Sergey
6 : : Email : YakushevS@list.ru
7 : : ****************************************************************************
8 : : * *
9 : : * This program is free software; you can redistribute it and/or modify *
10 : : * it under the terms of the GNU General Public License as published by *
11 : : * the Free Software Foundation; either version 2 of the License, or *
12 : : * (at your option) any later version. *
13 : : * *
14 : : ***************************************************************************/
15 : :
16 : : /**
17 : : * \file qgsvectorlayerdirector.cpp
18 : : * \brief implementation of QgsVectorLayerDirector
19 : : */
20 : :
21 : : #include "qgsvectorlayerdirector.h"
22 : : #include "qgsgraphbuilderinterface.h"
23 : :
24 : : #include "qgsfeatureiterator.h"
25 : : #include "qgsfeaturesource.h"
26 : : #include "qgsvectordataprovider.h"
27 : : #include "qgspoint.h"
28 : : #include "qgsgeometry.h"
29 : : #include "qgsdistancearea.h"
30 : : #include "qgswkbtypes.h"
31 : :
32 : : #include <QString>
33 : : #include <QtAlgorithms>
34 : :
35 : : #include <spatialindex/SpatialIndex.h>
36 : :
37 : : using namespace SpatialIndex;
38 : :
39 : 0 : struct TiePointInfo
40 : : {
41 : 0 : TiePointInfo() = default;
42 : 0 : TiePointInfo( int additionalPointId, QgsFeatureId featureId, const QgsPointXY &start, const QgsPointXY &end )
43 : 0 : : additionalPointId( additionalPointId )
44 : 0 : , mNetworkFeatureId( featureId )
45 : 0 : , mFirstPoint( start )
46 : 0 : , mLastPoint( end )
47 : 0 : {}
48 : :
49 : 0 : int additionalPointId = -1;
50 : : QgsPointXY mTiedPoint;
51 : 0 : double mLength = std::numeric_limits<double>::max();
52 : 0 : QgsFeatureId mNetworkFeatureId = -1;
53 : : QgsPointXY mFirstPoint;
54 : : QgsPointXY mLastPoint;
55 : : };
56 : :
57 : 0 : QgsVectorLayerDirector::QgsVectorLayerDirector( QgsFeatureSource *source,
58 : : int directionFieldId,
59 : : const QString &directDirectionValue,
60 : : const QString &reverseDirectionValue,
61 : : const QString &bothDirectionValue,
62 : : const Direction defaultDirection )
63 : 0 : : mSource( source )
64 : 0 : , mDirectionFieldId( directionFieldId )
65 : 0 : , mDirectDirectionValue( directDirectionValue )
66 : 0 : , mReverseDirectionValue( reverseDirectionValue )
67 : 0 : , mBothDirectionValue( bothDirectionValue )
68 : 0 : , mDefaultDirection( defaultDirection )
69 : 0 : {
70 : 0 : }
71 : :
72 : 0 : QString QgsVectorLayerDirector::name() const
73 : : {
74 : 0 : return QStringLiteral( "Vector line" );
75 : : }
76 : :
77 : 0 : QgsAttributeList QgsVectorLayerDirector::requiredAttributes() const
78 : : {
79 : 0 : QSet< int > attrs;
80 : :
81 : 0 : if ( mDirectionFieldId != -1 )
82 : 0 : attrs.insert( mDirectionFieldId );
83 : :
84 : 0 : for ( const QgsNetworkStrategy *strategy : mStrategies )
85 : : {
86 : 0 : attrs.unite( strategy->requiredAttributes() );
87 : : }
88 : 0 : return qgis::setToList( attrs );
89 : 0 : }
90 : :
91 : 0 : QgsVectorLayerDirector::Direction QgsVectorLayerDirector::directionForFeature( const QgsFeature &feature ) const
92 : : {
93 : 0 : if ( mDirectionFieldId < 0 )
94 : 0 : return mDefaultDirection;
95 : :
96 : 0 : QString str = feature.attribute( mDirectionFieldId ).toString();
97 : 0 : if ( str == mBothDirectionValue )
98 : : {
99 : 0 : return Direction::DirectionBoth;
100 : : }
101 : 0 : else if ( str == mDirectDirectionValue )
102 : : {
103 : 0 : return Direction::DirectionForward;
104 : : }
105 : 0 : else if ( str == mReverseDirectionValue )
106 : : {
107 : 0 : return Direction::DirectionBackward;
108 : : }
109 : : else
110 : : {
111 : 0 : return mDefaultDirection;
112 : : }
113 : 0 : }
114 : :
115 : : ///@cond PRIVATE
116 : 0 : class QgsNetworkVisitor : public SpatialIndex::IVisitor
117 : : {
118 : : public:
119 : 0 : explicit QgsNetworkVisitor( QVector< int > &pointIndexes )
120 : 0 : : mPoints( pointIndexes ) {}
121 : :
122 : 0 : void visitNode( const INode &n ) override
123 : 0 : { Q_UNUSED( n ) }
124 : :
125 : 0 : void visitData( const IData &d ) override
126 : : {
127 : 0 : mPoints.append( d.getIdentifier() );
128 : 0 : }
129 : :
130 : 0 : void visitData( std::vector<const IData *> &v ) override
131 : 0 : { Q_UNUSED( v ) }
132 : :
133 : : private:
134 : : QVector< int > &mPoints;
135 : : };
136 : :
137 : : ///@endcond
138 : :
139 : 0 : std::unique_ptr< SpatialIndex::ISpatialIndex > createVertexSpatialIndex( SpatialIndex::IStorageManager &storageManager )
140 : : {
141 : : // R-Tree parameters
142 : 0 : double fillFactor = 0.7;
143 : 0 : unsigned long indexCapacity = 10;
144 : 0 : unsigned long leafCapacity = 10;
145 : 0 : unsigned long dimension = 2;
146 : 0 : RTree::RTreeVariant variant = RTree::RV_RSTAR;
147 : :
148 : : // create R-tree
149 : : SpatialIndex::id_type indexId;
150 : 0 : std::unique_ptr< SpatialIndex::ISpatialIndex > iRTree( RTree::createNewRTree( storageManager, fillFactor, indexCapacity,
151 : 0 : leafCapacity, dimension, variant, indexId ) );
152 : 0 : return iRTree;
153 : 0 : }
154 : :
155 : 0 : int findClosestVertex( const QgsPointXY &point, SpatialIndex::ISpatialIndex *index, double tolerance )
156 : : {
157 : 0 : QVector< int > matching;
158 : 0 : QgsNetworkVisitor visitor( matching );
159 : :
160 : 0 : double pt1[2] = { point.x() - tolerance, point.y() - tolerance },
161 : 0 : pt2[2] = { point.x() + tolerance, point.y() + tolerance };
162 : 0 : SpatialIndex::Region searchRegion( pt1, pt2, 2 );
163 : :
164 : 0 : index->intersectsWithQuery( searchRegion, visitor );
165 : :
166 : 0 : return matching.empty() ? -1 : matching.at( 0 );
167 : 0 : }
168 : :
169 : 0 : void QgsVectorLayerDirector::makeGraph( QgsGraphBuilderInterface *builder, const QVector< QgsPointXY > &additionalPoints,
170 : : QVector< QgsPointXY > &snappedPoints, QgsFeedback *feedback ) const
171 : : {
172 : 0 : long featureCount = mSource->featureCount() * 2;
173 : 0 : int step = 0;
174 : :
175 : 0 : QgsCoordinateTransform ct;
176 : 0 : ct.setSourceCrs( mSource->sourceCrs() );
177 : 0 : if ( builder->coordinateTransformationEnabled() )
178 : : {
179 : 0 : ct.setDestinationCrs( builder->destinationCrs() );
180 : 0 : }
181 : :
182 : : // clear existing snapped points list, and resize to length of provided additional points
183 : 0 : snappedPoints = QVector< QgsPointXY >( additionalPoints.size(), QgsPointXY( 0.0, 0.0 ) );
184 : : // tie points = snapped location of specified additional points to network lines
185 : 0 : QVector< TiePointInfo > additionalTiePoints( additionalPoints.size() );
186 : :
187 : : // graph's vertices = all vertices in graph, with vertices within builder's tolerance collapsed together
188 : 0 : QVector< QgsPointXY > graphVertices;
189 : :
190 : : // spatial index for graph vertices
191 : 0 : std::unique_ptr< SpatialIndex::IStorageManager > iStorage( StorageManager::createNewMemoryStorageManager() );
192 : 0 : std::unique_ptr< SpatialIndex::ISpatialIndex > iRTree = createVertexSpatialIndex( *iStorage );
193 : :
194 : 0 : double tolerance = std::max( builder->topologyTolerance(), 1e-10 );
195 : 0 : auto findPointWithinTolerance = [&iRTree, tolerance]( const QgsPointXY & point )->int
196 : : {
197 : 0 : return findClosestVertex( point, iRTree.get(), tolerance );
198 : : };
199 : 0 : auto addPointToIndex = [&iRTree]( const QgsPointXY & point, int index )
200 : : {
201 : 0 : double coords[] = {point.x(), point.y()};
202 : 0 : iRTree->insertData( 0, nullptr, SpatialIndex::Point( coords, 2 ), index );
203 : 0 : };
204 : :
205 : : // first iteration - get all nodes from network, and snap additional points to network
206 : 0 : QgsFeatureIterator fit = mSource->getFeatures( QgsFeatureRequest().setNoAttributes() );
207 : 0 : QgsFeature feature;
208 : 0 : while ( fit.nextFeature( feature ) )
209 : : {
210 : 0 : if ( feedback && feedback->isCanceled() )
211 : 0 : return;
212 : :
213 : 0 : QgsMultiPolylineXY mpl;
214 : 0 : if ( QgsWkbTypes::flatType( feature.geometry().wkbType() ) == QgsWkbTypes::MultiLineString )
215 : 0 : mpl = feature.geometry().asMultiPolyline();
216 : 0 : else if ( QgsWkbTypes::flatType( feature.geometry().wkbType() ) == QgsWkbTypes::LineString )
217 : 0 : mpl.push_back( feature.geometry().asPolyline() );
218 : :
219 : 0 : for ( const QgsPolylineXY &line : std::as_const( mpl ) )
220 : : {
221 : 0 : QgsPointXY pt1, pt2;
222 : 0 : bool isFirstPoint = true;
223 : 0 : for ( const QgsPointXY &point : line )
224 : : {
225 : 0 : pt2 = ct.transform( point );
226 : :
227 : 0 : int pt2Idx = findPointWithinTolerance( pt2 ) ;
228 : 0 : if ( pt2Idx == -1 )
229 : : {
230 : : // no vertex already exists within tolerance - add to points, and index
231 : 0 : addPointToIndex( pt2, graphVertices.count() );
232 : 0 : graphVertices.push_back( pt2 );
233 : 0 : }
234 : : else
235 : : {
236 : : // vertex already exists within tolerance - use that
237 : 0 : pt2 = graphVertices.at( pt2Idx );
238 : : }
239 : :
240 : 0 : if ( !isFirstPoint )
241 : : {
242 : : // check if this line segment is a candidate for being closest to each additional point
243 : 0 : int i = 0;
244 : 0 : for ( const QgsPointXY &additionalPoint : additionalPoints )
245 : : {
246 : :
247 : 0 : QgsPointXY snappedPoint;
248 : 0 : double thisSegmentClosestDist = std::numeric_limits<double>::max();
249 : 0 : if ( pt1 == pt2 )
250 : : {
251 : 0 : thisSegmentClosestDist = additionalPoint.sqrDist( pt1 );
252 : 0 : snappedPoint = pt1;
253 : 0 : }
254 : : else
255 : : {
256 : 0 : thisSegmentClosestDist = additionalPoint.sqrDistToSegment( pt1.x(), pt1.y(),
257 : 0 : pt2.x(), pt2.y(), snappedPoint, 0 );
258 : : }
259 : :
260 : 0 : if ( thisSegmentClosestDist < additionalTiePoints[ i ].mLength )
261 : : {
262 : : // found a closer segment for this additional point
263 : 0 : TiePointInfo info( i, feature.id(), pt1, pt2 );
264 : 0 : info.mLength = thisSegmentClosestDist;
265 : 0 : info.mTiedPoint = snappedPoint;
266 : :
267 : 0 : additionalTiePoints[ i ] = info;
268 : 0 : snappedPoints[ i ] = info.mTiedPoint;
269 : 0 : }
270 : 0 : i++;
271 : : }
272 : 0 : }
273 : 0 : pt1 = pt2;
274 : 0 : isFirstPoint = false;
275 : : }
276 : : }
277 : 0 : if ( feedback )
278 : 0 : feedback->setProgress( 100.0 * static_cast< double >( ++step ) / featureCount );
279 : 0 : }
280 : :
281 : : // build a hash of feature ids to tie points which depend on this feature
282 : 0 : QHash< QgsFeatureId, QList< int > > tiePointNetworkFeatures;
283 : 0 : int i = 0;
284 : 0 : for ( TiePointInfo &info : additionalTiePoints )
285 : : {
286 : 0 : tiePointNetworkFeatures[ info.mNetworkFeatureId ] << i;
287 : 0 : i++;
288 : : }
289 : :
290 : : // add tied point to graph
291 : 0 : for ( int i = 0; i < snappedPoints.size(); ++i )
292 : : {
293 : : // check index to see if vertex exists within tolerance of tie point
294 : 0 : const QgsPointXY point = snappedPoints.at( i );
295 : 0 : int ptIdx = findPointWithinTolerance( point );
296 : 0 : if ( ptIdx == -1 )
297 : : {
298 : : // no vertex already within tolerance, add to index and network vertices
299 : 0 : addPointToIndex( point, graphVertices.count() );
300 : 0 : graphVertices.push_back( point );
301 : 0 : }
302 : : else
303 : : {
304 : : // otherwise snap tie point to vertex
305 : 0 : snappedPoints[ i ] = graphVertices.at( ptIdx );
306 : : }
307 : 0 : }
308 : : // also need to update tie points - they need to be matched for snapped points
309 : 0 : for ( int i = 0; i < additionalTiePoints.count(); ++i )
310 : : {
311 : 0 : additionalTiePoints[ i ].mTiedPoint = snappedPoints.at( additionalTiePoints.at( i ).additionalPointId );
312 : 0 : }
313 : :
314 : :
315 : : // begin graph construction
316 : :
317 : : // add vertices to graph
318 : : {
319 : 0 : int i = 0;
320 : 0 : for ( const QgsPointXY &point : graphVertices )
321 : : {
322 : 0 : builder->addVertex( i, point );
323 : 0 : i++;
324 : : }
325 : : }
326 : :
327 : 0 : fit = mSource->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes() ) );
328 : 0 : while ( fit.nextFeature( feature ) )
329 : : {
330 : 0 : if ( feedback && feedback->isCanceled() )
331 : 0 : return;
332 : :
333 : 0 : Direction direction = directionForFeature( feature );
334 : :
335 : : // begin features segments and add arc to the Graph;
336 : 0 : QgsMultiPolylineXY mpl;
337 : 0 : if ( QgsWkbTypes::flatType( feature.geometry().wkbType() ) == QgsWkbTypes::MultiLineString )
338 : 0 : mpl = feature.geometry().asMultiPolyline();
339 : 0 : else if ( QgsWkbTypes::flatType( feature.geometry().wkbType() ) == QgsWkbTypes::LineString )
340 : 0 : mpl.push_back( feature.geometry().asPolyline() );
341 : :
342 : 0 : for ( const QgsPolylineXY &line : std::as_const( mpl ) )
343 : : {
344 : 0 : QgsPointXY pt1, pt2;
345 : :
346 : 0 : bool isFirstPoint = true;
347 : 0 : for ( const QgsPointXY &point : line )
348 : : {
349 : 0 : pt2 = ct.transform( point );
350 : 0 : int pPt2idx = findPointWithinTolerance( pt2 );
351 : : Q_ASSERT_X( pPt2idx >= 0, "QgsVectorLayerDirectory::makeGraph", "encountered a vertex which was not present in graph" );
352 : 0 : pt2 = graphVertices.at( pPt2idx );
353 : :
354 : 0 : if ( !isFirstPoint )
355 : : {
356 : 0 : QMap< double, QgsPointXY > pointsOnArc;
357 : 0 : pointsOnArc[ 0.0 ] = pt1;
358 : 0 : pointsOnArc[ pt1.sqrDist( pt2 )] = pt2;
359 : :
360 : 0 : const QList< int > tiePointsForCurrentFeature = tiePointNetworkFeatures.value( feature.id() );
361 : 0 : for ( int tiePointIdx : tiePointsForCurrentFeature )
362 : : {
363 : 0 : const TiePointInfo &t = additionalTiePoints.at( tiePointIdx );
364 : 0 : if ( t.mFirstPoint == pt1 && t.mLastPoint == pt2 )
365 : : {
366 : 0 : pointsOnArc[ pt1.sqrDist( t.mTiedPoint )] = t.mTiedPoint;
367 : 0 : }
368 : : }
369 : :
370 : 0 : QgsPointXY arcPt1;
371 : 0 : QgsPointXY arcPt2;
372 : 0 : int pt1idx = -1;
373 : 0 : int pt2idx = -1;
374 : 0 : bool isFirstPoint = true;
375 : 0 : for ( auto arcPointIt = pointsOnArc.constBegin(); arcPointIt != pointsOnArc.constEnd(); ++arcPointIt )
376 : : {
377 : 0 : arcPt2 = arcPointIt.value();
378 : :
379 : 0 : pt2idx = findPointWithinTolerance( arcPt2 );
380 : : Q_ASSERT_X( pt2idx >= 0, "QgsVectorLayerDirectory::makeGraph", "encountered a vertex which was not present in graph" );
381 : 0 : arcPt2 = graphVertices.at( pt2idx );
382 : :
383 : 0 : if ( !isFirstPoint && arcPt1 != arcPt2 )
384 : : {
385 : 0 : double distance = builder->distanceArea()->measureLine( arcPt1, arcPt2 );
386 : 0 : QVector< QVariant > prop;
387 : 0 : prop.reserve( mStrategies.size() );
388 : 0 : for ( QgsNetworkStrategy *strategy : mStrategies )
389 : : {
390 : 0 : prop.push_back( strategy->cost( distance, feature ) );
391 : : }
392 : :
393 : 0 : if ( direction == Direction::DirectionForward ||
394 : 0 : direction == Direction::DirectionBoth )
395 : : {
396 : 0 : builder->addEdge( pt1idx, arcPt1, pt2idx, arcPt2, prop );
397 : 0 : }
398 : 0 : if ( direction == Direction::DirectionBackward ||
399 : 0 : direction == Direction::DirectionBoth )
400 : : {
401 : 0 : builder->addEdge( pt2idx, arcPt2, pt1idx, arcPt1, prop );
402 : 0 : }
403 : 0 : }
404 : 0 : pt1idx = pt2idx;
405 : 0 : arcPt1 = arcPt2;
406 : 0 : isFirstPoint = false;
407 : 0 : }
408 : 0 : }
409 : 0 : pt1 = pt2;
410 : 0 : isFirstPoint = false;
411 : : }
412 : : }
413 : 0 : if ( feedback )
414 : : {
415 : 0 : feedback->setProgress( 100.0 * static_cast< double >( ++step ) / featureCount );
416 : 0 : }
417 : :
418 : 0 : }
419 : 0 : }
|
