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 *
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 QgsWkbTypes::Type wkbType,
74 const QgsAbstractGeometry &geometry,
75 const QgsRectangle &envelope,
76 bool isRing )
77{
78 const unsigned int geometryType = QgsWkbTypes::singleType( QgsWkbTypes::flatType( wkbType ) );
79
80 // If the geometry is already minimal skip the generalization
81 const int minimumSize = geometryType == QgsWkbTypes::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 == QgsWkbTypes::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 == QgsWkbTypes::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 QgsWkbTypes::Type 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 QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( wkbType );
142
143 // Write the geometry
144 if ( flatType == QgsWkbTypes::LineString || flatType == QgsWkbTypes::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 if ( flatType == QgsWkbTypes::LineString )
154 {
155 // if we are making a linestring, we do it in an optimised way by directly constructing
156 // the final x/y vectors, which avoids calling the slower insertVertex method
157 lineStringX.reserve( numPoints );
158 lineStringY.reserve( numPoints );
159 }
160 else
161 {
162 output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.createEmptyWithSameType() ) );
163 }
164
165 double x = 0.0, y = 0.0, lastX = 0.0, lastY = 0.0;
166
167 if ( numPoints <= ( isaLinearRing ? 4 : 2 ) )
168 isGeneralizable = false;
169
170 bool isLongSegment;
171 bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
172
173 // Check whether the LinearRing is really closed.
174 if ( isaLinearRing )
175 {
176 isaLinearRing = qgsDoubleNear( srcCurve.xAt( 0 ), srcCurve.xAt( numPoints - 1 ) ) &&
177 qgsDoubleNear( srcCurve.yAt( 0 ), srcCurve.yAt( numPoints - 1 ) );
178 }
179
180 // Process each vertex...
181 switch ( simplifyAlgorithm )
182 {
183 case SnapToGrid:
184 {
185 const double gridOriginX = envelope.xMinimum();
186 const double gridOriginY = envelope.yMinimum();
187
188 // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
189 const float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
190
191 const double *xData = nullptr;
192 const double *yData = nullptr;
193 if ( flatType == QgsWkbTypes::LineString )
194 {
195 xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
196 yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
197 }
198
199 for ( int i = 0; i < numPoints; ++i )
200 {
201 if ( xData && yData )
202 {
203 x = *xData++;
204 y = *yData++;
205 }
206 else
207 {
208 x = srcCurve.xAt( i );
209 y = srcCurve.yAt( i );
210 }
211
212 if ( i == 0 ||
213 !isGeneralizable ||
214 !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
215 ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
216 {
217 if ( output )
218 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
219 else
220 {
221 lineStringX.append( x );
222 lineStringY.append( y );
223 }
224 lastX = x;
225 lastY = y;
226 }
227 }
228 break;
229 }
230
232 {
233 output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.snappedToGrid( map2pixelTol, map2pixelTol ) ) );
234 break;
235 }
236
237 case Visvalingam:
238 {
239 map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
240
241 EFFECTIVE_AREAS ea( srcCurve );
242
243 const int set_area = 0;
244 ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
245
246 for ( int i = 0; i < numPoints; ++i )
247 {
248 if ( ea.res_arealist[ i ] > map2pixelTol )
249 {
250 if ( output )
251 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
252 else
253 {
254 lineStringX.append( ea.inpts.at( i ).x() );
255 lineStringY.append( ea.inpts.at( i ).y() );
256 }
257 }
258 }
259 break;
260 }
261
262 case Distance:
263 {
264 map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
265
266 const double *xData = nullptr;
267 const double *yData = nullptr;
268 if ( flatType == QgsWkbTypes::LineString )
269 {
270 xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
271 yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
272 }
273
274 for ( int i = 0; i < numPoints; ++i )
275 {
276 if ( xData && yData )
277 {
278 x = *xData++;
279 y = *yData++;
280 }
281 else
282 {
283 x = srcCurve.xAt( i );
284 y = srcCurve.yAt( i );
285 }
286
287 isLongSegment = false;
288
289 if ( i == 0 ||
290 !isGeneralizable ||
291 ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
292 ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
293 {
294 if ( output )
295 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
296 else
297 {
298 lineStringX.append( x );
299 lineStringY.append( y );
300 }
301 lastX = x;
302 lastY = y;
303
304 hasLongSegments |= isLongSegment;
305 }
306 }
307 }
308 }
309
310 if ( !output )
311 {
312 output = std::make_unique< QgsLineString >( lineStringX, lineStringY );
313 }
314 if ( output->numPoints() < ( isaLinearRing ? 4 : 2 ) )
315 {
316 // we simplified the geometry too much!
317 if ( !hasLongSegments )
318 {
319 // approximate the geometry's shape by its bounding box
320 // (rect for linear ring / one segment for line string)
321 return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope, isaLinearRing );
322 }
323 else
324 {
325 // Bad luck! The simplified geometry is invalid and approximation by bounding box
326 // would create artifacts due to long segments.
327 // We will return the original geometry
328 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
329 }
330 }
331
332 if ( isaLinearRing )
333 {
334 // make sure we keep the linear ring closed
335 if ( !qgsDoubleNear( lastX, output->xAt( 0 ) ) || !qgsDoubleNear( lastY, output->yAt( 0 ) ) )
336 {
337 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( output->xAt( 0 ), output->yAt( 0 ) ) );
338 }
339 }
340
341 return std::move( output );
342 }
343 else if ( flatType == QgsWkbTypes::Polygon )
344 {
345 const QgsPolygon &srcPolygon = dynamic_cast<const QgsPolygon &>( geometry );
346 std::unique_ptr<QgsPolygon> polygon( new QgsPolygon() );
347 std::unique_ptr<QgsAbstractGeometry> extRing = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *srcPolygon.exteriorRing(), map2pixelTol, true );
348 polygon->setExteriorRing( qgsgeometry_cast<QgsCurve *>( extRing.release() ) );
349 for ( int i = 0; i < srcPolygon.numInteriorRings(); ++i )
350 {
351 const QgsCurve *sub = srcPolygon.interiorRing( i );
352 std::unique_ptr< QgsAbstractGeometry > ring = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, true );
353 polygon->addInteriorRing( qgsgeometry_cast<QgsCurve *>( ring.release() ) );
354 }
355 return std::move( polygon );
356 }
357 else if ( QgsWkbTypes::isMultiType( flatType ) )
358 {
359 const QgsGeometryCollection &srcCollection = dynamic_cast<const QgsGeometryCollection &>( geometry );
360 std::unique_ptr<QgsGeometryCollection> collection( srcCollection.createEmptyWithSameType() );
361 const int numGeoms = srcCollection.numGeometries();
362 collection->reserve( numGeoms );
363 for ( int i = 0; i < numGeoms; ++i )
364 {
365 const QgsAbstractGeometry *sub = srcCollection.geometryN( i );
366 std::unique_ptr< QgsAbstractGeometry > part = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, false );
368 }
369 return std::move( collection );
370 }
371 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
372}
373
375
377{
378 // Can replace the geometry by its BBOX ?
379 return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
380}
381
383{
384 if ( geometry.isNull() )
385 {
386 return QgsGeometry();
387 }
389 {
390 return geometry;
391 }
392
393 // Check whether the geometry can be simplified using the map2pixel context
394 const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry.wkbType() );
395 const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
396 if ( flatType == QgsWkbTypes::Point )
397 {
398 return geometry;
399 }
400
401 const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
402 const int numPoints = geometry.constGet()->nCoordinates();
403
404 if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
405 {
406 // No simplify simple geometries
407 return geometry;
408 }
409
410 const QgsRectangle envelope = geometry.boundingBox();
411 if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
412 {
413 //points are in average too far apart to lead to any significant simplification
414 return geometry;
415 }
416
417 return QgsGeometry( simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry.constGet(), mTolerance, false ) );
418}
419
421{
422 //
423 // IMPORTANT!!!!!!!
424 // We want to avoid any geometry cloning we possibly can here, which is why the
425 // "fail" paths always return nullptr
426 //
427
428 if ( !geometry )
429 {
430 return nullptr;
431 }
433 {
434 return nullptr;
435 }
436
437 // Check whether the geometry can be simplified using the map2pixel context
438 const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry->wkbType() );
439 const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
440 if ( flatType == QgsWkbTypes::Point )
441 {
442 return nullptr;
443 }
444
445 const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
446 const int numPoints = geometry->nCoordinates();
447
448 if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
449 {
450 // No simplify simple geometries
451 return nullptr;
452 }
453
454 const QgsRectangle envelope = geometry->boundingBox();
455 if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
456 {
457 //points are in average too far apart to lead to any significant simplification
458 return nullptr;
459 }
460
461 return simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry, mTolerance, false ).release();
462}
Abstract base class for all geometries.
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.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
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 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.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
Definition: qgsgeometry.h:166
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(Type type) SIP_HOLDGIL
Returns true if the WKB type is a multi type.
Definition: qgswkbtypes.h:862
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:70
static Type singleType(Type type) SIP_HOLDGIL
Returns the single type for a WKB type.
Definition: qgswkbtypes.h:157
static Type flatType(Type type) SIP_HOLDGIL
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:732
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:2527
Utility class for identifying a unique vertex within a geometry.
Definition: qgsvertexid.h:31