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