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