QGIS API Documentation  2.99.0-Master (90ae728)
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 
41 float QgsMapToPixelSimplifier::calculateLengthSquared2D( double x1, double y1, double x2, double y2 )
42 {
43  float vx = static_cast< float >( x2 - x1 );
44  float vy = static_cast< float >( y2 - y1 );
45 
46  return ( vx * vx ) + ( vy * vy );
47 }
48 
50 bool QgsMapToPixelSimplifier::equalSnapToGrid( double x1, double y1, double x2, double y2, double gridOriginX, double gridOriginY, float gridInverseSizeXY )
51 {
52  int grid_x1 = qRound(( x1 - gridOriginX ) * gridInverseSizeXY );
53  int grid_x2 = qRound(( x2 - gridOriginX ) * gridInverseSizeXY );
54  if ( grid_x1 != grid_x2 ) return false;
55 
56  int grid_y1 = qRound(( y1 - gridOriginY ) * gridInverseSizeXY );
57  int grid_y2 = qRound(( y2 - gridOriginY ) * gridInverseSizeXY );
58  if ( grid_y1 != grid_y2 ) return false;
59 
60  return true;
61 }
62 
64 // Helper simplification methods for Visvalingam method
65 
66 // It uses a refactored code of the liblwgeom implementation:
67 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.h
68 // https://github.com/postgis/postgis/blob/svn-trunk/liblwgeom/effectivearea.c
69 
70 #include "simplify/effectivearea.h"
71 
73 
76  QgsWkbTypes::Type wkbType,
77  const QgsAbstractGeometry& geometry,
78  const QgsRectangle &envelope )
79 {
80  unsigned int geometryType = QgsWkbTypes::singleType( QgsWkbTypes::flatType( wkbType ) );
81 
82  // If the geometry is already minimal skip the generalization
83  int minimumSize = geometryType == QgsWkbTypes::LineString ? 2 : 5;
84 
85  if ( geometry.nCoordinates() <= minimumSize )
86  {
87  return QgsGeometry( geometry.clone() );
88  }
89 
90  const double x1 = envelope.xMinimum();
91  const double y1 = envelope.yMinimum();
92  const double x2 = envelope.xMaximum();
93  const double y2 = envelope.yMaximum();
94 
95  // Write the generalized geometry
96  if ( geometryType == QgsWkbTypes::LineString )
97  {
98  QgsLineString* lineString = new QgsLineString();
99  lineString->addVertex( QgsPointV2( x1, y1 ) );
100  lineString->addVertex( QgsPointV2( x2, y2 ) );
101  return QgsGeometry( lineString );
102  }
103  else
104  {
105  return QgsGeometry::fromRect( envelope );
106  }
107 }
108 
109 template<class T>
110 static T* createEmptySameTypeGeom( const T& geom )
111 {
112  T* output( dynamic_cast<T*>( geom.clone() ) );
113  output->clear();
114  return output;
115 }
116 
118 QgsGeometry QgsMapToPixelSimplifier::simplifyGeometry(
119  int simplifyFlags,
121  QgsWkbTypes::Type wkbType,
122  const QgsAbstractGeometry& geometry,
123  const QgsRectangle &envelope, double map2pixelTol,
124  bool isaLinearRing )
125 {
126  bool isGeneralizable = true;
127 
128  // Can replace the geometry by its BBOX ?
130  isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
131  {
132  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope );
133  }
134 
136  isGeneralizable = false;
137 
138  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( wkbType );
139 
140  // Write the geometry
141  if ( flatType == QgsWkbTypes::LineString || flatType == QgsWkbTypes::CircularString )
142  {
143  const QgsCurve& srcCurve = dynamic_cast<const QgsCurve&>( geometry );
144  std::unique_ptr<QgsCurve> output( createEmptySameTypeGeom( srcCurve ) );
145  double x = 0.0, y = 0.0, lastX = 0.0, lastY = 0.0;
146  QgsRectangle r;
147  r.setMinimal();
148 
149  const int numPoints = srcCurve.numPoints();
150 
151  if ( numPoints <= ( isaLinearRing ? 4 : 2 ) )
152  isGeneralizable = false;
153 
154  bool isLongSegment;
155  bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
156 
157  // Check whether the LinearRing is really closed.
158  if ( isaLinearRing )
159  {
160  isaLinearRing = qgsDoubleNear( srcCurve.xAt( 0 ), srcCurve.xAt( numPoints - 1 ) ) &&
161  qgsDoubleNear( srcCurve.yAt( 0 ), srcCurve.yAt( numPoints - 1 ) );
162  }
163 
164  // Process each vertex...
165  switch ( simplifyAlgorithm )
166  {
167  case SnapToGrid:
168  {
169  double gridOriginX = envelope.xMinimum();
170  double gridOriginY = envelope.yMinimum();
171 
172  // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
173  float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
174 
175  for ( int i = 0; i < numPoints; ++i )
176  {
177  x = srcCurve.xAt( i );
178  y = srcCurve.yAt( i );
179 
180  if ( i == 0 ||
181  !isGeneralizable ||
182  !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
183  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
184  {
185  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPointV2( x, y ) );
186  lastX = x;
187  lastY = y;
188  }
189 
190  r.combineExtentWith( x, y );
191  }
192  break;
193  }
194 
195  case Visvalingam:
196  {
197  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
198 
199  EFFECTIVE_AREAS ea( srcCurve );
200 
201  int set_area = 0;
202  ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
203 
204  for ( int i = 0; i < numPoints; ++i )
205  {
206  if ( ea.res_arealist[ i ] > map2pixelTol )
207  {
208  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
209  }
210  }
211  break;
212  }
213 
214  case Distance:
215  {
216  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
217 
218  for ( int i = 0; i < numPoints; ++i )
219  {
220  x = srcCurve.xAt( i );
221  y = srcCurve.yAt( i );
222 
223  isLongSegment = false;
224 
225  if ( i == 0 ||
226  !isGeneralizable ||
227  ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
228  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
229  {
230  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPointV2( x, y ) );
231  lastX = x;
232  lastY = y;
233 
234  hasLongSegments |= isLongSegment;
235  }
236 
237  r.combineExtentWith( x, y );
238  }
239  }
240  }
241 
242  if ( output->numPoints() < ( isaLinearRing ? 4 : 2 ) )
243  {
244  // we simplified the geometry too much!
245  if ( !hasLongSegments )
246  {
247  // approximate the geometry's shape by its bounding box
248  // (rect for linear ring / one segment for line string)
249  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, r );
250  }
251  else
252  {
253  // Bad luck! The simplified geometry is invalid and approximation by bounding box
254  // would create artifacts due to long segments.
255  // We will return the original geometry
256  return QgsGeometry( geometry.clone() );
257  }
258  }
259 
260  if ( isaLinearRing )
261  {
262  // make sure we keep the linear ring closed
263  if ( !qgsDoubleNear( lastX, output->xAt( 0 ) ) || !qgsDoubleNear( lastY, output->yAt( 0 ) ) )
264  {
265  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPointV2( output->xAt( 0 ), output->yAt( 0 ) ) );
266  }
267  }
268 
269  return QgsGeometry( output.release() );
270  }
271  else if ( flatType == QgsWkbTypes::Polygon )
272  {
273  const QgsPolygonV2& srcPolygon = dynamic_cast<const QgsPolygonV2&>( geometry );
274  std::unique_ptr<QgsPolygonV2> polygon( new QgsPolygonV2() );
275  polygon->setExteriorRing( dynamic_cast<QgsCurve*>( simplifyGeometry( simplifyFlags, simplifyAlgorithm, srcPolygon.exteriorRing()->wkbType(), *srcPolygon.exteriorRing(), envelope, map2pixelTol, true ).geometry()->clone() ) );
276  for ( int i = 0; i < srcPolygon.numInteriorRings(); ++i )
277  {
278  const QgsCurve* sub = srcPolygon.interiorRing( i );
279  polygon->addInteriorRing( dynamic_cast<QgsCurve*>( simplifyGeometry( simplifyFlags, simplifyAlgorithm, sub->wkbType(), *sub, envelope, map2pixelTol, true ).geometry()->clone() ) );
280  }
281  return QgsGeometry( polygon.release() );
282  }
283  else if ( QgsWkbTypes::isMultiType( flatType ) )
284  {
285  const QgsGeometryCollection& srcCollection = dynamic_cast<const QgsGeometryCollection&>( geometry );
286  std::unique_ptr<QgsGeometryCollection> collection( createEmptySameTypeGeom( srcCollection ) );
287  const int numGeoms = srcCollection.numGeometries();
288  for ( int i = 0; i < numGeoms; ++i )
289  {
290  const QgsAbstractGeometry* sub = srcCollection.geometryN( i );
291  collection->addGeometry( simplifyGeometry( simplifyFlags, simplifyAlgorithm, sub->wkbType(), *sub, envelope, map2pixelTol, false ).geometry()->clone() );
292  }
293  return QgsGeometry( collection.release() );
294  }
295  return QgsGeometry( geometry.clone() );
296 }
297 
299 
301 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle& envelope, double map2pixelTol )
302 {
303  // Can replace the geometry by its BBOX ?
304  return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
305 }
306 
309 {
310  if ( geometry.isNull() )
311  {
312  return QgsGeometry();
313  }
315  {
316  return geometry;
317  }
318 
319  // Check whether the geometry can be simplified using the map2pixel context
320  const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry.wkbType() );
321  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
322  if ( flatType == QgsWkbTypes::Point )
323  {
324  return geometry;
325  }
326 
327  const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
328  const int numPoints = geometry.geometry()->nCoordinates();
329 
330  if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
331  {
332  // No simplify simple geometries
333  return geometry;
334  }
335 
336  const QgsRectangle envelope = geometry.boundingBox();
337  if ( qMax( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
338  {
339  //points are in average too far apart to lead to any significant simplification
340  return geometry;
341  }
342 
343  return simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, geometry.wkbType(), *geometry.geometry(), envelope, mTolerance, false );
344 }
const QgsCurve * interiorRing(int i) const
A rectangle specified with double values.
Definition: qgsrectangle.h:36
static Type singleType(Type type)
Returns the single type for a WKB type.
Definition: qgswkbtypes.h:145
void setMinimal()
Set a rectangle so that min corner is at max and max corner is at min.
void addVertex(const QgsPointV2 &pt)
Adds a new vertex to the end of the line string.
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:518
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
static QgsGeometry generalizeWkbGeometryByBoundingBox(QgsWkbTypes::Type wkbType, const QgsAbstractGeometry &geometry, const QgsRectangle &envelope)
Generalize the WKB-geometry using the BBOX of the original geometry.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:79
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:198
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:65
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
Width of the rectangle.
Definition: qgsrectangle.h:212
Point geometry type, with support for z-dimension and m-values.
Definition: qgspointv2.h:36
int simplifyFlags() const
Gets the simplification hints of the vector layer managed.
static T * createEmptySameTypeGeom(const T &geom)
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:33
Abstract base class for all geometries.
const QgsCurve * exteriorRing() const
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
double mTolerance
Distance tolerance for the simplification.
int numGeometries() const
Returns the number of geometries within the collection.
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:207
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:192
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:36
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
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:197
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:202
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:397
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
Height of the rectangle.
Definition: qgsrectangle.h:217