|
| 1 | +/*************************************************************************** |
| 2 | + qgsgeometrysnappersinglesource.cpp |
| 3 | + --------------------- |
| 4 | + Date : May 2018 |
| 5 | + Copyright : (C) 2018 by Martin Dobias |
| 6 | + Email : wonder dot sk at gmail dot com |
| 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 | +#include "qgsgeometrysnappersinglesource.h" |
| 17 | + |
| 18 | +#include "qgsfeatureiterator.h" |
| 19 | +#include "qgsfeaturesink.h" |
| 20 | +#include "qgsfeaturesource.h" |
| 21 | +#include "qgsfeedback.h" |
| 22 | +#include "qgsgeometrycollection.h" |
| 23 | +#include "qgsgeometryutils.h" |
| 24 | +#include "qgslinestring.h" |
| 25 | +#include "qgspolygon.h" |
| 26 | +#include "qgsspatialindex.h" |
| 27 | + |
| 28 | +//! record about vertex coordinates and index of anchor to which it is snapped |
| 29 | +struct AnchorPoint |
| 30 | +{ |
| 31 | + //! coordinates of the point |
| 32 | + double x, y; |
| 33 | + |
| 34 | + /** |
| 35 | + * Anchor information: |
| 36 | + * 0+ - index of anchor to which this point should be snapped |
| 37 | + * -1 - initial value (undefined) |
| 38 | + * -2 - this point is an anchor, i.e. do not snap this point (snap others to this point) |
| 39 | + */ |
| 40 | + int anchor; |
| 41 | +}; |
| 42 | + |
| 43 | + |
| 44 | +//! record about anchor being along a segment |
| 45 | +struct AnchorAlongSegment |
| 46 | +{ |
| 47 | + int anchor; //!< Index of the anchor point |
| 48 | + double along; //!< Distance of the anchor point along the segment |
| 49 | +}; |
| 50 | + |
| 51 | + |
| 52 | +static void buildSnapIndex( QgsFeatureIterator &fi, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, QgsFeedback *feedback, int &count, int totalCount ) |
| 53 | +{ |
| 54 | + QgsFeature f; |
| 55 | + int pntId = 0; |
| 56 | + |
| 57 | + while ( fi.nextFeature( f ) ) |
| 58 | + { |
| 59 | + if ( feedback->isCanceled() ) |
| 60 | + break; |
| 61 | + |
| 62 | + QgsGeometry g = f.geometry(); |
| 63 | + |
| 64 | + for ( auto it = g.vertices_begin(); it != g.vertices_end(); ++it ) |
| 65 | + { |
| 66 | + QgsPoint pt = *it; |
| 67 | + QgsRectangle rect( pt.x(), pt.y(), pt.x(), pt.y() ); |
| 68 | + |
| 69 | + QList<QgsFeatureId> ids = index.intersects( rect ); |
| 70 | + if ( ids.isEmpty() ) |
| 71 | + { |
| 72 | + // add to tree and to structure |
| 73 | + index.insertFeature( pntId, pt.boundingBox() ); |
| 74 | + |
| 75 | + AnchorPoint xp; |
| 76 | + xp.x = pt.x(); |
| 77 | + xp.y = pt.y(); |
| 78 | + xp.anchor = -1; |
| 79 | + pnts.append( xp ); |
| 80 | + pntId++; |
| 81 | + } |
| 82 | + } |
| 83 | + |
| 84 | + ++count; |
| 85 | + feedback->setProgress( 100. * count / totalCount ); |
| 86 | + } |
| 87 | +} |
| 88 | + |
| 89 | + |
| 90 | +static void assignAnchors( QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh ) |
| 91 | +{ |
| 92 | + double thresh2 = thresh * thresh; |
| 93 | + int nanchors = 0, ntosnap = 0; |
| 94 | + for ( int point = 0; point < pnts.count(); ++point ) |
| 95 | + { |
| 96 | + if ( pnts[point].anchor >= 0 ) |
| 97 | + continue; |
| 98 | + |
| 99 | + pnts[point].anchor = -2; // make it anchor |
| 100 | + nanchors++; |
| 101 | + |
| 102 | + // Find points in threshold |
| 103 | + double x = pnts[point].x, y = pnts[point].y; |
| 104 | + QgsRectangle rect( x - thresh, y - thresh, x + thresh, y + thresh ); |
| 105 | + |
| 106 | + const QList<QgsFeatureId> ids = index.intersects( rect ); |
| 107 | + for ( QgsFeatureId pointb : ids ) |
| 108 | + { |
| 109 | + if ( pointb == point ) |
| 110 | + continue; |
| 111 | + |
| 112 | + double dx = pnts[pointb].x - pnts[point].x; |
| 113 | + double dy = pnts[pointb].y - pnts[point].y; |
| 114 | + double dist2 = dx * dx + dy * dy; |
| 115 | + if ( dist2 > thresh2 ) |
| 116 | + continue; // outside threshold |
| 117 | + |
| 118 | + if ( pnts[pointb].anchor == -1 ) |
| 119 | + { |
| 120 | + // doesn't have an anchor yet |
| 121 | + pnts[pointb].anchor = point; |
| 122 | + ntosnap++; |
| 123 | + } |
| 124 | + else if ( pnts[pointb].anchor >= 0 ) |
| 125 | + { |
| 126 | + // check distance to previously assigned anchor |
| 127 | + double dx2 = pnts[pnts[pointb].anchor].x - pnts[pointb].x; |
| 128 | + double dy2 = pnts[pnts[pointb].anchor].y - pnts[pointb].y; |
| 129 | + double dist2_a = dx2 * dx2 + dy2 * dy2; |
| 130 | + if ( dist2 < dist2_a ) |
| 131 | + pnts[pointb].anchor = point; // replace old anchor |
| 132 | + } |
| 133 | + } |
| 134 | + } |
| 135 | +} |
| 136 | + |
| 137 | + |
| 138 | +static bool snapPoint( QgsPoint *pt, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts ) |
| 139 | +{ |
| 140 | + // Find point ( should always find one point ) |
| 141 | + QList<QgsFeatureId> fids = index.intersects( QgsRectangle( pt->x(), pt->y(), pt->x(), pt->y() ) ); |
| 142 | + Q_ASSERT( fids.count() == 1 ); |
| 143 | + |
| 144 | + int spoint = fids[0]; |
| 145 | + int anchor = pnts[spoint].anchor; |
| 146 | + |
| 147 | + if ( anchor >= 0 ) |
| 148 | + { |
| 149 | + // to be snapped |
| 150 | + pt->setX( pnts[anchor].x ); |
| 151 | + pt->setY( pnts[anchor].y ); |
| 152 | + return true; |
| 153 | + } |
| 154 | + |
| 155 | + return false; |
| 156 | +} |
| 157 | + |
| 158 | + |
| 159 | +static bool snapLineString( QgsLineString *linestring, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh ) |
| 160 | +{ |
| 161 | + QVector<QgsPoint> newPoints; |
| 162 | + QVector<int> anchors; // indexes of anchors for vertices |
| 163 | + double thresh2 = thresh * thresh; |
| 164 | + double minDistX, minDistY; // coordinates of the closest point on the segment line |
| 165 | + bool changed = false; |
| 166 | + |
| 167 | + // snap vertices |
| 168 | + for ( int v = 0; v < linestring->numPoints(); v++ ) |
| 169 | + { |
| 170 | + double x = linestring->xAt( v ); |
| 171 | + double y = linestring->yAt( v ); |
| 172 | + QgsRectangle rect( x, y, x, y ); |
| 173 | + |
| 174 | + // Find point ( should always find one point ) |
| 175 | + QList<QgsFeatureId> fids = index.intersects( rect ); |
| 176 | + Q_ASSERT( fids.count() == 1 ); |
| 177 | + |
| 178 | + int spoint = fids.first(); |
| 179 | + int anchor = pnts[spoint].anchor; |
| 180 | + if ( anchor >= 0 ) |
| 181 | + { |
| 182 | + // to be snapped |
| 183 | + linestring->setXAt( v, pnts[anchor].x ); |
| 184 | + linestring->setYAt( v, pnts[anchor].y ); |
| 185 | + anchors.append( anchor ); // point on new location |
| 186 | + changed = true; |
| 187 | + } |
| 188 | + else |
| 189 | + { |
| 190 | + anchors.append( spoint ); // old point |
| 191 | + } |
| 192 | + } |
| 193 | + |
| 194 | + // Snap all segments to anchors in threshold |
| 195 | + for ( int v = 0; v < linestring->numPoints() - 1; v++ ) |
| 196 | + { |
| 197 | + double x1 = linestring->xAt( v ); |
| 198 | + double x2 = linestring->xAt( v + 1 ); |
| 199 | + double y1 = linestring->yAt( v ); |
| 200 | + double y2 = linestring->yAt( v + 1 ); |
| 201 | + |
| 202 | + newPoints << linestring->pointN( v ); |
| 203 | + |
| 204 | + // Box |
| 205 | + double xmin = x1, xmax = x2, ymin = y1, ymax = y2; |
| 206 | + if ( xmin > xmax ) |
| 207 | + std::swap( xmin, xmax ); |
| 208 | + if ( ymin > ymax ) |
| 209 | + std::swap( ymin, ymax ); |
| 210 | + |
| 211 | + QgsRectangle rect( xmin - thresh, ymin - thresh, xmax + thresh, ymax + thresh ); |
| 212 | + |
| 213 | + // Find points |
| 214 | + const QList<QgsFeatureId> fids = index.intersects( rect ); |
| 215 | + |
| 216 | + QVector<AnchorAlongSegment> newVerticesAlongSegment; |
| 217 | + |
| 218 | + // Snap to anchor in threshold different from end points |
| 219 | + for ( QgsFeatureId fid : fids ) |
| 220 | + { |
| 221 | + int spoint = fid; |
| 222 | + |
| 223 | + if ( spoint == anchors[v] || spoint == anchors[v + 1] ) |
| 224 | + continue; // end point |
| 225 | + if ( pnts[spoint].anchor >= 0 ) |
| 226 | + continue; // point is not anchor |
| 227 | + |
| 228 | + // Check the distance |
| 229 | + double dist2 = QgsGeometryUtils::sqrDistToLine( pnts[spoint].x, pnts[spoint].y, x1, y1, x2, y2, minDistX, minDistY, 0 ); |
| 230 | + // skip points that are behind segment's endpoints or extremely close to them |
| 231 | + double dx1 = minDistX - x1, dx2 = minDistX - x2; |
| 232 | + double dy1 = minDistY - y1, dy2 = minDistY - y2; |
| 233 | + bool isOnSegment = !qgsDoubleNear( dx1 * dx1 + dy1 * dy1, 0 ) && !qgsDoubleNear( dx2 * dx2 + dy2 * dy2, 0 ); |
| 234 | + if ( isOnSegment && dist2 <= thresh2 ) |
| 235 | + { |
| 236 | + // an anchor is in the threshold |
| 237 | + AnchorAlongSegment item; |
| 238 | + item.anchor = spoint; |
| 239 | + item.along = QgsPointXY( x1, y1 ).distance( minDistX, minDistY ); |
| 240 | + newVerticesAlongSegment << item; |
| 241 | + } |
| 242 | + } |
| 243 | + |
| 244 | + if ( !newVerticesAlongSegment.isEmpty() ) |
| 245 | + { |
| 246 | + // sort by distance along the segment |
| 247 | + std::sort( newVerticesAlongSegment.begin(), newVerticesAlongSegment.end(), []( const AnchorAlongSegment & p1, const AnchorAlongSegment & p2 ) |
| 248 | + { |
| 249 | + return ( p1.along < p2.along ? -1 : ( p1.along > p2.along ) ); |
| 250 | + } ); |
| 251 | + |
| 252 | + // insert new vertices |
| 253 | + for ( int i = 0; i < newVerticesAlongSegment.count(); i++ ) |
| 254 | + { |
| 255 | + int anchor = newVerticesAlongSegment[i].anchor; |
| 256 | + newPoints << QgsPoint( pnts[anchor].x, pnts[anchor].y, 0 ); |
| 257 | + } |
| 258 | + changed = true; |
| 259 | + } |
| 260 | + } |
| 261 | + |
| 262 | + // append end point |
| 263 | + newPoints << linestring->pointN( linestring->numPoints() - 1 ); |
| 264 | + |
| 265 | + // replace linestring's points |
| 266 | + if ( changed ) |
| 267 | + linestring->setPoints( newPoints ); |
| 268 | + |
| 269 | + return changed; |
| 270 | +} |
| 271 | + |
| 272 | + |
| 273 | +static bool snapGeometry( QgsAbstractGeometry *g, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh ) |
| 274 | +{ |
| 275 | + bool changed = false; |
| 276 | + if ( QgsLineString *linestring = qgsgeometry_cast<QgsLineString *>( g ) ) |
| 277 | + { |
| 278 | + changed |= snapLineString( linestring, index, pnts, thresh ); |
| 279 | + } |
| 280 | + else if ( QgsPolygon *polygon = qgsgeometry_cast<QgsPolygon *>( g ) ) |
| 281 | + { |
| 282 | + if ( QgsLineString *exteriorRing = qgsgeometry_cast<QgsLineString *>( polygon->exteriorRing() ) ) |
| 283 | + changed |= snapLineString( exteriorRing, index, pnts, thresh ); |
| 284 | + for ( int i = 0; i < polygon->numInteriorRings(); ++i ) |
| 285 | + { |
| 286 | + if ( QgsLineString *interiorRing = qgsgeometry_cast<QgsLineString *>( polygon->interiorRing( i ) ) ) |
| 287 | + changed |= snapLineString( interiorRing, index, pnts, thresh ); |
| 288 | + } |
| 289 | + } |
| 290 | + else if ( QgsGeometryCollection *collection = qgsgeometry_cast<QgsGeometryCollection *>( g ) ) |
| 291 | + { |
| 292 | + for ( int i = 0; i < collection->numGeometries(); ++i ) |
| 293 | + changed |= snapGeometry( collection->geometryN( i ), index, pnts, thresh ); |
| 294 | + } |
| 295 | + else if ( QgsPoint *pt = qgsgeometry_cast<QgsPoint *>( g ) ) |
| 296 | + { |
| 297 | + changed |= snapPoint( pt, index, pnts ); |
| 298 | + } |
| 299 | + |
| 300 | + return changed; |
| 301 | +} |
| 302 | + |
| 303 | + |
| 304 | +int QgsGeometrySnapperSingleSource::run( const QgsFeatureSource &source, QgsFeatureSink &sink, double thresh, QgsFeedback *feedback ) |
| 305 | +{ |
| 306 | + // the logic here comes from GRASS implementation of Vect_snap_lines_list() |
| 307 | + |
| 308 | + int count = 0; |
| 309 | + int totalCount = source.featureCount() * 2; |
| 310 | + |
| 311 | + // step 1: record all point locations in a spatial index + extra data structure to keep |
| 312 | + // reference to which other point they have been snapped to (in the next phase). |
| 313 | + |
| 314 | + QgsSpatialIndex index; |
| 315 | + QVector<AnchorPoint> pnts; |
| 316 | + QgsFeatureRequest request; |
| 317 | + request.setSubsetOfAttributes( QgsAttributeList() ); |
| 318 | + QgsFeatureIterator fi = source.getFeatures( request ); |
| 319 | + buildSnapIndex( fi, index, pnts, feedback, count, totalCount ); |
| 320 | + |
| 321 | + if ( feedback->isCanceled() ) |
| 322 | + return 0; |
| 323 | + |
| 324 | + // step 2: go through all registered points and if not yet marked mark it as anchor and |
| 325 | + // assign this anchor to all not yet marked points in threshold |
| 326 | + |
| 327 | + assignAnchors( index, pnts, thresh ); |
| 328 | + |
| 329 | + // step 3: alignment of vertices and segments to the anchors |
| 330 | + // Go through all lines and: |
| 331 | + // 1) for all vertices: if not anchor snap it to its anchor |
| 332 | + // 2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course) |
| 333 | + |
| 334 | + int modified = 0; |
| 335 | + QgsFeature f; |
| 336 | + fi = source.getFeatures(); |
| 337 | + while ( fi.nextFeature( f ) ) |
| 338 | + { |
| 339 | + if ( feedback->isCanceled() ) |
| 340 | + break; |
| 341 | + |
| 342 | + QgsGeometry geom = f.geometry(); |
| 343 | + if ( snapGeometry( geom.get(), index, pnts, thresh ) ) |
| 344 | + { |
| 345 | + f.setGeometry( geom ); |
| 346 | + ++modified; |
| 347 | + } |
| 348 | + |
| 349 | + sink.addFeature( f, QgsFeatureSink::FastInsert ); |
| 350 | + |
| 351 | + ++count; |
| 352 | + feedback->setProgress( 100. * count / totalCount ); |
| 353 | + } |
| 354 | + |
| 355 | + return modified; |
| 356 | +} |
0 commit comments