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