QGIS API Documentation  3.23.0-Master (22c16f2067)
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  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 
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 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 
231  case SnappedToGridGlobal:
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 );
367  collection->addGeometry( part.release() );
368  }
369  return std::move( collection );
370  }
371  return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
372 }
373 
375 
376 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle &envelope, double map2pixelTol )
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 * 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.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
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.
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.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:125
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.)
Q_GADGET bool isNull
Definition: qgsgeometry.h:127
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:832
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:702
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:1491
Utility class for identifying a unique vertex within a geometry.
Definition: qgsvertexid.h:31