QGIS API Documentation  2.99.0-Master (9fdd060)
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"
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 QgsGeometry generalizeWkbGeometryByBoundingBox(
72  QgsWkbTypes::Type wkbType,
73  const QgsAbstractGeometry &geometry,
74  const QgsRectangle &envelope )
75 {
76  unsigned int geometryType = QgsWkbTypes::singleType( QgsWkbTypes::flatType( wkbType ) );
77
78  // If the geometry is already minimal skip the generalization
79  int minimumSize = geometryType == QgsWkbTypes::LineString ? 2 : 5;
80
81  if ( geometry.nCoordinates() <= minimumSize )
82  {
83  return QgsGeometry( geometry.clone() );
84  }
85
86  const double x1 = envelope.xMinimum();
87  const double y1 = envelope.yMinimum();
88  const double x2 = envelope.xMaximum();
89  const double y2 = envelope.yMaximum();
90
91  // Write the generalized geometry
92  if ( geometryType == QgsWkbTypes::LineString )
93  {
94  QgsLineString *lineString = new QgsLineString();
95  lineString->addVertex( QgsPoint( x1, y1 ) );
96  lineString->addVertex( QgsPoint( x2, y2 ) );
97  return QgsGeometry( lineString );
98  }
99  else
100  {
101  return QgsGeometry::fromRect( envelope );
102  }
103 }
104
105 template<class T>
106 static T *createEmptySameTypeGeom( const T &geom )
107 {
108  // TODO - this is inefficient - we clone the geometry's content only to throw it away immediately
109  T *output( qgsgeometry_cast<T *>( geom.clone() ) );
110  Q_ASSERT( output );
111  output->clear();
112  return output;
113 }
114
115 QgsGeometry QgsMapToPixelSimplifier::simplifyGeometry(
116  int simplifyFlags,
118  QgsWkbTypes::Type wkbType,
119  const QgsAbstractGeometry &geometry,
120  const QgsRectangle &envelope, double map2pixelTol,
121  bool isaLinearRing )
122 {
123  bool isGeneralizable = true;
124
125  // Can replace the geometry by its BBOX ?
127  isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
128  {
129  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope );
130  }
131
133  isGeneralizable = false;
134
135  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( wkbType );
136
137  // Write the geometry
138  if ( flatType == QgsWkbTypes::LineString || flatType == QgsWkbTypes::CircularString )
139  {
140  const QgsCurve &srcCurve = dynamic_cast<const QgsCurve &>( geometry );
141  std::unique_ptr<QgsCurve> output( createEmptySameTypeGeom( srcCurve ) );
142  double x = 0.0, y = 0.0, lastX = 0.0, lastY = 0.0;
143  QgsRectangle r;
144  r.setMinimal();
145
146  const int numPoints = srcCurve.numPoints();
147
148  if ( numPoints <= ( isaLinearRing ? 4 : 2 ) )
149  isGeneralizable = false;
150
151  bool isLongSegment;
152  bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
153
154  // Check whether the LinearRing is really closed.
155  if ( isaLinearRing )
156  {
157  isaLinearRing = qgsDoubleNear( srcCurve.xAt( 0 ), srcCurve.xAt( numPoints - 1 ) ) &&
158  qgsDoubleNear( srcCurve.yAt( 0 ), srcCurve.yAt( numPoints - 1 ) );
159  }
160
161  // Process each vertex...
162  switch ( simplifyAlgorithm )
163  {
164  case SnapToGrid:
165  {
166  double gridOriginX = envelope.xMinimum();
167  double gridOriginY = envelope.yMinimum();
168
169  // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
170  float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
171
172  for ( int i = 0; i < numPoints; ++i )
173  {
174  x = srcCurve.xAt( i );
175  y = srcCurve.yAt( i );
176
177  if ( i == 0 ||
178  !isGeneralizable ||
179  !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
180  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
181  {
182  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
183  lastX = x;
184  lastY = y;
185  }
186
187  r.combineExtentWith( x, y );
188  }
189  break;
190  }
191
192  case Visvalingam:
193  {
194  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
195
196  EFFECTIVE_AREAS ea( srcCurve );
197
198  int set_area = 0;
199  ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
200
201  for ( int i = 0; i < numPoints; ++i )
202  {
203  if ( ea.res_arealist[ i ] > map2pixelTol )
204  {
205  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
206  }
207  }
208  break;
209  }
210
211  case Distance:
212  {
213  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
214
215  for ( int i = 0; i < numPoints; ++i )
216  {
217  x = srcCurve.xAt( i );
218  y = srcCurve.yAt( i );
219
220  isLongSegment = false;
221
222  if ( i == 0 ||
223  !isGeneralizable ||
224  ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
225  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
226  {
227  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
228  lastX = x;
229  lastY = y;
230
231  hasLongSegments |= isLongSegment;
232  }
233
234  r.combineExtentWith( x, y );
235  }
236  }
237  }
238
239  if ( output->numPoints() < ( isaLinearRing ? 4 : 2 ) )
240  {
241  // we simplified the geometry too much!
242  if ( !hasLongSegments )
243  {
244  // approximate the geometry's shape by its bounding box
245  // (rect for linear ring / one segment for line string)
246  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, r );
247  }
248  else
249  {
250  // Bad luck! The simplified geometry is invalid and approximation by bounding box
251  // would create artifacts due to long segments.
252  // We will return the original geometry
253  return QgsGeometry( geometry.clone() );
254  }
255  }
256
257  if ( isaLinearRing )
258  {
259  // make sure we keep the linear ring closed
260  if ( !qgsDoubleNear( lastX, output->xAt( 0 ) ) || !qgsDoubleNear( lastY, output->yAt( 0 ) ) )
261  {
262  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( output->xAt( 0 ), output->yAt( 0 ) ) );
263  }
264  }
265
266  return QgsGeometry( output.release() );
267  }
268  else if ( flatType == QgsWkbTypes::Polygon )
269  {
270  const QgsPolygonV2 &srcPolygon = dynamic_cast<const QgsPolygonV2 &>( geometry );
271  std::unique_ptr<QgsPolygonV2> polygon( new QgsPolygonV2() );
272  polygon->setExteriorRing( qgsgeometry_cast<QgsCurve *>( simplifyGeometry( simplifyFlags, simplifyAlgorithm, srcPolygon.exteriorRing()->wkbType(), *srcPolygon.exteriorRing(), envelope, map2pixelTol, true ).geometry()->clone() ) );
273  for ( int i = 0; i < srcPolygon.numInteriorRings(); ++i )
274  {
275  const QgsCurve *sub = srcPolygon.interiorRing( i );
276  polygon->addInteriorRing( qgsgeometry_cast<QgsCurve *>( simplifyGeometry( simplifyFlags, simplifyAlgorithm, sub->wkbType(), *sub, envelope, map2pixelTol, true ).geometry()->clone() ) );
277  }
278  return QgsGeometry( polygon.release() );
279  }
280  else if ( QgsWkbTypes::isMultiType( flatType ) )
281  {
282  const QgsGeometryCollection &srcCollection = dynamic_cast<const QgsGeometryCollection &>( geometry );
283  std::unique_ptr<QgsGeometryCollection> collection( createEmptySameTypeGeom( srcCollection ) );
284  const int numGeoms = srcCollection.numGeometries();
285  for ( int i = 0; i < numGeoms; ++i )
286  {
287  const QgsAbstractGeometry *sub = srcCollection.geometryN( i );
288  collection->addGeometry( simplifyGeometry( simplifyFlags, simplifyAlgorithm, sub->wkbType(), *sub, envelope, map2pixelTol, false ).geometry()->clone() );
289  }
290  return QgsGeometry( collection.release() );
291  }
292  return QgsGeometry( geometry.clone() );
293 }
294
296
297 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle &envelope, double map2pixelTol )
298 {
299  // Can replace the geometry by its BBOX ?
300  return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
301 }
302
304 {
305  if ( geometry.isNull() )
306  {
307  return QgsGeometry();
308  }
310  {
311  return geometry;
312  }
313
314  // Check whether the geometry can be simplified using the map2pixel context
315  const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry.wkbType() );
316  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
317  if ( flatType == QgsWkbTypes::Point )
318  {
319  return geometry;
320  }
321
322  const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
323  const int numPoints = geometry.geometry()->nCoordinates();
324
325  if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
326  {
327  // No simplify simple geometries
328  return geometry;
329  }
330
331  const QgsRectangle envelope = geometry.boundingBox();
332  if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
333  {
334  //points are in average too far apart to lead to any significant simplification
335  return geometry;
336  }
337
338  return simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, geometry.wkbType(), *geometry.geometry(), envelope, mTolerance, false );
339 }
const QgsCurve * interiorRing(int i) const
A rectangle specified with double values.
Definition: qgsrectangle.h:39
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.
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.)
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:94
No simplification can be applied.
The simplification uses a grid (similar to ST_SnapToGrid) to remove duplicate points.
bool insertVertex(double x, double y, int beforeVertex)
Insert a new vertex before the given vertex index, ring and item (first number is index 0) If the req...
SimplifyAlgorithm
Types of simplification algorithms that can be used.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:227
Polygon geometry type.
Definition: qgspolygon.h:31
static QgsGeometry fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
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.
Adds a new vertex to the end of the line string.
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:131
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...
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.
const QgsCurve * exteriorRing() const
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.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:119
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:104
void combineExtentWith(const QgsRectangle &rect)
Expand the rectangle so that covers both the original rectangle and the given rectangle.
virtual QgsGeometry simplify(const QgsGeometry &geometry) const override
Returns a simplified version the specified geometry.
QgsMapToPixelSimplifier(int simplifyFlags, double tolerance, SimplifyAlgorithm simplifyAlgorithm=Distance)
Constructor.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:41
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
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:109
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:114
int mSimplifyFlags
Current simplification flags.
SimplifyAlgorithm mSimplifyAlgorithm
Current algorithm.
The geometries can be simplified using the current map2pixel context state.
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.
QgsAbstractGeometry * geometry() const
Returns the underlying geometry store.
The simplification uses the distance between points to remove duplicate points.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:138