QGIS API Documentation  3.0.2-Girona (307d082)
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 = 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  return QgsGeometry( qgis::make_unique< QgsLineString >( QVector<double>() << x1 << x2, QVector<double>() << y1 << y2 ) );
95  }
96  else
97  {
98  return QgsGeometry::fromRect( envelope );
99  }
100 }
101 
102 QgsGeometry QgsMapToPixelSimplifier::simplifyGeometry(
103  int simplifyFlags,
105  QgsWkbTypes::Type wkbType,
106  const QgsAbstractGeometry &geometry,
107  const QgsRectangle &envelope, double map2pixelTol,
108  bool isaLinearRing )
109 {
110  bool isGeneralizable = true;
111 
112  // Can replace the geometry by its BBOX ?
114  isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
115  {
116  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope );
117  }
118 
120  isGeneralizable = false;
121 
122  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( wkbType );
123 
124  // Write the geometry
125  if ( flatType == QgsWkbTypes::LineString || flatType == QgsWkbTypes::CircularString )
126  {
127  const QgsCurve &srcCurve = dynamic_cast<const QgsCurve &>( geometry );
128  const int numPoints = srcCurve.numPoints();
129 
130  std::unique_ptr<QgsCurve> output;
131 
132  QVector< double > lineStringX;
133  QVector< double > lineStringY;
134  if ( flatType == QgsWkbTypes::LineString )
135  {
136  // if we are making a linestring, we do it in an optimised way by directly constructing
137  // the final x/y vectors, which avoids calling the slower insertVertex method
138  lineStringX.reserve( numPoints );
139  lineStringY.reserve( numPoints );
140  }
141  else
142  {
143  output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.createEmptyWithSameType() ) );
144  }
145 
146  double x = 0.0, y = 0.0, lastX = 0.0, lastY = 0.0;
147  QgsRectangle r;
148  r.setMinimal();
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  if ( output )
185  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
186  else
187  {
188  lineStringX.append( x );
189  lineStringY.append( y );
190  }
191  lastX = x;
192  lastY = y;
193  }
194 
195  r.combineExtentWith( x, y );
196  }
197  break;
198  }
199 
200  case Visvalingam:
201  {
202  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
203 
204  EFFECTIVE_AREAS ea( srcCurve );
205 
206  int set_area = 0;
207  ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
208 
209  for ( int i = 0; i < numPoints; ++i )
210  {
211  if ( ea.res_arealist[ i ] > map2pixelTol )
212  {
213  if ( output )
214  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
215  else
216  {
217  lineStringX.append( ea.inpts.at( i ).x() );
218  lineStringY.append( ea.inpts.at( i ).y() );
219  }
220  }
221  }
222  break;
223  }
224 
225  case Distance:
226  {
227  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
228 
229  for ( int i = 0; i < numPoints; ++i )
230  {
231  x = srcCurve.xAt( i );
232  y = srcCurve.yAt( i );
233 
234  isLongSegment = false;
235 
236  if ( i == 0 ||
237  !isGeneralizable ||
238  ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
239  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
240  {
241  if ( output )
242  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
243  else
244  {
245  lineStringX.append( x );
246  lineStringY.append( y );
247  }
248  lastX = x;
249  lastY = y;
250 
251  hasLongSegments |= isLongSegment;
252  }
253 
254  r.combineExtentWith( x, y );
255  }
256  }
257  }
258 
259  if ( !output )
260  {
261  output = qgis::make_unique< QgsLineString >( lineStringX, lineStringY );
262  }
263  if ( output->numPoints() < ( isaLinearRing ? 4 : 2 ) )
264  {
265  // we simplified the geometry too much!
266  if ( !hasLongSegments )
267  {
268  // approximate the geometry's shape by its bounding box
269  // (rect for linear ring / one segment for line string)
270  return generalizeWkbGeometryByBoundingBox( wkbType, geometry, r );
271  }
272  else
273  {
274  // Bad luck! The simplified geometry is invalid and approximation by bounding box
275  // would create artifacts due to long segments.
276  // We will return the original geometry
277  return QgsGeometry( geometry.clone() );
278  }
279  }
280 
281  if ( isaLinearRing )
282  {
283  // make sure we keep the linear ring closed
284  if ( !qgsDoubleNear( lastX, output->xAt( 0 ) ) || !qgsDoubleNear( lastY, output->yAt( 0 ) ) )
285  {
286  output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( output->xAt( 0 ), output->yAt( 0 ) ) );
287  }
288  }
289 
290  return QgsGeometry( output.release() );
291  }
292  else if ( flatType == QgsWkbTypes::Polygon )
293  {
294  const QgsPolygon &srcPolygon = dynamic_cast<const QgsPolygon &>( geometry );
295  std::unique_ptr<QgsPolygon> polygon( new QgsPolygon() );
296  polygon->setExteriorRing( qgsgeometry_cast<QgsCurve *>( simplifyGeometry( simplifyFlags, simplifyAlgorithm, srcPolygon.exteriorRing()->wkbType(), *srcPolygon.exteriorRing(), envelope, map2pixelTol, true ).constGet()->clone() ) );
297  for ( int i = 0; i < srcPolygon.numInteriorRings(); ++i )
298  {
299  const QgsCurve *sub = srcPolygon.interiorRing( i );
300  polygon->addInteriorRing( qgsgeometry_cast<QgsCurve *>( simplifyGeometry( simplifyFlags, simplifyAlgorithm, sub->wkbType(), *sub, envelope, map2pixelTol, true ).constGet()->clone() ) );
301  }
302  return QgsGeometry( polygon.release() );
303  }
304  else if ( QgsWkbTypes::isMultiType( flatType ) )
305  {
306  const QgsGeometryCollection &srcCollection = dynamic_cast<const QgsGeometryCollection &>( geometry );
307  std::unique_ptr<QgsGeometryCollection> collection( srcCollection.createEmptyWithSameType() );
308  const int numGeoms = srcCollection.numGeometries();
309  for ( int i = 0; i < numGeoms; ++i )
310  {
311  const QgsAbstractGeometry *sub = srcCollection.geometryN( i );
312  collection->addGeometry( simplifyGeometry( simplifyFlags, simplifyAlgorithm, sub->wkbType(), *sub, envelope, map2pixelTol, false ).constGet()->clone() );
313  }
314  return QgsGeometry( collection.release() );
315  }
316  return QgsGeometry( geometry.clone() );
317 }
318 
320 
321 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle &envelope, double map2pixelTol )
322 {
323  // Can replace the geometry by its BBOX ?
324  return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
325 }
326 
328 {
329  if ( geometry.isNull() )
330  {
331  return QgsGeometry();
332  }
334  {
335  return geometry;
336  }
337 
338  // Check whether the geometry can be simplified using the map2pixel context
339  const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry.wkbType() );
340  const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
341  if ( flatType == QgsWkbTypes::Point )
342  {
343  return geometry;
344  }
345 
346  const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
347  const int numPoints = geometry.constGet()->nCoordinates();
348 
349  if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
350  {
351  // No simplify simple geometries
352  return geometry;
353  }
354 
355  const QgsRectangle envelope = geometry.boundingBox();
356  if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
357  {
358  //points are in average too far apart to lead to any significant simplification
359  return geometry;
360  }
361 
362  return simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, geometry.wkbType(), *geometry.constGet(), envelope, mTolerance, false );
363 }
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:111
No simplification can be applied.
virtual QgsAbstractGeometry * createEmptyWithSameType() const =0
Creates a new geometry with the same class and same WKB type as the original and transfers ownership...
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:251
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.
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:142
QgsGeometryCollection * createEmptyWithSameType() const override
Creates a new geometry with the same class and same WKB type as the original and transfers ownership...
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.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:130
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:115
void combineExtentWith(const QgsRectangle &rect)
Expand the rectangle so that covers both the original rectangle and the given rectangle.
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.
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:120
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:125
int mSimplifyFlags
Current simplification flags.
SimplifyAlgorithm mSimplifyAlgorithm
Current algorithm.
The geometries can be simplified using the current map2pixel context state.
Polygon geometry type.
Definition: qgspolygon.h:31
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.
The simplification uses the distance between points to remove duplicate points.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:149