QGIS API Documentation  2.13.0-Master
qgstininterpolator.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgstininterpolator.cpp
3  ----------------------
4  begin : March 10, 2008
5  copyright : (C) 2008 by Marco Hugentobler
6  email : marco dot hugentobler at karto dot baug dot ethz dot ch
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgstininterpolator.h"
20 #include "DualEdgeTriangulation.h"
21 #include "NormVecDecorator.h"
23 #include "Point3D.h"
24 #include "qgsfeature.h"
25 #include "qgsgeometry.h"
26 #include "qgsvectorlayer.h"
27 #include "qgswkbptr.h"
28 #include <QProgressDialog>
29 
30 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData>& inputData, TIN_INTERPOLATION interpolation, bool showProgressDialog )
31  : QgsInterpolator( inputData )
32  , mTriangulation( nullptr )
33  , mTriangleInterpolator( nullptr )
34  , mIsInitialized( false )
35  , mShowProgressDialog( showProgressDialog )
36  , mExportTriangulationToFile( false )
37  , mInterpolation( interpolation )
38 {
39 }
40 
42 {
43  delete mTriangulation;
44  delete mTriangleInterpolator;
45 }
46 
47 int QgsTINInterpolator::interpolatePoint( double x, double y, double& result )
48 {
49  if ( !mIsInitialized )
50  {
51  initialize();
52  }
53 
54  if ( !mTriangleInterpolator )
55  {
56  return 1;
57  }
58 
59  Point3D r;
60  if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
61  {
62  return 2;
63  }
64  result = r.getZ();
65  return 0;
66 }
67 
68 void QgsTINInterpolator::initialize()
69 {
70  DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( 100000, nullptr );
71  if ( mInterpolation == CloughTocher )
72  {
74  dec->addTriangulation( theDualEdgeTriangulation );
75  mTriangulation = dec;
76  }
77  else
78  {
79  mTriangulation = theDualEdgeTriangulation;
80  }
81 
82  //get number of features if we use a progress bar
83  int nFeatures = 0;
84  int nProcessedFeatures = 0;
85  if ( mShowProgressDialog )
86  {
87  Q_FOREACH ( const LayerData& layer, mLayerData )
88  {
89  if ( layer.vectorLayer )
90  {
91  nFeatures += layer.vectorLayer->featureCount();
92  }
93  }
94  }
95 
96  QProgressDialog* theProgressDialog = nullptr;
97  if ( mShowProgressDialog )
98  {
99  theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, nullptr );
100  theProgressDialog->setWindowModality( Qt::WindowModal );
101  }
102 
103 
104  QgsFeature f;
105  Q_FOREACH ( const LayerData& layer, mLayerData )
106  {
107  if ( layer.vectorLayer )
108  {
109  QgsAttributeList attList;
110  if ( !layer.zCoordInterpolation )
111  {
112  attList.push_back( layer.interpolationAttribute );
113  }
114 
115  QgsFeatureIterator fit = layer.vectorLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ) );
116 
117  while ( fit.nextFeature( f ) )
118  {
119  if ( mShowProgressDialog )
120  {
121  if ( theProgressDialog->wasCanceled() )
122  {
123  break;
124  }
125  theProgressDialog->setValue( nProcessedFeatures );
126  }
127  insertData( &f, layer.zCoordInterpolation, layer.interpolationAttribute, layer.mInputType );
128  ++nProcessedFeatures;
129  }
130  }
131  }
132 
133  delete theProgressDialog;
134 
135  if ( mInterpolation == CloughTocher )
136  {
137  CloughTocherInterpolator* ctInterpolator = new CloughTocherInterpolator();
138  NormVecDecorator* dec = dynamic_cast<NormVecDecorator*>( mTriangulation );
139  if ( dec )
140  {
141  QProgressDialog* progressDialog = nullptr;
142  if ( mShowProgressDialog ) //show a progress dialog because it can take a long time...
143  {
144  progressDialog = new QProgressDialog();
145  progressDialog->setLabelText( QObject::tr( "Estimating normal derivatives..." ) );
146  }
147  dec->estimateFirstDerivatives( progressDialog );
148  delete progressDialog;
149  ctInterpolator->setTriangulation( dec );
150  dec->setTriangleInterpolator( ctInterpolator );
151  mTriangleInterpolator = ctInterpolator;
152  }
153  }
154  else //linear
155  {
156  mTriangleInterpolator = new LinTriangleInterpolator( theDualEdgeTriangulation );
157  }
158  mIsInitialized = true;
159 
160  //debug
161  if ( mExportTriangulationToFile )
162  {
163  theDualEdgeTriangulation->saveAsShapefile( mTriangulationFilePath );
164  }
165 }
166 
167 int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputType type )
168 {
169  if ( !f )
170  {
171  return 1;
172  }
173 
174  const QgsGeometry* g = f->constGeometry();
175  {
176  if ( !g )
177  {
178  return 2;
179  }
180  }
181 
182  //check attribute value
183  double attributeValue = 0;
184  bool attributeConversionOk = false;
185  if ( !zCoord )
186  {
187  QVariant attributeVariant = f->attribute( attr );
188  if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
189  {
190  return 3;
191  }
192  attributeValue = attributeVariant.toDouble( &attributeConversionOk );
193  if ( !attributeConversionOk || qIsNaN( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
194  {
195  return 4;
196  }
197  }
198 
199  //parse WKB. It is ugly, but we cannot use the methods with QgsPoint because they don't contain z-values for 25D types
200  bool hasZValue = false;
201  double x, y, z;
202  QgsConstWkbPtr currentWkbPtr( g->asWkb() + 1 + sizeof( int ) );
203  //maybe a structure or break line
204  Line3D* line = nullptr;
205 
206  QGis::WkbType wkbType = g->wkbType();
207  switch ( wkbType )
208  {
209  case QGis::WKBPoint25D:
210  hasZValue = true;
211  FALLTHROUGH;
212  case QGis::WKBPoint:
213  {
214  currentWkbPtr >> x >> y;
215  if ( zCoord && hasZValue )
216  {
217  currentWkbPtr >> z;
218  }
219  else
220  {
221  z = attributeValue;
222  }
223  Point3D* thePoint = new Point3D( x, y, z );
224  if ( mTriangulation->addPoint( thePoint ) == -100 )
225  {
226  return -1;
227  }
228  break;
229  }
231  hasZValue = true;
232  FALLTHROUGH;
233  case QGis::WKBMultiPoint:
234  {
235  int nPoints;
236  currentWkbPtr >> nPoints;
237  for ( int index = 0; index < nPoints; ++index )
238  {
239  currentWkbPtr += 1 + sizeof( int );
240  currentWkbPtr >> x >> y;
241  if ( hasZValue ) //skip z-coordinate for 25D geometries
242  {
243  currentWkbPtr >> z;
244  }
245  else
246  {
247  z = attributeValue;
248  }
249  }
250  break;
251  }
253  hasZValue = true;
254  FALLTHROUGH;
255  case QGis::WKBLineString:
256  {
257  if ( type != POINTS )
258  {
259  line = new Line3D();
260  }
261  int nPoints;
262  currentWkbPtr >> nPoints;
263  for ( int index = 0; index < nPoints; ++index )
264  {
265  currentWkbPtr >> x >> y;
266  if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
267  {
268  currentWkbPtr >> z;
269  }
270  else
271  {
272  z = attributeValue;
273  }
274 
275  if ( type == POINTS )
276  {
277  //todo: handle error code -100
278  mTriangulation->addPoint( new Point3D( x, y, z ) );
279  }
280  else
281  {
282  line->insertPoint( new Point3D( x, y, z ) );
283  }
284  }
285 
286  if ( type != POINTS )
287  {
288  mTriangulation->addLine( line, type == BREAK_LINES );
289  }
290  break;
291  }
293  hasZValue = true;
294  FALLTHROUGH;
296  {
297  int nLines;
298  currentWkbPtr >> nLines;
299  for ( int index = 0; index < nLines; ++index )
300  {
301  if ( type != POINTS )
302  {
303  line = new Line3D();
304  }
305  int nPoints;
306  currentWkbPtr >> nPoints;
307  for ( int index2 = 0; index2 < nPoints; ++index2 )
308  {
309  currentWkbPtr >> x >> y;
310  if ( hasZValue ) //skip z-coordinate for 25D geometries
311  {
312  currentWkbPtr >> z;
313  }
314  else
315  {
316  z = attributeValue;
317  }
318 
319  if ( type == POINTS )
320  {
321  //todo: handle error code -100
322  mTriangulation->addPoint( new Point3D( x, y, z ) );
323  }
324  else
325  {
326  line->insertPoint( new Point3D( x, y, z ) );
327  }
328  }
329  if ( type != POINTS )
330  {
331  mTriangulation->addLine( line, type == BREAK_LINES );
332  }
333  }
334  break;
335  }
336  case QGis::WKBPolygon25D:
337  hasZValue = true;
338  FALLTHROUGH;
339  case QGis::WKBPolygon:
340  {
341  int nRings;
342  currentWkbPtr >> nRings;
343  for ( int index = 0; index < nRings; ++index )
344  {
345  if ( type != POINTS )
346  {
347  line = new Line3D();
348  }
349 
350  int nPoints;
351  currentWkbPtr >> nPoints;
352  for ( int index2 = 0; index2 < nPoints; ++index2 )
353  {
354  currentWkbPtr >> x >> y;
355  if ( hasZValue ) //skip z-coordinate for 25D geometries
356  {
357  currentWkbPtr >> z;
358  }
359  else
360  {
361  z = attributeValue;
362  }
363  if ( type == POINTS )
364  {
365  //todo: handle error code -100
366  mTriangulation->addPoint( new Point3D( x, y, z ) );
367  }
368  else
369  {
370  line->insertPoint( new Point3D( x, y, z ) );
371  }
372  }
373 
374  if ( type != POINTS )
375  {
376  mTriangulation->addLine( line, type == BREAK_LINES );
377  }
378  }
379  break;
380  }
381 
383  hasZValue = true;
384  FALLTHROUGH;
386  {
387  int nPolys;
388  currentWkbPtr >> nPolys;
389  for ( int index = 0; index < nPolys; ++index )
390  {
391  currentWkbPtr += 1 + sizeof( int );
392  int nRings;
393  currentWkbPtr >> nRings;
394  for ( int index2 = 0; index2 < nRings; ++index2 )
395  {
396  if ( type != POINTS )
397  {
398  line = new Line3D();
399  }
400  int nPoints;
401  currentWkbPtr >> nPoints;
402  for ( int index3 = 0; index3 < nPoints; ++index3 )
403  {
404  currentWkbPtr >> x >> y;
405  if ( hasZValue ) //skip z-coordinate for 25D geometries
406  {
407  currentWkbPtr >> z;
408  }
409  else
410  {
411  z = attributeValue;
412  }
413  if ( type == POINTS )
414  {
415  //todo: handle error code -100
416  mTriangulation->addPoint( new Point3D( x, y, z ) );
417  }
418  else
419  {
420  line->insertPoint( new Point3D( x, y, z ) );
421  }
422  }
423  if ( type != POINTS )
424  {
425  mTriangulation->addLine( line, type == BREAK_LINES );
426  }
427  }
428  }
429  break;
430  }
431  default:
432  //should not happen...
433  break;
434  }
435 
436  return 0;
437 }
438 
Decorator class which adds the functionality of estimating normals at the data points.
Wrapper for iterator of features from vector data provider or vector layer.
static unsigned index
QList< LayerData > mLayerData
virtual int addPoint(Point3D *p)=0
Adds a point to the triangulation Ownership is transferred to this class.
Interface class for interpolations.
QgsVectorLayer * vectorLayer
virtual void addTriangulation(Triangulation *t)
Adds an association to a triangulation.
Definition: TriDecorator.h:77
void push_back(const T &value)
virtual void setTriangulation(NormVecDecorator *tin)
LinTriangleInterpolator is a class which interpolates linearly on a triangulation.
void setWindowModality(Qt::WindowModality windowModality)
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())
Query the provider for features specified in request.
void setLabelText(const QString &text)
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:76
WkbType
Used for symbology operations.
Definition: qgis.h:57
QgsTINInterpolator(const QList< LayerData > &inputData, TIN_INTERPOLATION interpolation=Linear, bool showProgressDialog=false)
This class represents a line.
Definition: Line3D.h:24
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:187
double getZ() const
Returns the z-coordinate of the point.
Definition: Point3D.h:87
QString tr(const char *sourceText, const char *disambiguation, int n)
DualEdgeTriangulation is an implementation of a triangulation class based on the dual edge data struc...
long featureCount(QgsSymbolV2 *symbol)
Number of features rendered with specified symbol.
void setValue(int progress)
InputType
Describes the type of input data.
Point3D is a class to represent a three dimensional point.
Definition: Point3D.h:23
#define FALLTHROUGH
Definition: qgis.h:431
virtual bool saveAsShapefile(const QString &fileName) const override
Saves the triangulation as a (line) shapefile.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
QGis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
virtual void addLine(Line3D *line, bool breakline)=0
Adds a line (e.g.
void setTriangleInterpolator(TriangleInterpolator *inter) override
Sets an interpolator.
bool estimateFirstDerivatives(QProgressDialog *d=nullptr)
This method adds the functionality of estimating normals at the data points.
void insertPoint(Point3D *p)
Inserts a node behind the current position and sets the current position to this new node...
A layer together with the information about interpolation attribute / z-coordinate interpolation and ...
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:271
int interpolatePoint(double x, double y, double &result) override
Calculates interpolation value for map coordinates x, y.
This is an implementation of a Clough-Tocher interpolator based on a triangular tessellation.
const QgsGeometry * constGeometry() const
Gets a const pointer to the geometry object associated with this feature.
Definition: qgsfeature.cpp:82
virtual bool calcPoint(double x, double y, Point3D *result)=0
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point...
bool isValid() const
double toDouble(bool *ok) const
bool nextFeature(QgsFeature &f)
const unsigned char * asWkb() const
Returns the buffer containing this geometry in WKB format.