QGIS API Documentation  3.37.0-Master (a5b4d9743e8)
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"
28 #include "qgsgeometrycollection.h"
29 #include "qgsvertexid.h"
30 
31 QgsMapToPixelSimplifier::QgsMapToPixelSimplifier( int simplifyFlags, double tolerance, SimplifyAlgorithm simplifyAlgorithm )
32  : mSimplifyFlags( simplifyFlags )
33  , mSimplifyAlgorithm( simplifyAlgorithm )
34  , mTolerance( tolerance )
35 {
36 }
37 
39 // Helper simplification methods
40 
41 float 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 
49 bool 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 
72 static 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 
122 std::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 
261  case SnappedToGridGlobal:
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 
412 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle &envelope, double map2pixelTol )
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:182
@ LineString
LineString.
@ Polygon
Polygon.
@ CircularString
CircularString.
Abstract base class for all geometries.
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.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
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
Returns the WKB type of the geometry.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
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.
Geometry collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from 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.
int numGeometries() const
Returns the number of geometries within the collection.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:162
Q_GADGET bool isNull
Definition: qgsgeometry.h:164
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.
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:33
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:201
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:211
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:236
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:196
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:206
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:243
static bool isMultiType(Qgis::WkbType type)
Returns true if the WKB type is a multi type.
Definition: qgswkbtypes.h:758
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.
Definition: qgswkbtypes.h:628
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:5172
Utility class for identifying a unique vertex within a geometry.
Definition: qgsvertexid.h:30