QGIS API Documentation 3.39.0-Master (bca3cdb6021)
Loading...
Searching...
No Matches
qgsmaptopixelgeometrysimplifier.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsmaptopixelgeometrysimplifier.cpp
3 ---------------------
4 begin : December 2013
5 copyright : (C) 2013 by Alvaro Huarte
6 email : http://wiki.osgeo.org/wiki/Alvaro_Huarte
7
8 ***************************************************************************
9 * *
10 * This program is free software; you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 ***************************************************************************/
16
17#include <limits>
18#include <memory>
19
21#include "qgsrectangle.h"
22#include "qgsgeometry.h"
23#include "qgslinestring.h"
24#include "qgspolygon.h"
26#include "qgsvertexid.h"
27
28QgsMapToPixelSimplifier::QgsMapToPixelSimplifier( int simplifyFlags, double tolerance, Qgis::VectorSimplificationAlgorithm simplifyAlgorithm )
29 : mSimplifyFlags( simplifyFlags )
30 , mSimplifyAlgorithm( simplifyAlgorithm )
31 , mTolerance( tolerance )
32{
33}
34
36// Helper simplification methods
37
38float QgsMapToPixelSimplifier::calculateLengthSquared2D( double x1, double y1, double x2, double y2 )
39{
40 const float vx = static_cast< float >( x2 - x1 );
41 const float vy = static_cast< float >( y2 - y1 );
42
43 return ( vx * vx ) + ( vy * vy );
44}
45
46bool QgsMapToPixelSimplifier::equalSnapToGrid( double x1, double y1, double x2, double y2, double gridOriginX, double gridOriginY, float gridInverseSizeXY )
47{
48 const int grid_x1 = std::round( ( x1 - gridOriginX ) * gridInverseSizeXY );
49 const int grid_x2 = std::round( ( x2 - gridOriginX ) * gridInverseSizeXY );
50 if ( grid_x1 != grid_x2 ) return false;
51
52 const int grid_y1 = std::round( ( y1 - gridOriginY ) * gridInverseSizeXY );
53 const int grid_y2 = std::round( ( y2 - gridOriginY ) * gridInverseSizeXY );
54 return grid_y1 == grid_y2;
55}
56
58// Helper simplification methods for Visvalingam method
59
60// It uses a refactored code of the liblwgeom implementation:
61// https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.h
62// https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.c
63
64#include "simplify/effectivearea.h"
65
67
69static std::unique_ptr< QgsAbstractGeometry > generalizeWkbGeometryByBoundingBox(
70 Qgis::WkbType wkbType,
71 const QgsAbstractGeometry &geometry,
72 const QgsRectangle &envelope,
73 bool isRing )
74{
75 const Qgis::WkbType geometryType = QgsWkbTypes::singleType( QgsWkbTypes::flatType( wkbType ) );
76
77 // If the geometry is already minimal skip the generalization
78 const int minimumSize = geometryType == Qgis::WkbType::LineString ? 2 : 5;
79
80 if ( geometry.nCoordinates() <= minimumSize )
81 {
82 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
83 }
84
85 const double x1 = envelope.xMinimum();
86 const double y1 = envelope.yMinimum();
87 const double x2 = envelope.xMaximum();
88 const double y2 = envelope.yMaximum();
89
90 // Write the generalized geometry
91 if ( geometryType == Qgis::WkbType::LineString && !isRing )
92 {
93 return std::make_unique< QgsLineString >( QVector<double>() << x1 << x2, QVector<double>() << y1 << y2 );
94 }
95 else
96 {
97 std::unique_ptr< QgsLineString > ext = std::make_unique< QgsLineString >(
98 QVector< double >() << x1
99 << x2
100 << x2
101 << x1
102 << x1,
103 QVector< double >() << y1
104 << y1
105 << y2
106 << y2
107 << y1 );
108 if ( geometryType == Qgis::WkbType::LineString )
109 return std::move( ext );
110 else
111 {
112 std::unique_ptr< QgsPolygon > polygon = std::make_unique< QgsPolygon >();
113 polygon->setExteriorRing( ext.release() );
114 return std::move( polygon );
115 }
116 }
117}
118
119std::unique_ptr< QgsAbstractGeometry > QgsMapToPixelSimplifier::simplifyGeometry( int simplifyFlags,
120 Qgis::VectorSimplificationAlgorithm simplifyAlgorithm,
121 const QgsAbstractGeometry &geometry, double map2pixelTol,
122 bool isaLinearRing )
123{
124 bool isGeneralizable = true;
125 const Qgis::WkbType wkbType = geometry.wkbType();
126
127 // Can replace the geometry by its BBOX ?
128 const QgsRectangle envelope = geometry.boundingBox();
130 isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
131 {
132 return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope, isaLinearRing );
133 }
134
136 isGeneralizable = false;
137
138 const Qgis::WkbType flatType = QgsWkbTypes::flatType( wkbType );
139
140 // Write the geometry
141 if ( flatType == Qgis::WkbType::LineString || flatType == Qgis::WkbType::CircularString )
142 {
143 const QgsCurve &srcCurve = dynamic_cast<const QgsCurve &>( geometry );
144 const int numPoints = srcCurve.numPoints();
145
146 std::unique_ptr<QgsCurve> output;
147
148 QVector< double > lineStringX;
149 QVector< double > lineStringY;
150 QVector< double > lineStringZ;
151 QVector< double > lineStringM;
152 if ( flatType == Qgis::WkbType::LineString )
153 {
154 // if we are making a linestring, we do it in an optimised way by directly constructing
155 // the final x/y vectors, which avoids calling the slower insertVertex method
156 lineStringX.reserve( numPoints );
157 lineStringY.reserve( numPoints );
158
159 if ( geometry.is3D() )
160 lineStringZ.reserve( numPoints );
161
162 if ( geometry.isMeasure() )
163 lineStringM.reserve( numPoints );
164 }
165 else
166 {
167 output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.createEmptyWithSameType() ) );
168 }
169
170 double x = 0.0, y = 0.0, z = 0.0, m = 0.0, lastX = 0.0, lastY = 0.0;
171
172 if ( numPoints <= ( isaLinearRing ? 4 : 2 ) )
173 isGeneralizable = false;
174
175 bool isLongSegment;
176 bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
177 const bool is3D = geometry.is3D();
178 const bool isMeasure = geometry.isMeasure();
179
180 // Check whether the LinearRing is really closed.
181 if ( isaLinearRing )
182 {
183 isaLinearRing = qgsDoubleNear( srcCurve.xAt( 0 ), srcCurve.xAt( numPoints - 1 ) ) &&
184 qgsDoubleNear( srcCurve.yAt( 0 ), srcCurve.yAt( numPoints - 1 ) );
185 }
186
187 // Process each vertex...
188 switch ( simplifyAlgorithm )
189 {
191 {
192 const double gridOriginX = envelope.xMinimum();
193 const double gridOriginY = envelope.yMinimum();
194
195 // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
196 const float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
197
198 const double *xData = nullptr;
199 const double *yData = nullptr;
200 const double *zData = nullptr;
201 const double *mData = nullptr;
202 if ( flatType == Qgis::WkbType::LineString )
203 {
204 xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
205 yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
206
207 if ( is3D )
208 zData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->zData();
209
210 if ( isMeasure )
211 mData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->mData();
212 }
213
214 for ( int i = 0; i < numPoints; ++i )
215 {
216 if ( xData && yData )
217 {
218 x = *xData++;
219 y = *yData++;
220 }
221 else
222 {
223 x = srcCurve.xAt( i );
224 y = srcCurve.yAt( i );
225 }
226
227 if ( is3D )
228 z = zData ? *zData++ : srcCurve.zAt( i );
229
230 if ( isMeasure )
231 m = mData ? *mData++ : srcCurve.mAt( i );
232
233 if ( i == 0 ||
234 !isGeneralizable ||
235 !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
236 ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
237 {
238 if ( output )
239 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y, z, m ) );
240 else
241 {
242 lineStringX.append( x );
243 lineStringY.append( y );
244
245 if ( is3D )
246 lineStringZ.append( z );
247
248 if ( isMeasure )
249 lineStringM.append( m );
250 }
251 lastX = x;
252 lastY = y;
253 }
254 }
255 break;
256 }
257
259 {
260 output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.snappedToGrid( map2pixelTol, map2pixelTol ) ) );
261 break;
262 }
263
265 {
266 map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
267
268 EFFECTIVE_AREAS ea( srcCurve );
269
270 const int set_area = 0;
271 ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
272
273 for ( int i = 0; i < numPoints; ++i )
274 {
275 if ( ea.res_arealist[ i ] > map2pixelTol )
276 {
277 if ( output )
278 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
279 else
280 {
281 lineStringX.append( ea.inpts.at( i ).x() );
282 lineStringY.append( ea.inpts.at( i ).y() );
283
284 if ( is3D )
285 lineStringZ.append( ea.inpts.at( i ).z() );
286
287 if ( isMeasure )
288 lineStringM.append( ea.inpts.at( i ).m() );
289 }
290 }
291 }
292 break;
293 }
294
296 {
297 map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
298
299 const double *xData = nullptr;
300 const double *yData = nullptr;
301 if ( flatType == Qgis::WkbType::LineString )
302 {
303 xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
304 yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
305 }
306
307 for ( int i = 0; i < numPoints; ++i )
308 {
309 if ( xData && yData )
310 {
311 x = *xData++;
312 y = *yData++;
313 }
314 else
315 {
316 x = srcCurve.xAt( i );
317 y = srcCurve.yAt( i );
318 }
319
320 isLongSegment = false;
321
322 if ( i == 0 ||
323 !isGeneralizable ||
324 ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
325 ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
326 {
327 if ( output )
328 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
329 else
330 {
331 lineStringX.append( x );
332 lineStringY.append( y );
333 }
334 lastX = x;
335 lastY = y;
336
337 hasLongSegments |= isLongSegment;
338 }
339 }
340 }
341 }
342
343 if ( !output )
344 {
345 output = std::make_unique< QgsLineString >( lineStringX, lineStringY, lineStringZ, lineStringM );
346 }
347 if ( output->numPoints() < ( isaLinearRing ? 4 : 2 ) )
348 {
349 // we simplified the geometry too much!
350 if ( !hasLongSegments )
351 {
352 // approximate the geometry's shape by its bounding box
353 // (rect for linear ring / one segment for line string)
354 return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope, isaLinearRing );
355 }
356 else
357 {
358 // Bad luck! The simplified geometry is invalid and approximation by bounding box
359 // would create artifacts due to long segments.
360 // We will return the original geometry
361 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
362 }
363 }
364
365 if ( isaLinearRing )
366 {
367 // make sure we keep the linear ring closed
368 if ( !qgsDoubleNear( lastX, output->xAt( 0 ) ) || !qgsDoubleNear( lastY, output->yAt( 0 ) ) )
369 {
370 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( output->xAt( 0 ), output->yAt( 0 ) ) );
371 }
372 }
373
374 return std::move( output );
375 }
376 else if ( flatType == Qgis::WkbType::Polygon )
377 {
378 const QgsPolygon &srcPolygon = dynamic_cast<const QgsPolygon &>( geometry );
379 std::unique_ptr<QgsPolygon> polygon( new QgsPolygon() );
380 std::unique_ptr<QgsAbstractGeometry> extRing = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *srcPolygon.exteriorRing(), map2pixelTol, true );
381 polygon->setExteriorRing( qgsgeometry_cast<QgsCurve *>( extRing.release() ) );
382 for ( int i = 0; i < srcPolygon.numInteriorRings(); ++i )
383 {
384 const QgsCurve *sub = srcPolygon.interiorRing( i );
385 std::unique_ptr< QgsAbstractGeometry > ring = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, true );
386 polygon->addInteriorRing( qgsgeometry_cast<QgsCurve *>( ring.release() ) );
387 }
388 return std::move( polygon );
389 }
390 else if ( QgsWkbTypes::isMultiType( flatType ) )
391 {
392 const QgsGeometryCollection &srcCollection = dynamic_cast<const QgsGeometryCollection &>( geometry );
393 std::unique_ptr<QgsGeometryCollection> collection( srcCollection.createEmptyWithSameType() );
394 const int numGeoms = srcCollection.numGeometries();
395 collection->reserve( numGeoms );
396 for ( int i = 0; i < numGeoms; ++i )
397 {
398 const QgsAbstractGeometry *sub = srcCollection.geometryN( i );
399 std::unique_ptr< QgsAbstractGeometry > part = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, false );
400 collection->addGeometry( part.release() );
401 }
402 return std::move( collection );
403 }
404 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
405}
406
408
410{
411 // Can replace the geometry by its BBOX ?
412 return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
413}
414
416{
417 if ( geometry.isNull() )
418 {
419 return QgsGeometry();
420 }
422 {
423 return geometry;
424 }
425
426 // Check whether the geometry can be simplified using the map2pixel context
427 const Qgis::WkbType singleType = QgsWkbTypes::singleType( geometry.wkbType() );
428 const Qgis::WkbType flatType = QgsWkbTypes::flatType( singleType );
429 if ( flatType == Qgis::WkbType::Point )
430 {
431 return geometry;
432 }
433
434 const bool isaLinearRing = flatType == Qgis::WkbType::Polygon;
435 const int numPoints = geometry.constGet()->nCoordinates();
436
437 if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
438 {
439 // No simplify simple geometries
440 return geometry;
441 }
442
443 const QgsRectangle envelope = geometry.boundingBox();
444 if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
445 {
446 //points are in average too far apart to lead to any significant simplification
447 return geometry;
448 }
449
450 return QgsGeometry( simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry.constGet(), mTolerance, false ) );
451}
452
454{
455 //
456 // IMPORTANT!!!!!!!
457 // We want to avoid any geometry cloning we possibly can here, which is why the
458 // "fail" paths always return nullptr
459 //
460
461 if ( !geometry )
462 {
463 return nullptr;
464 }
466 {
467 return nullptr;
468 }
469
470 // Check whether the geometry can be simplified using the map2pixel context
471 const Qgis::WkbType singleType = QgsWkbTypes::singleType( geometry->wkbType() );
472 const Qgis::WkbType flatType = QgsWkbTypes::flatType( singleType );
473 if ( flatType == Qgis::WkbType::Point )
474 {
475 return nullptr;
476 }
477
478 const bool isaLinearRing = flatType == Qgis::WkbType::Polygon;
479 const int numPoints = geometry->nCoordinates();
480
481 if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
482 {
483 // No simplify simple geometries
484 return nullptr;
485 }
486
487 const QgsRectangle envelope = geometry->boundingBox();
488 if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
489 {
490 //points are in average too far apart to lead to any significant simplification
491 return nullptr;
492 }
493
494 return simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry, mTolerance, false ).release();
495}
VectorSimplificationAlgorithm
Simplification algorithms for vector features.
Definition qgis.h:2781
@ Distance
The simplification uses the distance between points to remove duplicate points.
@ SnappedToGridGlobal
Snap to a global grid based on the tolerance. Good for consistent results for incoming vertices,...
@ SnapToGrid
The simplification uses a grid (similar to ST_SnapToGrid) to remove duplicate points.
@ Visvalingam
The simplification gives each point in a line an importance weighting, so that least important points...
WkbType
The WKB type describes the number of dimensions a geometry has.
Definition qgis.h:256
@ LineString
LineString.
@ Polygon
Polygon.
@ CircularString
CircularString.
Abstract base class for all geometries.
virtual QgsAbstractGeometry * snappedToGrid(double hSpacing, double vSpacing, double dSpacing=0, double mSpacing=0, bool removeRedundantPoints=false) const =0
Makes a new geometry with all the points or vertices snapped to the closest point of the grid.
bool isMeasure() const
Returns true if the geometry contains m values.
virtual QgsRectangle boundingBox() const
Returns the minimal bounding box for the geometry.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
virtual QgsAbstractGeometry * createEmptyWithSameType() const =0
Creates a new geometry with the same class and same WKB type as the original and transfers ownership.
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
Qgis::WkbType wkbType() const
Returns the WKB type of the geometry.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
Abstract base class for curved geometry type.
Definition qgscurve.h:35
virtual int numPoints() const =0
Returns the number of points in the curve.
virtual double xAt(int index) const =0
Returns the x-coordinate of the specified node in the line string.
virtual double zAt(int index) const =0
Returns the z-coordinate of the specified node in the line string.
virtual double mAt(int index) const =0
Returns the m-coordinate of the specified node in the line string.
virtual double yAt(int index) const =0
Returns the y-coordinate of the specified node in the line string.
QgsGeometryCollection * createEmptyWithSameType() const override
Creates a new geometry with the same class and same WKB type as the original and transfers ownership.
int numGeometries() const
Returns the number of geometries within the collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
A geometry is the spatial representation of a feature.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
Qgis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
double mTolerance
Distance tolerance for the simplification.
static float calculateLengthSquared2D(double x1, double y1, double x2, double y2)
Returns the squared 2D-distance of the vector defined by the two points specified.
static bool isGeneralizableByMapBoundingBox(const QgsRectangle &envelope, double map2pixelTol)
Returns whether the envelope can be replaced by its BBOX when is applied the specified map2pixel cont...
int simplifyFlags() const
Gets the simplification hints of the vector layer managed.
Qgis::VectorSimplificationAlgorithm mSimplifyAlgorithm
Current algorithm.
@ NoFlags
No simplification can be applied.
@ SimplifyEnvelope
The geometries can be fully simplified by its BoundingBox.
@ SimplifyGeometry
The geometries can be simplified using the current map2pixel context state.
int mSimplifyFlags
Current simplification flags.
static bool equalSnapToGrid(double x1, double y1, double x2, double y2, double gridOriginX, double gridOriginY, float gridInverseSizeXY)
Returns whether the points belong to the same grid.
QgsGeometry simplify(const QgsGeometry &geometry) const override
Returns a simplified version the specified geometry.
Qgis::VectorSimplificationAlgorithm simplifyAlgorithm() const
Gets the local simplification algorithm of the vector layer managed.
QgsMapToPixelSimplifier(int simplifyFlags, double tolerance, Qgis::VectorSimplificationAlgorithm simplifyAlgorithm=Qgis::VectorSimplificationAlgorithm::Distance)
Constructor.
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
Polygon geometry type.
Definition qgspolygon.h:33
A rectangle specified with double values.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
double width() const
Returns the width of the rectangle.
double xMaximum() const
Returns the x maximum value (right side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
double height() const
Returns the height of the rectangle.
static bool isMultiType(Qgis::WkbType type)
Returns true if the WKB type is a multi type.
static Qgis::WkbType singleType(Qgis::WkbType type)
Returns the single type for a WKB type.
Definition qgswkbtypes.h:53
static Qgis::WkbType flatType(Qgis::WkbType type)
Returns the flat type for a WKB type.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition qgis.h:5883
Utility class for identifying a unique vertex within a geometry.
Definition qgsvertexid.h:30