QGIS API Documentation
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 
23 
24 QgsMapToPixelSimplifier::QgsMapToPixelSimplifier( int simplifyFlags, double tolerance )
25  : mSimplifyFlags( simplifyFlags )
26  , mTolerance( tolerance )
27 {
28 }
29 
31 {
32 }
33 
35 // Helper simplification methods
36 
38 float QgsMapToPixelSimplifier::calculateLengthSquared2D( double x1, double y1, double x2, double y2 )
39 {
40  float vx = static_cast< float >( x2 - x1 );
41  float vy = static_cast< float >( y2 - y1 );
42 
43  return vx*vx + vy*vy;
44 }
45 
47 inline static QgsRectangle calculateBoundingBox( QGis::WkbType wkbType, QgsConstWkbPtr wkbPtr, int numPoints )
48 {
49  QgsRectangle r;
50  r.setMinimal();
51 
52  int skipZM = ( QGis::wkbDimensions( wkbType ) - 2 ) * sizeof( double );
53  Q_ASSERT( skipZM >= 0 );
54 
55  for ( int i = 0; i < numPoints; ++i )
56  {
57  double x, y;
58  wkbPtr >> x >> y;
59  wkbPtr += skipZM;
60  r.combineExtentWith( x, y );
61  }
62 
63  return r;
64 }
65 
68  QGis::WkbType wkbType,
69  QgsConstWkbPtr sourceWkbPtr,
70  QgsWkbPtr targetWkbPtr,
71  int &targetWkbSize,
72  const QgsRectangle &envelope, bool writeHeader )
73 {
74  QgsWkbPtr savedTargetWkb( targetWkbPtr );
75  unsigned int geometryType = QGis::singleType( QGis::flatType( wkbType ) );
76 
77  int skipZM = ( QGis::wkbDimensions( wkbType ) - 2 ) * sizeof( double );
78  Q_ASSERT( skipZM >= 0 );
79 
80  // If the geometry is already minimal skip the generalization
81  int minimumSize = geometryType == QGis::WKBLineString ? 4 + 2 * ( 2 * sizeof( double ) + skipZM ) : 8 + 5 * ( 2 * sizeof( double ) + skipZM );
82 
83  if ( writeHeader )
84  minimumSize += 5;
85 
86  if ( sourceWkbPtr.remaining() <= minimumSize )
87  {
88  targetWkbSize = 0;
89  return false;
90  }
91 
92  double x1 = envelope.xMinimum();
93  double y1 = envelope.yMinimum();
94  double x2 = envelope.xMaximum();
95  double y2 = envelope.yMaximum();
96 
97  // Write the main header of the geometry
98  if ( writeHeader )
99  {
100  targetWkbPtr << ( char ) QgsApplication::endian() << geometryType;
101 
102  if ( geometryType == QGis::WKBPolygon ) // numRings
103  {
104  targetWkbPtr << 1;
105  }
106  }
107 
108  // Write the generalized geometry
109  if ( geometryType == QGis::WKBLineString )
110  {
111  targetWkbPtr << 2 << x1 << y1 << x2 << y2;
112  }
113  else
114  {
115  targetWkbPtr << 5 << x1 << y1 << x2 << y1 << x2 << y2 << x1 << y2 << x1 << y1;
116  }
117 
118  targetWkbSize += targetWkbPtr - savedTargetWkb;
119 
120  return true;
121 }
122 
124 bool QgsMapToPixelSimplifier::simplifyWkbGeometry(
125  int simplifyFlags, QGis::WkbType wkbType,
126  QgsConstWkbPtr sourceWkbPtr,
127  QgsWkbPtr targetWkbPtr,
128  int &targetWkbSize,
129  const QgsRectangle &envelope, double map2pixelTol,
130  bool writeHeader, bool isaLinearRing )
131 {
132  bool isGeneralizable = true;
133  bool result = false;
134 
135  // Save initial WKB settings to use when the simplification creates invalid geometries
136  QgsConstWkbPtr sourcePrevWkbPtr( sourceWkbPtr );
137  QgsWkbPtr targetPrevWkbPtr( targetWkbPtr );
138  int targetWkbPrevSize = targetWkbSize;
139 
140  // Can replace the geometry by its BBOX ?
142  isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
143  {
144  isGeneralizable = generalizeWkbGeometryByBoundingBox( wkbType, sourceWkbPtr, targetWkbPtr, targetWkbSize, envelope, writeHeader );
145  if ( isGeneralizable )
146  return true;
147  }
148 
150  isGeneralizable = false;
151 
152  // Write the main header of the geometry
153  if ( writeHeader )
154  {
155  QgsWKBTypes::Type geometryType = sourceWkbPtr.readHeader();
156 
157  targetWkbPtr << ( char ) QgsApplication::endian() << QgsWKBTypes::flatType( geometryType );
158 
159  targetWkbSize += targetWkbPtr - targetPrevWkbPtr;
160  }
161 
162  unsigned int flatType = QGis::flatType( wkbType );
163 
164  // Write the geometry
165  if ( flatType == QGis::WKBLineString || isaLinearRing )
166  {
167  QgsWkbPtr savedTargetWkbPtr( targetWkbPtr );
168  double x = 0.0, y = 0.0, lastX = 0, lastY = 0;
169  QgsRectangle r;
170  r.setMinimal();
171 
172  int skipZM = ( QGis::wkbDimensions( wkbType ) - 2 ) * sizeof( double );
173  Q_ASSERT( skipZM >= 0 );
174 
175  int numPoints;
176  sourceWkbPtr >> numPoints;
177 
178  if ( numPoints <= ( isaLinearRing ? 5 : 2 ) )
179  isGeneralizable = false;
180 
181  QgsWkbPtr numPtr( targetWkbPtr );
182 
183  int numTargetPoints = 0;
184  targetWkbPtr << numTargetPoints;
185  targetWkbSize += 4;
186 
187  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
188 
189  bool isLongSegment;
190  bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
191 
192  // Check whether the LinearRing is really closed.
193  if ( isaLinearRing )
194  {
195  QgsConstWkbPtr checkPtr( sourceWkbPtr );
196 
197  double x1, y1, x2, y2;
198 
199  checkPtr >> x1 >> y1;
200  checkPtr += skipZM + ( numPoints - 2 ) * ( 2 * sizeof( double ) + skipZM );
201  checkPtr >> x2 >> y2;
202 
203  isaLinearRing = qgsDoubleNear( x1, x2 ) && qgsDoubleNear( y1, y2 );
204  }
205 
206  // Process each vertex...
207  for ( int i = 0; i < numPoints; ++i )
208  {
209  sourceWkbPtr >> x >> y;
210  sourceWkbPtr += skipZM;
211 
212  isLongSegment = false;
213 
214  if ( i == 0 ||
215  !isGeneralizable ||
216  ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
217  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
218  {
219  targetWkbPtr << x << y;
220  lastX = x;
221  lastY = y;
222  numTargetPoints++;
223 
224  hasLongSegments |= isLongSegment;
225  }
226 
227  r.combineExtentWith( x, y );
228  }
229 
230  QgsWkbPtr nextPointPtr( targetWkbPtr );
231 
232  targetWkbPtr = savedTargetWkbPtr;
233  targetWkbPtr += sizeof( int );
234 
235  if ( numTargetPoints < ( isaLinearRing ? 4 : 2 ) )
236  {
237  // we simplified the geometry too much!
238  if ( !hasLongSegments )
239  {
240  // approximate the geometry's shape by its bounding box
241  // (rect for linear ring / one segment for line string)
242  QgsWkbPtr tempWkbPtr( targetWkbPtr );
243  int targetWkbTempSize = targetWkbSize;
244 
245  sourceWkbPtr = sourcePrevWkbPtr;
246  targetWkbPtr = targetPrevWkbPtr;
247  targetWkbSize = targetWkbPrevSize;
248  if ( generalizeWkbGeometryByBoundingBox( wkbType, sourceWkbPtr, targetWkbPtr, targetWkbSize, r, writeHeader ) )
249  return true;
250 
251  targetWkbPtr = tempWkbPtr;
252  targetWkbSize = targetWkbTempSize;
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. Worst of all, we may have overwritten
258  // the original coordinates by the simplified ones (source and target WKB ptr can be the same)
259  // so we cannot even undo the changes here. We will return invalid geometry and hope that
260  // other pieces of QGIS will survive that :-/
261  }
262  }
263 
264  if ( isaLinearRing )
265  {
266  // make sure we keep the linear ring closed
267  targetWkbPtr << x << y;
268  if ( !qgsDoubleNear( lastX, x ) || !qgsDoubleNear( lastY, y ) )
269  {
270  nextPointPtr << x << y;
271  numTargetPoints++;
272  }
273  }
274 
275  numPtr << numTargetPoints;
276  targetWkbSize += numTargetPoints * sizeof( double ) * 2;
277 
278  result = numPoints != numTargetPoints;
279  }
280  else if ( flatType == QGis::WKBPolygon )
281  {
282  int numRings;
283  sourceWkbPtr >> numRings;
284  targetWkbPtr << numRings;
285  targetWkbSize += 4;
286 
287  for ( int i = 0; i < numRings; ++i )
288  {
289  int numPoints_i;
290  sourceWkbPtr >> numPoints_i;
291 
292  QgsRectangle envelope_i = numRings == 1 ? envelope : calculateBoundingBox( wkbType, sourceWkbPtr, numPoints_i );
293 
294  sourceWkbPtr -= sizeof( int );
295 
296  int sourceWkbSize_i = sizeof( int ) + numPoints_i * QGis::wkbDimensions( wkbType ) * sizeof( double );
297  int targetWkbSize_i = 0;
298 
299  result |= simplifyWkbGeometry( simplifyFlags, wkbType, sourceWkbPtr, targetWkbPtr, targetWkbSize_i, envelope_i, map2pixelTol, false, true );
300  sourceWkbPtr += sourceWkbSize_i;
301  targetWkbPtr += targetWkbSize_i;
302 
303  targetWkbSize += targetWkbSize_i;
304  }
305  }
306  else if ( flatType == QGis::WKBMultiLineString || flatType == QGis::WKBMultiPolygon )
307  {
308  int numGeoms;
309  sourceWkbPtr >> numGeoms;
310  targetWkbPtr << numGeoms;
311  targetWkbSize += 4;
312 
313  QgsConstWkbPtr sourceWkbPtr2( sourceWkbPtr );
314 
315  for ( int i = 0; i < numGeoms; ++i )
316  {
317  int sourceWkbSize_i = 0;
318  int targetWkbSize_i = 0;
319 
320  sourceWkbPtr2.readHeader();
321 
322  // ... calculate the wkb-size of the current child complex geometry
323  if ( flatType == QGis::WKBMultiLineString )
324  {
325  int numPoints_i;
326  sourceWkbPtr2 >> numPoints_i;
327 
328  int wkbSize_i = numPoints_i * QGis::wkbDimensions( wkbType ) * sizeof( double );
329 
330  sourceWkbSize_i += 9 + wkbSize_i;
331  sourceWkbPtr2 += wkbSize_i;
332  }
333  else
334  {
335  int numPrings_i;
336  sourceWkbPtr2 >> numPrings_i;
337  sourceWkbSize_i = 1 + 2 * sizeof( int );
338 
339  for ( int j = 0; j < numPrings_i; ++j )
340  {
341  int numPoints_i;
342  sourceWkbPtr2 >> numPoints_i;
343 
344  int wkbSize_i = numPoints_i * QGis::wkbDimensions( wkbType ) * sizeof( double );
345 
346  sourceWkbSize_i += 4 + wkbSize_i;
347  sourceWkbPtr2 += wkbSize_i;
348  }
349  }
350  result |= simplifyWkbGeometry( simplifyFlags, QGis::singleType( wkbType ), sourceWkbPtr, targetWkbPtr, targetWkbSize_i, envelope, map2pixelTol, true, false );
351  sourceWkbPtr += sourceWkbSize_i;
352  targetWkbPtr += targetWkbSize_i;
353 
354  targetWkbSize += targetWkbSize_i;
355  }
356  }
357 
358  return result;
359 }
360 
362 
364 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle& envelope, double map2pixelTol )
365 {
366  // Can replace the geometry by its BBOX ?
367  return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
368 }
369 
372 {
373  QgsGeometry* g = new QgsGeometry();
374 
375  int wkbSize = geometry->wkbSize();
376  unsigned char* wkb = new unsigned char[ wkbSize ];
377  memcpy( wkb, geometry->asWkb(), wkbSize );
378  g->fromWkb( wkb, wkbSize );
380 
381  return g;
382 }
383 
386 {
387  int finalWkbSize = 0;
388 
389  // Check whether the geometry can be simplified using the map2pixel context
390  QGis::GeometryType geometryType = geometry->type();
391  if ( !( geometryType == QGis::Line || geometryType == QGis::Polygon ) )
392  return false;
393 
394  QgsRectangle envelope = geometry->boundingBox();
395  QGis::WkbType wkbType = geometry->wkbType();
396 
397  QgsConstWkbPtr wkbPtr( geometry->asWkb(), geometry->wkbSize() );
398 
399  unsigned char* targetWkb = new unsigned char[wkbPtr.remaining()];
400  memcpy( targetWkb, wkbPtr, wkbPtr.remaining() );
401  QgsWkbPtr targetWkbPtr( targetWkb, wkbPtr.remaining() );
402 
403  try
404  {
405  if ( simplifyWkbGeometry( simplifyFlags, wkbType, wkbPtr, targetWkbPtr, finalWkbSize, envelope, tolerance ) )
406  {
407  unsigned char *finalWkb = new unsigned char[finalWkbSize];
408  memcpy( finalWkb, targetWkb, finalWkbSize );
409  geometry->fromWkb( finalWkb, finalWkbSize );
410  delete [] targetWkb;
411  return true;
412  }
413  }
414  catch ( const QgsWkbException &e )
415  {
416  Q_UNUSED( e );
417  QgsDebugMsg( QString( "Exception thrown by simplifier: %1" ) .arg( e.what() ) );
418  }
419  delete [] targetWkb;
420  return false;
421 }
422 
425 {
426  return simplifyGeometry( geometry, mSimplifyFlags, mTolerance );
427 }
static WkbType flatType(WkbType type)
Map 2d+ to 2d type.
Definition: qgis.cpp:395
A rectangle specified with double values.
Definition: qgsrectangle.h:35
void setMinimal()
Set a rectangle so that min corner is at max and max corner is at min.
GeometryType
Definition: qgis.h:111
static int wkbDimensions(WkbType type)
Definition: qgis.cpp:426
int remaining() const
Definition: qgswkbptr.h:121
void fromWkb(unsigned char *wkb, int length)
Set the geometry, feeding in the buffer containing OGC Well-Known Binary and the buffer&#39;s length...
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:197
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
virtual bool simplifyGeometry(QgsGeometry *geometry) const override
Simplifies the specified geometry.
QgsRectangle boundingBox() const
Returns the bounding box of this feature.
QGis::GeometryType type() const
Returns type of the geometry as a QGis::GeometryType.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:76
virtual QgsGeometry * simplify(QgsGeometry *geometry) const override
Returns a simplified version the specified geometry.
WkbType
Used for symbology operations.
Definition: qgis.h:57
static endian_t endian()
Returns whether this machine uses big or little endian.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:348
QgsWKBTypes::Type readHeader() const
Definition: qgswkbptr.cpp:37
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...
void combineExtentWith(QgsRectangle *rect)
expand the rectangle so that covers both the original rectangle and the given rectangle ...
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:202
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:187
static QgsRectangle calculateBoundingBox(QGis::WkbType wkbType, QgsConstWkbPtr wkbPtr, int numPoints)
Returns the BBOX of the specified WKB-point stream.
The geometries can be fully simplified by its BoundingBox.
QgsMapToPixelSimplifier(int simplifyFlags, double tolerance)
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.
QGis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
static WkbType singleType(WkbType type)
Map multi to single type.
Definition: qgis.cpp:353
double mTolerance
Distance tolerance for the simplification.
QString what() const
Definition: qgsexception.h:36
int wkbSize() const
Returns the size of the WKB in asWkb().
static bool generalizeWkbGeometryByBoundingBox(QGis::WkbType wkbType, QgsConstWkbPtr sourceWkbPtr, QgsWkbPtr targetWkbPtr, int &targetWkbSize, const QgsRectangle &envelope, bool writeHeader)
Generalize the WKB-geometry using the BBOX of the original geometry.
static Type flatType(Type type)
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:366
int mSimplifyFlags
Current simplification flags.
The geometries can be simplified using the current map2pixel context state.
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:207
const unsigned char * asWkb() const
Returns the buffer containing this geometry in WKB format.
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:192
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212