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