QGIS API Documentation  3.2.0-Bonn (bc43194)
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,
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  QgsRectangle r;
166  r.setMinimal();
167 
168  if ( numPoints <= ( isaLinearRing ? 4 : 2 ) )
169  isGeneralizable = false;
170 
171  bool isLongSegment;
172  bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
173 
174  // Check whether the LinearRing is really closed.
175  if ( isaLinearRing )
176  {
177  isaLinearRing = qgsDoubleNear( srcCurve.xAt( 0 ), srcCurve.xAt( numPoints - 1 ) ) &&
178  qgsDoubleNear( srcCurve.yAt( 0 ), srcCurve.yAt( numPoints - 1 ) );
179  }
180 
181  // Process each vertex...
182  switch ( simplifyAlgorithm )
183  {
184  case SnapToGrid:
185  {
186  double gridOriginX = envelope.xMinimum();
187  double gridOriginY = envelope.yMinimum();
188 
189  // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
190  float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
191 
192  const double *xData = nullptr;
193  const double *yData = nullptr;
194  if ( flatType == QgsWkbTypes::LineString )
195  {
196  xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
197  yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
198  }
199 
200  for ( int i = 0; i < numPoints; ++i )
201  {
202  if ( xData && yData )
203  {
204  x = *xData++;
205  y = *yData++;
206  }
207  else
208  {
209  x = srcCurve.xAt( i );
210  y = srcCurve.yAt( i );
211  }
212 
213  if ( i == 0 ||
214  !isGeneralizable ||
215  !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
216  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
217  {
218  if ( output )
219  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
220  else
221  {
222  lineStringX.append( x );
223  lineStringY.append( y );
224  }
225  lastX = x;
226  lastY = y;
227  }
228 
229  r.combineExtentWith( x, y );
230  }
231  break;
232  }
233 
234  case Visvalingam:
235  {
236  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
237 
238  EFFECTIVE_AREAS ea( srcCurve );
239 
240  int set_area = 0;
241  ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
242 
243  for ( int i = 0; i < numPoints; ++i )
244  {
245  if ( ea.res_arealist[ i ] > map2pixelTol )
246  {
247  if ( output )
248  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
249  else
250  {
251  lineStringX.append( ea.inpts.at( i ).x() );
252  lineStringY.append( ea.inpts.at( i ).y() );
253  }
254  }
255  }
256  break;
257  }
258 
259  case Distance:
260  {
261  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
262 
263  const double *xData = nullptr;
264  const double *yData = nullptr;
265  if ( flatType == QgsWkbTypes::LineString )
266  {
267  xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
268  yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
269  }
270 
271  for ( int i = 0; i < numPoints; ++i )
272  {
273  if ( xData && yData )
274  {
275  x = *xData++;
276  y = *yData++;
277  }
278  else
279  {
280  x = srcCurve.xAt( i );
281  y = srcCurve.yAt( i );
282  }
283 
284  isLongSegment = false;
285 
286  if ( i == 0 ||
287  !isGeneralizable ||
288  ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
289  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
290  {
291  if ( output )
292  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
293  else
294  {
295  lineStringX.append( x );
296  lineStringY.append( y );
297  }
298  lastX = x;
299  lastY = y;
300 
301  hasLongSegments |= isLongSegment;
302  }
303 
304  r.combineExtentWith( x, y );
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, r, 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  for ( int i = 0; i < numGeoms; ++i )
362  {
363  const QgsAbstractGeometry *sub = srcCollection.geometryN( i );
364  std::unique_ptr< QgsAbstractGeometry > part = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, false );
365  collection->addGeometry( part.release() );
366  }
367  return std::move( collection );
368  }
369  return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
370 }
371 
373 
374 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle &envelope, double map2pixelTol )
375 {
376  // Can replace the geometry by its BBOX ?
377  return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
378 }
379 
381 {
382  if ( geometry.isNull() )
383  {
384  return QgsGeometry();
385  }
387  {
388  return geometry;
389  }
390 
391  // Check whether the geometry can be simplified using the map2pixel context
392  const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry.wkbType() );
393  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
394  if ( flatType == QgsWkbTypes::Point )
395  {
396  return geometry;
397  }
398 
399  const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
400  const int numPoints = geometry.constGet()->nCoordinates();
401 
402  if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
403  {
404  // No simplify simple geometries
405  return geometry;
406  }
407 
408  const QgsRectangle envelope = geometry.boundingBox();
409  if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
410  {
411  //points are in average too far apart to lead to any significant simplification
412  return geometry;
413  }
414 
415  return QgsGeometry( simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry.constGet(), mTolerance, false ) );
416 }
A rectangle specified with double values.
Definition: qgsrectangle.h:40
static Type singleType(Type type)
Returns the single type for a WKB type.
Definition: qgswkbtypes.h:152
void setMinimal()
Set a rectangle so that min corner is at max and max corner is at min.
Definition: qgsrectangle.h:150
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
static bool isMultiType(Type type)
Returns true if the WKB type is a multi type.
Definition: qgswkbtypes.h:557
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:251
const QgsCurve * interiorRing(int i) const
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:104
No simplification can be applied.
virtual QgsAbstractGeometry * createEmptyWithSameType() const =0
Creates a new geometry with the same class and same WKB type as the original and transfers ownership...
The simplification uses a grid (similar to ST_SnapToGrid) to remove duplicate points.
virtual QgsRectangle boundingBox() const =0
Returns the minimal bounding box for the geometry.
SimplifyAlgorithm
Types of simplification algorithms that can be used.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
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 numInteriorRings() const
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:67
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.
SimplifyAlgorithm simplifyAlgorithm() const
Gets the local simplification algorithm of the vector layer managed.
Utility class for identifying a unique vertex within a geometry.
Geometry collection.
The geometries can be fully simplified by its BoundingBox.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:201
QgsGeometryCollection * createEmptyWithSameType() const override
Creates a new geometry with the same class and same WKB type as the original and transfers ownership...
int simplifyFlags() const
Gets the simplification hints of the vector layer managed.
The simplification gives each point in a line an importance weighting, so that least important points...
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
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.
virtual double xAt(int index) const =0
Returns the x-coordinate of the specified node in the line string.
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
Abstract base class for all geometries.
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
double mTolerance
Distance tolerance for the simplification.
int numGeometries() const
Returns the number of geometries within the collection.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:176
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:161
void combineExtentWith(const QgsRectangle &rect)
Expands the rectangle so that it covers both the original rectangle and the given rectangle...
Definition: qgsrectangle.h:352
QgsGeometry simplify(const QgsGeometry &geometry) const override
Returns a simplified version the specified geometry.
QgsMapToPixelSimplifier(int simplifyFlags, double tolerance, SimplifyAlgorithm simplifyAlgorithm=Distance)
Constructor.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
virtual double yAt(int index) const =0
Returns the y-coordinate of the specified node in the line string.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:166
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:171
int mSimplifyFlags
Current simplification flags.
SimplifyAlgorithm mSimplifyAlgorithm
Current algorithm.
The geometries can be simplified using the current map2pixel context state.
Polygon geometry type.
Definition: qgspolygon.h:31
const QgsCurve * exteriorRing() const
virtual int nCoordinates() const
Returns the number of nodes contained in the geometry.
static Type flatType(Type type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:427
virtual int numPoints() const =0
Returns the number of points in the curve.
The simplification uses the distance between points to remove duplicate points.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:208