QGIS API Documentation  2.99.0-Master (23ddace)
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 = qRound( ( x1 - gridOriginX ) * gridInverseSizeXY );
51  int grid_x2 = qRound( ( x2 - gridOriginX ) * gridInverseSizeXY );
52  if ( grid_x1 != grid_x2 ) return false;
53 
54  int grid_y1 = qRound( ( y1 - gridOriginY ) * gridInverseSizeXY );
55  int grid_y2 = qRound( ( y2 - gridOriginY ) * gridInverseSizeXY );
56  if ( grid_y1 != grid_y2 ) return false;
57 
58  return true;
59 }
60 
62 // Helper simplification methods for Visvalingam method
63 
64 // It uses a refactored code of the liblwgeom implementation:
65 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.h
66 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.c
67 
68 #include "simplify/effectivearea.h"
69 
71 
73 static QgsGeometry generalizeWkbGeometryByBoundingBox(
74  QgsWkbTypes::Type wkbType,
75  const QgsAbstractGeometry &geometry,
76  const QgsRectangle &envelope )
77 {
78  unsigned int geometryType = QgsWkbTypes::singleType( QgsWkbTypes::flatType( wkbType ) );
79 
80  // If the geometry is already minimal skip the generalization
81  int minimumSize = geometryType == QgsWkbTypes::LineString ? 2 : 5;
82 
83  if ( geometry.nCoordinates() <= minimumSize )
84  {
85  return QgsGeometry( 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 )
95  {
96  QgsLineString *lineString = new QgsLineString();
97  lineString->addVertex( QgsPoint( x1, y1 ) );
98  lineString->addVertex( QgsPoint( x2, y2 ) );
99  return QgsGeometry( lineString );
100  }
101  else
102  {
103  return QgsGeometry::fromRect( envelope );
104  }
105 }
106 
107 template<class T>
108 static T *createEmptySameTypeGeom( const T &geom )
109 {
110  T *output( dynamic_cast<T *>( geom.clone() ) );
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( dynamic_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( dynamic_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 ( qMax( 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:38
static Type singleType(Type type)
Returns the single type for a WKB type.
Definition: qgswkbtypes.h:150
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:550
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:96
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:203
Polygon geometry type.
Definition: qgspolygon.h:30
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:66
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.
void addVertex(const QgsPoint &pt)
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:118
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:34
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:36
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:106
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:91
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:40
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:96
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:101
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:423
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:125