QGIS API Documentation  3.18.1-Zürich (202f1bf7e5)
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 
30 QgsMapToPixelSimplifier::QgsMapToPixelSimplifier( int simplifyFlags, double tolerance, SimplifyAlgorithm simplifyAlgorithm )
31  : mSimplifyFlags( simplifyFlags )
32  , mSimplifyAlgorithm( simplifyAlgorithm )
33  , mTolerance( tolerance )
34 {
35 }
36 
38 // Helper simplification methods
39 
40 float QgsMapToPixelSimplifier::calculateLengthSquared2D( double x1, double y1, double x2, double y2 )
41 {
42  float vx = static_cast< float >( x2 - x1 );
43  float vy = static_cast< float >( y2 - y1 );
44 
45  return ( vx * vx ) + ( vy * vy );
46 }
47 
48 bool QgsMapToPixelSimplifier::equalSnapToGrid( double x1, double y1, double x2, double y2, double gridOriginX, double gridOriginY, float gridInverseSizeXY )
49 {
50  int grid_x1 = std::round( ( x1 - gridOriginX ) * gridInverseSizeXY );
51  int grid_x2 = std::round( ( x2 - gridOriginX ) * gridInverseSizeXY );
52  if ( grid_x1 != grid_x2 ) return false;
53 
54  int grid_y1 = std::round( ( y1 - gridOriginY ) * gridInverseSizeXY );
55  int grid_y2 = std::round( ( y2 - gridOriginY ) * gridInverseSizeXY );
56  return grid_y1 == grid_y2;
57 }
58 
60 // Helper simplification methods for Visvalingam method
61 
62 // It uses a refactored code of the liblwgeom implementation:
63 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.h
64 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.c
65 
66 #include "simplify/effectivearea.h"
67 
69 
71 static std::unique_ptr< QgsAbstractGeometry > generalizeWkbGeometryByBoundingBox(
72  QgsWkbTypes::Type wkbType,
73  const QgsAbstractGeometry &geometry,
74  const QgsRectangle &envelope,
75  bool isRing )
76 {
77  unsigned int geometryType = QgsWkbTypes::singleType( QgsWkbTypes::flatType( wkbType ) );
78 
79  // If the geometry is already minimal skip the generalization
80  int minimumSize = geometryType == QgsWkbTypes::LineString ? 2 : 5;
81 
82  if ( geometry.nCoordinates() <= minimumSize )
83  {
84  return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
85  }
86 
87  const double x1 = envelope.xMinimum();
88  const double y1 = envelope.yMinimum();
89  const double x2 = envelope.xMaximum();
90  const double y2 = envelope.yMaximum();
91 
92  // Write the generalized geometry
93  if ( geometryType == QgsWkbTypes::LineString && !isRing )
94  {
95  return qgis::make_unique< QgsLineString >( QVector<double>() << x1 << x2, QVector<double>() << y1 << y2 );
96  }
97  else
98  {
99  std::unique_ptr< QgsLineString > ext = qgis::make_unique< QgsLineString >(
100  QVector< double >() << x1
101  << x2
102  << x2
103  << x1
104  << x1,
105  QVector< double >() << y1
106  << y1
107  << y2
108  << y2
109  << y1 );
110  if ( geometryType == QgsWkbTypes::LineString )
111  return std::move( ext );
112  else
113  {
114  std::unique_ptr< QgsPolygon > polygon = qgis::make_unique< QgsPolygon >();
115  polygon->setExteriorRing( ext.release() );
116  return std::move( polygon );
117  }
118  }
119 }
120 
121 std::unique_ptr< QgsAbstractGeometry > QgsMapToPixelSimplifier::simplifyGeometry( int simplifyFlags,
122  SimplifyAlgorithm simplifyAlgorithm,
123  const QgsAbstractGeometry &geometry, double map2pixelTol,
124  bool isaLinearRing )
125 {
126  bool isGeneralizable = true;
127  QgsWkbTypes::Type wkbType = geometry.wkbType();
128 
129  // Can replace the geometry by its BBOX ?
130  QgsRectangle envelope = geometry.boundingBox();
132  isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
133  {
134  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope, isaLinearRing );
135  }
136 
138  isGeneralizable = false;
139 
140  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( wkbType );
141 
142  // Write the geometry
143  if ( flatType == QgsWkbTypes::LineString || flatType == QgsWkbTypes::CircularString )
144  {
145  const QgsCurve &srcCurve = dynamic_cast<const QgsCurve &>( geometry );
146  const int numPoints = srcCurve.numPoints();
147 
148  std::unique_ptr<QgsCurve> output;
149 
150  QVector< double > lineStringX;
151  QVector< double > lineStringY;
152  if ( flatType == QgsWkbTypes::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  else
160  {
161  output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.createEmptyWithSameType() ) );
162  }
163 
164  double x = 0.0, y = 0.0, lastX = 0.0, lastY = 0.0;
165 
166  if ( numPoints <= ( isaLinearRing ? 4 : 2 ) )
167  isGeneralizable = false;
168 
169  bool isLongSegment;
170  bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
171 
172  // Check whether the LinearRing is really closed.
173  if ( isaLinearRing )
174  {
175  isaLinearRing = qgsDoubleNear( srcCurve.xAt( 0 ), srcCurve.xAt( numPoints - 1 ) ) &&
176  qgsDoubleNear( srcCurve.yAt( 0 ), srcCurve.yAt( numPoints - 1 ) );
177  }
178 
179  // Process each vertex...
180  switch ( simplifyAlgorithm )
181  {
182  case SnapToGrid:
183  {
184  double gridOriginX = envelope.xMinimum();
185  double gridOriginY = envelope.yMinimum();
186 
187  // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
188  float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
189 
190  const double *xData = nullptr;
191  const double *yData = nullptr;
192  if ( flatType == QgsWkbTypes::LineString )
193  {
194  xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
195  yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
196  }
197 
198  for ( int i = 0; i < numPoints; ++i )
199  {
200  if ( xData && yData )
201  {
202  x = *xData++;
203  y = *yData++;
204  }
205  else
206  {
207  x = srcCurve.xAt( i );
208  y = srcCurve.yAt( i );
209  }
210 
211  if ( i == 0 ||
212  !isGeneralizable ||
213  !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
214  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
215  {
216  if ( output )
217  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
218  else
219  {
220  lineStringX.append( x );
221  lineStringY.append( y );
222  }
223  lastX = x;
224  lastY = y;
225  }
226  }
227  break;
228  }
229 
230  case SnappedToGridGlobal:
231  {
232  output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.snappedToGrid( map2pixelTol, map2pixelTol ) ) );
233  break;
234  }
235 
236  case Visvalingam:
237  {
238  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
239 
240  EFFECTIVE_AREAS ea( srcCurve );
241 
242  int set_area = 0;
243  ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
244 
245  for ( int i = 0; i < numPoints; ++i )
246  {
247  if ( ea.res_arealist[ i ] > map2pixelTol )
248  {
249  if ( output )
250  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
251  else
252  {
253  lineStringX.append( ea.inpts.at( i ).x() );
254  lineStringY.append( ea.inpts.at( i ).y() );
255  }
256  }
257  }
258  break;
259  }
260 
261  case Distance:
262  {
263  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
264 
265  const double *xData = nullptr;
266  const double *yData = nullptr;
267  if ( flatType == QgsWkbTypes::LineString )
268  {
269  xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
270  yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
271  }
272 
273  for ( int i = 0; i < numPoints; ++i )
274  {
275  if ( xData && yData )
276  {
277  x = *xData++;
278  y = *yData++;
279  }
280  else
281  {
282  x = srcCurve.xAt( i );
283  y = srcCurve.yAt( i );
284  }
285 
286  isLongSegment = false;
287 
288  if ( i == 0 ||
289  !isGeneralizable ||
290  ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
291  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
292  {
293  if ( output )
294  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
295  else
296  {
297  lineStringX.append( x );
298  lineStringY.append( y );
299  }
300  lastX = x;
301  lastY = y;
302 
303  hasLongSegments |= isLongSegment;
304  }
305  }
306  }
307  }
308 
309  if ( !output )
310  {
311  output = qgis::make_unique< QgsLineString >( lineStringX, lineStringY );
312  }
313  if ( output->numPoints() < ( isaLinearRing ? 4 : 2 ) )
314  {
315  // we simplified the geometry too much!
316  if ( !hasLongSegments )
317  {
318  // approximate the geometry's shape by its bounding box
319  // (rect for linear ring / one segment for line string)
320  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope, isaLinearRing );
321  }
322  else
323  {
324  // Bad luck! The simplified geometry is invalid and approximation by bounding box
325  // would create artifacts due to long segments.
326  // We will return the original geometry
327  return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
328  }
329  }
330 
331  if ( isaLinearRing )
332  {
333  // make sure we keep the linear ring closed
334  if ( !qgsDoubleNear( lastX, output->xAt( 0 ) ) || !qgsDoubleNear( lastY, output->yAt( 0 ) ) )
335  {
336  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( output->xAt( 0 ), output->yAt( 0 ) ) );
337  }
338  }
339 
340  return std::move( output );
341  }
342  else if ( flatType == QgsWkbTypes::Polygon )
343  {
344  const QgsPolygon &srcPolygon = dynamic_cast<const QgsPolygon &>( geometry );
345  std::unique_ptr<QgsPolygon> polygon( new QgsPolygon() );
346  std::unique_ptr<QgsAbstractGeometry> extRing = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *srcPolygon.exteriorRing(), map2pixelTol, true );
347  polygon->setExteriorRing( qgsgeometry_cast<QgsCurve *>( extRing.release() ) );
348  for ( int i = 0; i < srcPolygon.numInteriorRings(); ++i )
349  {
350  const QgsCurve *sub = srcPolygon.interiorRing( i );
351  std::unique_ptr< QgsAbstractGeometry > ring = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, true );
352  polygon->addInteriorRing( qgsgeometry_cast<QgsCurve *>( ring.release() ) );
353  }
354  return std::move( polygon );
355  }
356  else if ( QgsWkbTypes::isMultiType( flatType ) )
357  {
358  const QgsGeometryCollection &srcCollection = dynamic_cast<const QgsGeometryCollection &>( geometry );
359  std::unique_ptr<QgsGeometryCollection> collection( srcCollection.createEmptyWithSameType() );
360  const int numGeoms = srcCollection.numGeometries();
361  collection->reserve( numGeoms );
362  for ( int i = 0; i < numGeoms; ++i )
363  {
364  const QgsAbstractGeometry *sub = srcCollection.geometryN( i );
365  std::unique_ptr< QgsAbstractGeometry > part = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, false );
366  collection->addGeometry( part.release() );
367  }
368  return std::move( collection );
369  }
370  return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
371 }
372 
374 
375 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle &envelope, double map2pixelTol )
376 {
377  // Can replace the geometry by its BBOX ?
378  return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
379 }
380 
382 {
383  if ( geometry.isNull() )
384  {
385  return QgsGeometry();
386  }
388  {
389  return geometry;
390  }
391 
392  // Check whether the geometry can be simplified using the map2pixel context
393  const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry.wkbType() );
394  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
395  if ( flatType == QgsWkbTypes::Point )
396  {
397  return geometry;
398  }
399 
400  const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
401  const int numPoints = geometry.constGet()->nCoordinates();
402 
403  if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
404  {
405  // No simplify simple geometries
406  return geometry;
407  }
408 
409  const QgsRectangle envelope = geometry.boundingBox();
410  if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
411  {
412  //points are in average too far apart to lead to any significant simplification
413  return geometry;
414  }
415 
416  return QgsGeometry( simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry.constGet(), mTolerance, false ) );
417 }
418 
420 {
421  //
422  // IMPORTANT!!!!!!!
423  // We want to avoid any geometry cloning we possibly can here, which is why the
424  // "fail" paths always return nullptr
425  //
426 
427  if ( !geometry )
428  {
429  return nullptr;
430  }
432  {
433  return nullptr;
434  }
435 
436  // Check whether the geometry can be simplified using the map2pixel context
437  const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry->wkbType() );
438  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
439  if ( flatType == QgsWkbTypes::Point )
440  {
441  return nullptr;
442  }
443 
444  const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
445  const int numPoints = geometry->nCoordinates();
446 
447  if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
448  {
449  // No simplify simple geometries
450  return nullptr;
451  }
452 
453  const QgsRectangle envelope = geometry->boundingBox();
454  if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
455  {
456  //points are in average too far apart to lead to any significant simplification
457  return nullptr;
458  }
459 
460  return simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry, mTolerance, false ).release();
461 }
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:124
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:126
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:38
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:172
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:162
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:167
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:177
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:209
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:202
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:316
Utility class for identifying a unique vertex within a geometry.