QGIS API Documentation  2.13.0-Master
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 <stdexcept>
20 #include "qgsapplication.h"
21 #include "qgslogger.h"
22 
23 class QgsParserException: public std::runtime_error
24 {
25  public:
27  : std::runtime_error( msg.toStdString() )
28  {}
29 };
30 
32 {
33  public:
35  : QgsParserException( QString( "Premature end of WKB: " ) + msg )
36  {}
37 };
38 
39 QgsMapToPixelSimplifier::QgsMapToPixelSimplifier( int simplifyFlags, double tolerance )
40  : mSimplifyFlags( simplifyFlags )
41  , mTolerance( tolerance )
42 {
43 }
44 
46 {
47 }
48 
50 // Helper simplification methods
51 
53 float QgsMapToPixelSimplifier::calculateLengthSquared2D( double x1, double y1, double x2, double y2 )
54 {
55  float vx = static_cast< float >( x2 - x1 );
56  float vy = static_cast< float >( y2 - y1 );
57 
58  return vx*vx + vy*vy;
59 }
60 
62 inline static QgsRectangle calculateBoundingBox( QGis::WkbType wkbType, const unsigned char* wkb, int numPoints )
63 {
64  double x, y;
65  QgsRectangle r;
66  r.setMinimal();
67 
68  int sizeOfDoubleX = sizeof( double );
69  int sizeOfDoubleY = ( QGis::wkbDimensions( wkbType ) - 1 ) * sizeof( double );
70 
71  for ( int i = 0; i < numPoints; ++i )
72  {
73  memcpy( &x, wkb, sizeof( double ) );
74  wkb += sizeOfDoubleX;
75  memcpy( &y, wkb, sizeof( double ) );
76  wkb += sizeOfDoubleY;
77  r.combineExtentWith( x, y );
78  }
79 
80  return r;
81 }
82 
85  QGis::WkbType wkbType,
86  const unsigned char *sourceWkb, int sourceWkbSize,
87  unsigned char *targetWkb, int &targetWkbSize,
88  const QgsRectangle &envelope, bool writeHeader )
89 {
90  Q_UNUSED( sourceWkb );
91  unsigned char* wkb2 = targetWkb;
92  unsigned int geometryType = QGis::singleType( QGis::flatType( wkbType ) );
93 
94  int sizeOfDoubleX = sizeof( double );
95  int sizeOfDoubleY = sizeof( double ) * ( QGis::wkbDimensions( wkbType ) - 1 );
96 
97  // If the geometry is already minimal skip the generalization
98  int minimumSize = geometryType == QGis::WKBLineString ? 4 + 2 * ( sizeOfDoubleX + sizeOfDoubleY ) : 8 + 5 * ( sizeOfDoubleX + sizeOfDoubleY );
99 
100  if ( writeHeader )
101  minimumSize += 5;
102 
103  if ( sourceWkbSize <= minimumSize )
104  {
105  targetWkbSize = 0;
106  return false;
107  }
108 
109  double x1 = envelope.xMinimum();
110  double y1 = envelope.yMinimum();
111  double x2 = envelope.xMaximum();
112  double y2 = envelope.yMaximum();
113 
114  // Write the main header of the geometry
115  if ( writeHeader )
116  {
117  char byteOrder = QgsApplication::endian(); // byteOrder
118  memcpy( targetWkb, &byteOrder, 1 );
119  targetWkb += 1;
120 
121  memcpy( targetWkb, &geometryType, 4 ); // type
122  targetWkb += 4;
123 
124  if ( geometryType == QGis::WKBPolygon ) // numRings
125  {
126  int numRings = 1;
127  memcpy( targetWkb, &numRings, 4 );
128  targetWkb += 4;
129  }
130  }
131 
132  // Write the generalized geometry
133  if ( geometryType == QGis::WKBLineString )
134  {
135  int numPoints = 2;
136  memcpy( targetWkb, &numPoints, 4 ); // numPoints;
137  targetWkb += 4;
138 
139  memcpy( targetWkb, &x1, sizeof( double ) );
140  targetWkb += sizeof( double );
141  memcpy( targetWkb, &y1, sizeof( double ) );
142  targetWkb += sizeof( double );
143 
144  memcpy( targetWkb, &x2, sizeof( double ) );
145  targetWkb += sizeof( double );
146  memcpy( targetWkb, &y2, sizeof( double ) );
147  targetWkb += sizeof( double );
148  }
149  else
150  {
151  int numPoints = 5;
152  memcpy( targetWkb, &numPoints, 4 ); // numPoints;
153  targetWkb += 4;
154 
155  memcpy( targetWkb, &x1, sizeof( double ) );
156  targetWkb += sizeof( double );
157  memcpy( targetWkb, &y1, sizeof( double ) );
158  targetWkb += sizeof( double );
159 
160  memcpy( targetWkb, &x2, sizeof( double ) );
161  targetWkb += sizeof( double );
162  memcpy( targetWkb, &y1, sizeof( double ) );
163  targetWkb += sizeof( double );
164 
165  memcpy( targetWkb, &x2, sizeof( double ) );
166  targetWkb += sizeof( double );
167  memcpy( targetWkb, &y2, sizeof( double ) );
168  targetWkb += sizeof( double );
169 
170  memcpy( targetWkb, &x1, sizeof( double ) );
171  targetWkb += sizeof( double );
172  memcpy( targetWkb, &y2, sizeof( double ) );
173  targetWkb += sizeof( double );
174 
175  memcpy( targetWkb, &x1, sizeof( double ) );
176  targetWkb += sizeof( double );
177  memcpy( targetWkb, &y1, sizeof( double ) );
178  targetWkb += sizeof( double );
179  }
180  targetWkbSize += targetWkb - wkb2;
181 
182  return true;
183 }
184 
186 bool QgsMapToPixelSimplifier::simplifyWkbGeometry(
187  int simplifyFlags, QGis::WkbType wkbType,
188  const unsigned char *sourceWkb, int sourceWkbSize,
189  unsigned char *targetWkb, int &targetWkbSize,
190  const QgsRectangle &envelope, double map2pixelTol,
191  bool writeHeader, bool isaLinearRing )
192 {
193  bool isGeneralizable = true;
194  bool result = false;
195 
196  // Save initial WKB settings to use when the simplification creates invalid geometries
197  const unsigned char* sourcePrevWkb = sourceWkb;
198  unsigned char* targetPrevWkb = targetWkb;
199  int targetWkbPrevSize = targetWkbSize;
200 
201  const unsigned char* endOfSourceWkb = sourceWkb + sourceWkbSize;
202 
203  // Can replace the geometry by its BBOX ?
205  isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
206  {
207  isGeneralizable = generalizeWkbGeometryByBoundingBox( wkbType, sourceWkb, sourceWkbSize, targetWkb, targetWkbSize, envelope, writeHeader );
208  if ( isGeneralizable )
209  return true;
210  }
211 
213  isGeneralizable = false;
214 
215  // Write the main header of the geometry
216  if ( writeHeader )
217  {
218  if ( sourceWkbSize < 5 )
219  throw QgsParserException( QString( "Premature end of WKB reading header " ) );
220 
221  targetWkb[0] = sourceWkb[0]; // byteOrder
222  sourceWkb += 1;
223  targetWkb += 1;
224 
225  int geometryType;
226  memcpy( &geometryType, sourceWkb, 4 );
227  int flatType = QGis::flatType( static_cast< QGis::WkbType >( geometryType ) );
228  memcpy( targetWkb, &flatType, 4 ); // type
229  sourceWkb += 4;
230  targetWkb += 4;
231 
232  targetWkbSize += 5;
233  }
234 
235  const unsigned char* wkb1 = sourceWkb;
236  unsigned char* wkb2 = targetWkb;
237  unsigned int flatType = QGis::flatType( wkbType );
238 
239  // Write the geometry
240  if ( flatType == QGis::WKBLineString || isaLinearRing )
241  {
242  double x, y, lastX = 0, lastY = 0;
243  QgsRectangle r;
244  r.setMinimal();
245 
246  int sizeOfDoubleX = sizeof( double );
247  int sizeOfDoubleY = ( QGis::wkbDimensions( wkbType ) - 1 ) * sizeof( double );
248 
249  if ( sourceWkb + 4 >= endOfSourceWkb )
250  throw QgsShortWkbException( "reading numPoints" );
251 
252  int numPoints;
253  memcpy( &numPoints, sourceWkb, 4 );
254  sourceWkb += 4;
255  if ( numPoints <= ( isaLinearRing ? 5 : 2 ) )
256  isGeneralizable = false;
257 
258  int numTargetPoints = 0;
259  memcpy( targetWkb, &numTargetPoints, 4 );
260  targetWkb += 4;
261  targetWkbSize += 4;
262 
263  double* ptr = reinterpret_cast< double* >( targetWkb );
264  map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
265 
266  bool isLongSegment;
267  bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
268 
269  // Check whether the LinearRing is really closed.
270  if ( isaLinearRing )
271  {
272  double x1, y1, x2, y2;
273 
274  const unsigned char* startWkbX = sourceWkb;
275  const unsigned char* startWkbY = startWkbX + sizeOfDoubleX;
276  const unsigned char* finalWkbX = sourceWkb + ( numPoints - 1 ) * ( sizeOfDoubleX + sizeOfDoubleY );
277  const unsigned char* finalWkbY = finalWkbX + sizeOfDoubleX;
278 
279  if ( finalWkbY + sizeof( double ) > endOfSourceWkb )
280  throw QgsShortWkbException( "reading last point" );
281 
282 
283  memcpy( &x1, startWkbX, sizeof( double ) );
284  memcpy( &y1, startWkbY, sizeof( double ) );
285  memcpy( &x2, finalWkbX, sizeof( double ) );
286  memcpy( &y2, finalWkbY, sizeof( double ) );
287 
288  isaLinearRing = qgsDoubleNear( x1, x2 ) && qgsDoubleNear( y1, y2 );
289  }
290 
291  // Process each vertex...
292  for ( int i = 0; i < numPoints; ++i )
293  {
294  if ( sourceWkb + sizeOfDoubleX + sizeOfDoubleY > endOfSourceWkb )
295  {
296  throw QgsParserException( QString( "Premature end of WKB reading point %1/%2" ) .arg( i + 1 ) .arg( numPoints ) );
297  }
298 
299  memcpy( &x, sourceWkb, sizeof( double ) );
300  sourceWkb += sizeOfDoubleX;
301  memcpy( &y, sourceWkb, sizeof( double ) );
302  sourceWkb += sizeOfDoubleY;
303 
304  isLongSegment = false;
305 
306  if ( i == 0 ||
307  !isGeneralizable ||
308  ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
309  ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
310  {
311  memcpy( ptr, &x, sizeof( double ) );
312  lastX = x;
313  ptr++;
314  memcpy( ptr, &y, sizeof( double ) );
315  lastY = y;
316  ptr++;
317  numTargetPoints++;
318 
319  hasLongSegments |= isLongSegment;
320  }
321 
322  r.combineExtentWith( x, y );
323  }
324  targetWkb = wkb2 + 4;
325 
326  if ( numTargetPoints < ( isaLinearRing ? 4 : 2 ) )
327  {
328  // we simplified the geometry too much!
329  if ( !hasLongSegments )
330  {
331  // approximate the geometry's shape by its bounding box
332  // (rect for linear ring / one segment for line string)
333  unsigned char* targetTempWkb = targetWkb;
334  int targetWkbTempSize = targetWkbSize;
335 
336  sourceWkb = sourcePrevWkb;
337  targetWkb = targetPrevWkb;
338  targetWkbSize = targetWkbPrevSize;
339  if ( generalizeWkbGeometryByBoundingBox( wkbType, sourceWkb, sourceWkbSize, targetWkb, targetWkbSize, r, writeHeader ) )
340  return true;
341 
342  targetWkb = targetTempWkb;
343  targetWkbSize = targetWkbTempSize;
344  }
345  else
346  {
347  // Bad luck! The simplified geometry is invalid and approximation by bounding box
348  // would create artifacts due to long segments. Worst of all, we may have overwritten
349  // the original coordinates by the simplified ones (source and target WKB ptr can be the same)
350  // so we cannot even undo the changes here. We will return invalid geometry and hope that
351  // other pieces of QGIS will survive that :-/
352  }
353  }
354  if ( isaLinearRing )
355  {
356  // make sure we keep the linear ring closed
357  memcpy( &x, targetWkb + 0, sizeof( double ) );
358  memcpy( &y, targetWkb + sizeof( double ), sizeof( double ) );
359  if ( !qgsDoubleNear( lastX, x ) || !qgsDoubleNear( lastY, y ) )
360  {
361  memcpy( ptr, &x, sizeof( double ) );
362  ptr++;
363  memcpy( ptr, &y, sizeof( double ) );
364  ptr++;
365  numTargetPoints++;
366  }
367  }
368  targetWkbSize += numTargetPoints * sizeof( double ) * 2;
369  targetWkb = wkb2;
370 
371  memcpy( targetWkb, &numTargetPoints, 4 );
372  result = numPoints != numTargetPoints;
373  }
374  else if ( flatType == QGis::WKBPolygon )
375  {
376  int numRings;
377  memcpy( &numRings, sourceWkb, 4 );
378  sourceWkb += 4;
379 
380  memcpy( targetWkb, &numRings, 4 );
381  targetWkb += 4;
382  targetWkbSize += 4;
383 
384  for ( int i = 0; i < numRings; ++i )
385  {
386  int numPoints_i;
387  memcpy( &numPoints_i, sourceWkb, 4 );
388  QgsRectangle envelope_i = numRings == 1 ? envelope : calculateBoundingBox( wkbType, sourceWkb + 4, numPoints_i );
389 
390  int sourceWkbSize_i = 4 + numPoints_i * QGis::wkbDimensions( wkbType ) * sizeof( double );
391  int targetWkbSize_i = 0;
392 
393  result |= simplifyWkbGeometry( simplifyFlags, wkbType, sourceWkb, sourceWkbSize_i, targetWkb, targetWkbSize_i, envelope_i, map2pixelTol, false, true );
394  sourceWkb += sourceWkbSize_i;
395  targetWkb += targetWkbSize_i;
396 
397  targetWkbSize += targetWkbSize_i;
398  }
399  }
400  else if ( flatType == QGis::WKBMultiLineString || flatType == QGis::WKBMultiPolygon )
401  {
402  int numGeoms;
403  memcpy( &numGeoms, sourceWkb, 4 );
404  sourceWkb += 4;
405  wkb1 += 4;
406 
407  memcpy( targetWkb, &numGeoms, 4 );
408  targetWkb += 4;
409  targetWkbSize += 4;
410 
411  for ( int i = 0; i < numGeoms; ++i )
412  {
413  int sourceWkbSize_i = 0;
414  int targetWkbSize_i = 0;
415 
416  // ... calculate the wkb-size of the current child complex geometry
417  if ( flatType == QGis::WKBMultiLineString )
418  {
419  int numPoints_i;
420  memcpy( &numPoints_i, wkb1 + 5, 4 );
421  int wkbSize_i = 4 + numPoints_i * QGis::wkbDimensions( wkbType ) * sizeof( double );
422 
423  sourceWkbSize_i += 5 + wkbSize_i;
424  wkb1 += 5 + wkbSize_i;
425  }
426  else
427  {
428  int numPrings_i;
429  memcpy( &numPrings_i, wkb1 + 5, 4 );
430  sourceWkbSize_i = 9;
431  wkb1 += 9;
432 
433  for ( int j = 0; j < numPrings_i; ++j )
434  {
435  int numPoints_i;
436  memcpy( &numPoints_i, wkb1, 4 );
437  int wkbSize_i = 4 + numPoints_i * QGis::wkbDimensions( wkbType ) * sizeof( double );
438 
439  sourceWkbSize_i += wkbSize_i;
440  wkb1 += wkbSize_i;
441  }
442  }
443  result |= simplifyWkbGeometry( simplifyFlags, QGis::singleType( wkbType ), sourceWkb, endOfSourceWkb - sourceWkb, targetWkb, targetWkbSize_i, envelope, map2pixelTol, true, false );
444  sourceWkb += sourceWkbSize_i;
445  targetWkb += targetWkbSize_i;
446 
447  targetWkbSize += targetWkbSize_i;
448  }
449  }
450  return result;
451 }
452 
454 
456 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle& envelope, double map2pixelTol )
457 {
458  // Can replace the geometry by its BBOX ?
459  return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
460 }
461 
464 {
465  QgsGeometry* g = new QgsGeometry();
466 
467  int wkbSize = geometry->wkbSize();
468  unsigned char* wkb = reinterpret_cast< unsigned char* >( malloc( wkbSize ) );
469  memcpy( wkb, geometry->asWkb(), wkbSize );
470  g->fromWkb( wkb, wkbSize );
472 
473  return g;
474 }
475 
478 {
479  int finalWkbSize = 0;
480 
481  // Check whether the geometry can be simplified using the map2pixel context
482  QGis::GeometryType geometryType = geometry->type();
483  if ( !( geometryType == QGis::Line || geometryType == QGis::Polygon ) )
484  return false;
485 
486  QgsRectangle envelope = geometry->boundingBox();
487  QGis::WkbType wkbType = geometry->wkbType();
488 
489  const unsigned char* wkb = geometry->asWkb();
490  int wkbSize = geometry->wkbSize();
491 
492  unsigned char* targetWkb = new unsigned char[wkbSize];
493  memcpy( targetWkb, wkb, wkbSize );
494 
495  try
496  {
497  if ( simplifyWkbGeometry( simplifyFlags, wkbType, wkb, wkbSize, targetWkb, finalWkbSize, envelope, tolerance ) )
498  {
499  unsigned char* finalWkb = new unsigned char[finalWkbSize];
500  memcpy( finalWkb, targetWkb, finalWkbSize );
501  geometry->fromWkb( finalWkb, finalWkbSize );
502  delete [] targetWkb;
503  return true;
504  }
505  }
506  catch ( const QgsParserException &e )
507  {
508  QgsDebugMsg( QString( "Exception thrown by simplifier: %1" ) .arg( e.what() ) );
509  }
510  delete [] targetWkb;
511  return false;
512 }
513 
516 {
517  return simplifyGeometry( geometry, mSimplifyFlags, mTolerance );
518 }
static WkbType flatType(WkbType type)
Map 2d+ to 2d type.
Definition: qgis.cpp:352
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:383
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:196
#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)
Definition: qgis.h:285
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:201
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:186
The geometries can be fully simplified by its BoundingBox.
QgsMapToPixelSimplifier(int simplifyFlags, double tolerance)
static QgsRectangle calculateBoundingBox(QGis::WkbType wkbType, const unsigned char *wkb, int numPoints)
Returns the BBOX of the specified WKB-point stream.
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:310
double mTolerance
Distance tolerance for the simplification.
static bool generalizeWkbGeometryByBoundingBox(QGis::WkbType wkbType, const unsigned char *sourceWkb, int sourceWkbSize, unsigned char *targetWkb, int &targetWkbSize, const QgsRectangle &envelope, bool writeHeader)
Generalize the WKB-geometry using the BBOX of the original geometry.
int wkbSize() const
Returns the size of the WKB in asWkb().
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:206
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:191
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:211