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
20 #include "qgsmaptopixelgeometrysimplifier.h"
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
QgsMapToPixelSimplifier(int simplifyFlags,double tolerance,SimplifyAlgorithm simplifyAlgorithm)30 QgsMapToPixelSimplifier::QgsMapToPixelSimplifier( int simplifyFlags, double tolerance, SimplifyAlgorithm simplifyAlgorithm )
31 : mSimplifyFlags( simplifyFlags )
32 , mSimplifyAlgorithm( simplifyAlgorithm )
33 , mTolerance( tolerance )
34 {
35 }
36
37 //////////////////////////////////////////////////////////////////////////////////////////////
38 // Helper simplification methods
39
calculateLengthSquared2D(double x1,double y1,double x2,double y2)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
equalSnapToGrid(double x1,double y1,double x2,double y2,double gridOriginX,double gridOriginY,float gridInverseSizeXY)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
59 //////////////////////////////////////////////////////////////////////////////////////////////
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
68 //////////////////////////////////////////////////////////////////////////////////////////////
69
70 //! Generalize the WKB-geometry using the BBOX of the original geometry
generalizeWkbGeometryByBoundingBox(QgsWkbTypes::Type wkbType,const QgsAbstractGeometry & geometry,const QgsRectangle & envelope,bool isRing)71 static std::unique_ptr< QgsAbstractGeometry > generalizeWkbGeometryByBoundingBox(
72 QgsWkbTypes::Type wkbType,
73 const QgsAbstractGeometry &geometry,
74 const QgsRectangle &envelope,
75 bool isRing )
76 {
77 unsigned int geometryType = QgsWkbTypes::singleType( QgsWkbTypes::flatType( wkbType ) );
78
79 // If the geometry is already minimal skip the generalization
80 int minimumSize = geometryType == QgsWkbTypes::LineString ? 2 : 5;
81
82 if ( geometry.nCoordinates() <= minimumSize )
83 {
84 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
85 }
86
87 const double x1 = envelope.xMinimum();
88 const double y1 = envelope.yMinimum();
89 const double x2 = envelope.xMaximum();
90 const double y2 = envelope.yMaximum();
91
92 // Write the generalized geometry
93 if ( geometryType == QgsWkbTypes::LineString && !isRing )
94 {
95 return qgis::make_unique< QgsLineString >( QVector<double>() << x1 << x2, QVector<double>() << y1 << y2 );
96 }
97 else
98 {
99 std::unique_ptr< QgsLineString > ext = qgis::make_unique< QgsLineString >(
100 QVector< double >() << x1
101 << x2
102 << x2
103 << x1
104 << x1,
105 QVector< double >() << y1
106 << y1
107 << y2
108 << y2
109 << y1 );
110 if ( geometryType == QgsWkbTypes::LineString )
111 return std::move( ext );
112 else
113 {
114 std::unique_ptr< QgsPolygon > polygon = qgis::make_unique< QgsPolygon >();
115 polygon->setExteriorRing( ext.release() );
116 return std::move( polygon );
117 }
118 }
119 }
120
simplifyGeometry(int simplifyFlags,SimplifyAlgorithm simplifyAlgorithm,const QgsAbstractGeometry & geometry,double map2pixelTol,bool isaLinearRing)121 std::unique_ptr< QgsAbstractGeometry > QgsMapToPixelSimplifier::simplifyGeometry( int simplifyFlags,
122 SimplifyAlgorithm simplifyAlgorithm,
123 const QgsAbstractGeometry &geometry, double map2pixelTol,
124 bool isaLinearRing )
125 {
126 bool isGeneralizable = true;
127 QgsWkbTypes::Type wkbType = geometry.wkbType();
128
129 // Can replace the geometry by its BBOX ?
130 QgsRectangle envelope = geometry.boundingBox();
131 if ( ( simplifyFlags & QgsMapToPixelSimplifier::SimplifyEnvelope ) &&
132 isGeneralizableByMapBoundingBox( envelope, map2pixelTol ) )
133 {
134 return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope, isaLinearRing );
135 }
136
137 if ( !( simplifyFlags & QgsMapToPixelSimplifier::SimplifyGeometry ) )
138 isGeneralizable = false;
139
140 const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( wkbType );
141
142 // Write the geometry
143 if ( flatType == QgsWkbTypes::LineString || flatType == QgsWkbTypes::CircularString )
144 {
145 const QgsCurve &srcCurve = dynamic_cast<const QgsCurve &>( geometry );
146 const int numPoints = srcCurve.numPoints();
147
148 std::unique_ptr<QgsCurve> output;
149
150 QVector< double > lineStringX;
151 QVector< double > lineStringY;
152 if ( flatType == QgsWkbTypes::LineString )
153 {
154 // if we are making a linestring, we do it in an optimised way by directly constructing
155 // the final x/y vectors, which avoids calling the slower insertVertex method
156 lineStringX.reserve( numPoints );
157 lineStringY.reserve( numPoints );
158 }
159 else
160 {
161 output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.createEmptyWithSameType() ) );
162 }
163
164 double x = 0.0, y = 0.0, lastX = 0.0, lastY = 0.0;
165
166 if ( numPoints <= ( isaLinearRing ? 4 : 2 ) )
167 isGeneralizable = false;
168
169 bool isLongSegment;
170 bool hasLongSegments = false; //-> To avoid replace the simplified geometry by its BBOX when there are 'long' segments.
171
172 // Check whether the LinearRing is really closed.
173 if ( isaLinearRing )
174 {
175 isaLinearRing = qgsDoubleNear( srcCurve.xAt( 0 ), srcCurve.xAt( numPoints - 1 ) ) &&
176 qgsDoubleNear( srcCurve.yAt( 0 ), srcCurve.yAt( numPoints - 1 ) );
177 }
178
179 // Process each vertex...
180 switch ( simplifyAlgorithm )
181 {
182 case SnapToGrid:
183 {
184 double gridOriginX = envelope.xMinimum();
185 double gridOriginY = envelope.yMinimum();
186
187 // Use a factor for the maximum displacement distance for simplification, similar as GeoServer does
188 float gridInverseSizeXY = map2pixelTol != 0 ? ( float )( 1.0f / ( 0.8 * map2pixelTol ) ) : 0.0f;
189
190 const double *xData = nullptr;
191 const double *yData = nullptr;
192 if ( flatType == QgsWkbTypes::LineString )
193 {
194 xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
195 yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
196 }
197
198 for ( int i = 0; i < numPoints; ++i )
199 {
200 if ( xData && yData )
201 {
202 x = *xData++;
203 y = *yData++;
204 }
205 else
206 {
207 x = srcCurve.xAt( i );
208 y = srcCurve.yAt( i );
209 }
210
211 if ( i == 0 ||
212 !isGeneralizable ||
213 !equalSnapToGrid( x, y, lastX, lastY, gridOriginX, gridOriginY, gridInverseSizeXY ) ||
214 ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
215 {
216 if ( output )
217 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
218 else
219 {
220 lineStringX.append( x );
221 lineStringY.append( y );
222 }
223 lastX = x;
224 lastY = y;
225 }
226 }
227 break;
228 }
229
230 case SnappedToGridGlobal:
231 {
232 output.reset( qgsgeometry_cast< QgsCurve * >( srcCurve.snappedToGrid( map2pixelTol, map2pixelTol ) ) );
233 break;
234 }
235
236 case Visvalingam:
237 {
238 map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'Area' calculations.
239
240 EFFECTIVE_AREAS ea( srcCurve );
241
242 int set_area = 0;
243 ptarray_calc_areas( &ea, isaLinearRing ? 4 : 2, set_area, map2pixelTol );
244
245 for ( int i = 0; i < numPoints; ++i )
246 {
247 if ( ea.res_arealist[ i ] > map2pixelTol )
248 {
249 if ( output )
250 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), ea.inpts.at( i ) );
251 else
252 {
253 lineStringX.append( ea.inpts.at( i ).x() );
254 lineStringY.append( ea.inpts.at( i ).y() );
255 }
256 }
257 }
258 break;
259 }
260
261 case Distance:
262 {
263 map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.
264
265 const double *xData = nullptr;
266 const double *yData = nullptr;
267 if ( flatType == QgsWkbTypes::LineString )
268 {
269 xData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->xData();
270 yData = qgsgeometry_cast< const QgsLineString * >( &srcCurve )->yData();
271 }
272
273 for ( int i = 0; i < numPoints; ++i )
274 {
275 if ( xData && yData )
276 {
277 x = *xData++;
278 y = *yData++;
279 }
280 else
281 {
282 x = srcCurve.xAt( i );
283 y = srcCurve.yAt( i );
284 }
285
286 isLongSegment = false;
287
288 if ( i == 0 ||
289 !isGeneralizable ||
290 ( isLongSegment = ( calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol ) ) ||
291 ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
292 {
293 if ( output )
294 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( x, y ) );
295 else
296 {
297 lineStringX.append( x );
298 lineStringY.append( y );
299 }
300 lastX = x;
301 lastY = y;
302
303 hasLongSegments |= isLongSegment;
304 }
305 }
306 }
307 }
308
309 if ( !output )
310 {
311 output = qgis::make_unique< QgsLineString >( lineStringX, lineStringY );
312 }
313 if ( output->numPoints() < ( isaLinearRing ? 4 : 2 ) )
314 {
315 // we simplified the geometry too much!
316 if ( !hasLongSegments )
317 {
318 // approximate the geometry's shape by its bounding box
319 // (rect for linear ring / one segment for line string)
320 return generalizeWkbGeometryByBoundingBox( wkbType, geometry, envelope, isaLinearRing );
321 }
322 else
323 {
324 // Bad luck! The simplified geometry is invalid and approximation by bounding box
325 // would create artifacts due to long segments.
326 // We will return the original geometry
327 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
328 }
329 }
330
331 if ( isaLinearRing )
332 {
333 // make sure we keep the linear ring closed
334 if ( !qgsDoubleNear( lastX, output->xAt( 0 ) ) || !qgsDoubleNear( lastY, output->yAt( 0 ) ) )
335 {
336 output->insertVertex( QgsVertexId( 0, 0, output->numPoints() ), QgsPoint( output->xAt( 0 ), output->yAt( 0 ) ) );
337 }
338 }
339
340 return std::move( output );
341 }
342 else if ( flatType == QgsWkbTypes::Polygon )
343 {
344 const QgsPolygon &srcPolygon = dynamic_cast<const QgsPolygon &>( geometry );
345 std::unique_ptr<QgsPolygon> polygon( new QgsPolygon() );
346 std::unique_ptr<QgsAbstractGeometry> extRing = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *srcPolygon.exteriorRing(), map2pixelTol, true );
347 polygon->setExteriorRing( qgsgeometry_cast<QgsCurve *>( extRing.release() ) );
348 for ( int i = 0; i < srcPolygon.numInteriorRings(); ++i )
349 {
350 const QgsCurve *sub = srcPolygon.interiorRing( i );
351 std::unique_ptr< QgsAbstractGeometry > ring = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, true );
352 polygon->addInteriorRing( qgsgeometry_cast<QgsCurve *>( ring.release() ) );
353 }
354 return std::move( polygon );
355 }
356 else if ( QgsWkbTypes::isMultiType( flatType ) )
357 {
358 const QgsGeometryCollection &srcCollection = dynamic_cast<const QgsGeometryCollection &>( geometry );
359 std::unique_ptr<QgsGeometryCollection> collection( srcCollection.createEmptyWithSameType() );
360 const int numGeoms = srcCollection.numGeometries();
361 collection->reserve( numGeoms );
362 for ( int i = 0; i < numGeoms; ++i )
363 {
364 const QgsAbstractGeometry *sub = srcCollection.geometryN( i );
365 std::unique_ptr< QgsAbstractGeometry > part = simplifyGeometry( simplifyFlags, simplifyAlgorithm, *sub, map2pixelTol, false );
366 collection->addGeometry( part.release() );
367 }
368 return std::move( collection );
369 }
370 return std::unique_ptr< QgsAbstractGeometry >( geometry.clone() );
371 }
372
373 //////////////////////////////////////////////////////////////////////////////////////////////
374
isGeneralizableByMapBoundingBox(const QgsRectangle & envelope,double map2pixelTol)375 bool QgsMapToPixelSimplifier::isGeneralizableByMapBoundingBox( const QgsRectangle &envelope, double map2pixelTol )
376 {
377 // Can replace the geometry by its BBOX ?
378 return envelope.width() < map2pixelTol && envelope.height() < map2pixelTol;
379 }
380
simplify(const QgsGeometry & geometry) const381 QgsGeometry QgsMapToPixelSimplifier::simplify( const QgsGeometry &geometry ) const
382 {
383 if ( geometry.isNull() )
384 {
385 return QgsGeometry();
386 }
387 if ( mSimplifyFlags == QgsMapToPixelSimplifier::NoFlags )
388 {
389 return geometry;
390 }
391
392 // Check whether the geometry can be simplified using the map2pixel context
393 const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry.wkbType() );
394 const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
395 if ( flatType == QgsWkbTypes::Point )
396 {
397 return geometry;
398 }
399
400 const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
401 const int numPoints = geometry.constGet()->nCoordinates();
402
403 if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
404 {
405 // No simplify simple geometries
406 return geometry;
407 }
408
409 const QgsRectangle envelope = geometry.boundingBox();
410 if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
411 {
412 //points are in average too far apart to lead to any significant simplification
413 return geometry;
414 }
415
416 return QgsGeometry( simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry.constGet(), mTolerance, false ) );
417 }
418
simplify(const QgsAbstractGeometry * geometry) const419 QgsAbstractGeometry *QgsMapToPixelSimplifier::simplify( const QgsAbstractGeometry *geometry ) const
420 {
421 //
422 // IMPORTANT!!!!!!!
423 // We want to avoid any geometry cloning we possibly can here, which is why the
424 // "fail" paths always return nullptr
425 //
426
427 if ( !geometry )
428 {
429 return nullptr;
430 }
431 if ( mSimplifyFlags == QgsMapToPixelSimplifier::NoFlags )
432 {
433 return nullptr;
434 }
435
436 // Check whether the geometry can be simplified using the map2pixel context
437 const QgsWkbTypes::Type singleType = QgsWkbTypes::singleType( geometry->wkbType() );
438 const QgsWkbTypes::Type flatType = QgsWkbTypes::flatType( singleType );
439 if ( flatType == QgsWkbTypes::Point )
440 {
441 return nullptr;
442 }
443
444 const bool isaLinearRing = flatType == QgsWkbTypes::Polygon;
445 const int numPoints = geometry->nCoordinates();
446
447 if ( numPoints <= ( isaLinearRing ? 6 : 3 ) )
448 {
449 // No simplify simple geometries
450 return nullptr;
451 }
452
453 const QgsRectangle envelope = geometry->boundingBox();
454 if ( std::max( envelope.width(), envelope.height() ) / numPoints > mTolerance * 2.0 )
455 {
456 //points are in average too far apart to lead to any significant simplification
457 return nullptr;
458 }
459
460 return simplifyGeometry( mSimplifyFlags, mSimplifyAlgorithm, *geometry, mTolerance, false ).release();
461 }
462