Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
27 #include "qgsvectorlayer.h"
28 #include <QProgressDialog>
29 
30 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData>& inputData, TIN_INTERPOLATION interpolation, bool showProgressDialog )
31  : QgsInterpolator( inputData )
32  , mTriangulation( 0 )
33  , mTriangleInterpolator( 0 )
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 
69 {
70  DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( 100000, 0 );
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  QList<LayerData>::iterator layerDataIt = mLayerData.begin();
88  for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
89  {
90  if ( layerDataIt->vectorLayer )
91  {
92  nFeatures += layerDataIt->vectorLayer->featureCount();
93  }
94  }
95  }
96 
97  QProgressDialog* theProgressDialog = 0;
98  if ( mShowProgressDialog )
99  {
100  theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, 0 );
101  theProgressDialog->setWindowModality( Qt::WindowModal );
102  }
103 
104 
105  QgsFeature f;
106  QList<LayerData>::iterator layerDataIt = mLayerData.begin();
107  for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
108  {
109  if ( layerDataIt->vectorLayer )
110  {
111  QgsAttributeList attList;
112  if ( !layerDataIt->zCoordInterpolation )
113  {
114  attList.push_back( layerDataIt->interpolationAttribute );
115  }
116  layerDataIt->vectorLayer->select( attList );
117  while ( layerDataIt->vectorLayer->nextFeature( f ) )
118  {
119  if ( mShowProgressDialog )
120  {
121  if ( theProgressDialog->wasCanceled() )
122  {
123  break;
124  }
125  theProgressDialog->setValue( nProcessedFeatures );
126  }
127  insertData( &f, layerDataIt->zCoordInterpolation, layerDataIt->interpolationAttribute, layerDataIt->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 = 0;
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
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  QgsGeometry* g = f->geometry();
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  QgsAttributeMap attMap = f->attributeMap();
188  QgsAttributeMap::const_iterator att_it = attMap.find( attr );
189  if ( att_it == attMap.end() ) //attribute not found, something must be wrong (e.g. NULL value)
190  {
191  return 3;
192  }
193  attributeValue = att_it.value().toDouble( &attributeConversionOk );
194  if ( !attributeConversionOk || qIsNaN( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
195  {
196  return 4;
197  }
198  }
199 
200  //parse WKB. It is ugly, but we cannot use the methods with QgsPoint because they don't contain z-values for 25D types
201  bool hasZValue = false;
202  double x, y, z;
203  unsigned char* currentWkbPtr = g->asWkb();
204  //maybe a structure or break line
205  Line3D* line = 0;
206 
207  QGis::WkbType wkbType = g->wkbType();
208  switch ( wkbType )
209  {
210  case QGis::WKBPoint25D:
211  hasZValue = true;
212  case QGis::WKBPoint:
213  {
214  currentWkbPtr += ( 1 + sizeof( int ) );
215  x = *(( double * )( currentWkbPtr ) );
216  currentWkbPtr += sizeof( double );
217  y = *(( double * )( currentWkbPtr ) );
218  if ( zCoord && hasZValue )
219  {
220  currentWkbPtr += sizeof( double );
221  z = *(( double * )( currentWkbPtr ) );
222  }
223  else
224  {
225  z = attributeValue;
226  }
227  Point3D* thePoint = new Point3D( x, y, z );
228  if ( mTriangulation->addPoint( thePoint ) == -100 )
229  {
230  return -1;
231  }
232  break;
233  }
235  hasZValue = true;
236  case QGis::WKBMultiPoint:
237  {
238  currentWkbPtr += ( 1 + sizeof( int ) );
239  int* npoints = ( int* )currentWkbPtr;
240  currentWkbPtr += sizeof( int );
241  for ( int index = 0; index < *npoints; ++index )
242  {
243  currentWkbPtr += ( 1 + sizeof( int ) ); //skip endian and point type
244  x = *(( double* )currentWkbPtr );
245  currentWkbPtr += sizeof( double );
246  y = *(( double* )currentWkbPtr );
247  currentWkbPtr += sizeof( double );
248  if ( hasZValue ) //skip z-coordinate for 25D geometries
249  {
250  z = *(( double* )currentWkbPtr );
251  currentWkbPtr += sizeof( double );
252  }
253  else
254  {
255  z = attributeValue;
256  }
257  }
258  break;
259  }
261  hasZValue = true;
262  case QGis::WKBLineString:
263  {
264  if ( type != POINTS )
265  {
266  line = new Line3D();
267  }
268  currentWkbPtr += ( 1 + sizeof( int ) );
269  int* npoints = ( int* )currentWkbPtr;
270  currentWkbPtr += sizeof( int );
271  for ( int index = 0; index < *npoints; ++index )
272  {
273  x = *(( double * )( currentWkbPtr ) );
274  currentWkbPtr += sizeof( double );
275  y = *(( double * )( currentWkbPtr ) );
276  currentWkbPtr += sizeof( double );
277  if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
278  {
279  z = *(( double * )( currentWkbPtr ) );
280  }
281  else
282  {
283  z = attributeValue;
284  }
285  if ( hasZValue )
286  {
287  currentWkbPtr += sizeof( double );
288  }
289 
290  if ( type == POINTS )
291  {
292  //todo: handle error code -100
293  mTriangulation->addPoint( new Point3D( x, y, z ) );
294  }
295  else
296  {
297  line->insertPoint( new Point3D( x, y, z ) );
298  }
299  }
300 
301  if ( type != POINTS )
302  {
303  mTriangulation->addLine( line, type == BREAK_LINES );
304  }
305  break;
306  }
308  hasZValue = true;
310  {
311  currentWkbPtr += ( 1 + sizeof( int ) );
312  int* nlines = ( int* )currentWkbPtr;
313  int* npoints = 0;
314  currentWkbPtr += sizeof( int );
315  for ( int index = 0; index < *nlines; ++index )
316  {
317  if ( type != POINTS )
318  {
319  line = new Line3D();
320  }
321  currentWkbPtr += ( sizeof( int ) + 1 );
322  npoints = ( int* )currentWkbPtr;
323  currentWkbPtr += sizeof( int );
324  for ( int index2 = 0; index2 < *npoints; ++index2 )
325  {
326  x = *(( double* )currentWkbPtr );
327  currentWkbPtr += sizeof( double );
328  y = *(( double* )currentWkbPtr );
329  currentWkbPtr += sizeof( double );
330 
331  if ( hasZValue ) //skip z-coordinate for 25D geometries
332  {
333  z = *(( double* ) currentWkbPtr );
334  currentWkbPtr += sizeof( double );
335  }
336  else
337  {
338  z = attributeValue;
339  }
340 
341  if ( type == POINTS )
342  {
343  //todo: handle error code -100
344  mTriangulation->addPoint( new Point3D( x, y, z ) );
345  }
346  else
347  {
348  line->insertPoint( new Point3D( x, y, z ) );
349  }
350  }
351  if ( type != POINTS )
352  {
353  mTriangulation->addLine( line, type == BREAK_LINES );
354  }
355  }
356  break;
357  }
358  case QGis::WKBPolygon25D:
359  hasZValue = true;
360  case QGis::WKBPolygon:
361  {
362  currentWkbPtr += ( 1 + sizeof( int ) );
363  int* nrings = ( int* )currentWkbPtr;
364  currentWkbPtr += sizeof( int );
365  int* npoints;
366  for ( int index = 0; index < *nrings; ++index )
367  {
368  if ( type != POINTS )
369  {
370  line = new Line3D();
371  }
372 
373  npoints = ( int* )currentWkbPtr;
374  currentWkbPtr += sizeof( int );
375  for ( int index2 = 0; index2 < *npoints; ++index2 )
376  {
377  x = *(( double* )currentWkbPtr );
378  currentWkbPtr += sizeof( double );
379  y = *(( double* )currentWkbPtr );
380  currentWkbPtr += sizeof( double );
381  if ( hasZValue ) //skip z-coordinate for 25D geometries
382  {
383  z = *(( double* )currentWkbPtr );;
384  currentWkbPtr += sizeof( double );
385  }
386  else
387  {
388  z = attributeValue;
389  }
390  if ( type == POINTS )
391  {
392  //todo: handle error code -100
393  mTriangulation->addPoint( new Point3D( x, y, z ) );
394  }
395  else
396  {
397  line->insertPoint( new Point3D( x, y, z ) );
398  }
399  }
400 
401  if ( type != POINTS )
402  {
403  mTriangulation->addLine( line, type == BREAK_LINES );
404  }
405  }
406  break;
407  }
408 
410  hasZValue = true;
412  {
413  currentWkbPtr += ( 1 + sizeof( int ) );
414  int* npolys = ( int* )currentWkbPtr;
415  int* nrings;
416  int* npoints;
417  currentWkbPtr += sizeof( int );
418  for ( int index = 0; index < *npolys; ++index )
419  {
420  currentWkbPtr += ( 1 + sizeof( int ) ); //skip endian and polygon type
421  nrings = ( int* )currentWkbPtr;
422  currentWkbPtr += sizeof( int );
423  for ( int index2 = 0; index2 < *nrings; ++index2 )
424  {
425  if ( type != POINTS )
426  {
427  line = new Line3D();
428  }
429  npoints = ( int* )currentWkbPtr;
430  currentWkbPtr += sizeof( int );
431  for ( int index3 = 0; index3 < *npoints; ++index3 )
432  {
433  x = *(( double* )currentWkbPtr );
434  currentWkbPtr += sizeof( double );
435  y = *(( double* )currentWkbPtr );
436  currentWkbPtr += sizeof( double );
437  if ( hasZValue ) //skip z-coordinate for 25D geometries
438  {
439  z = *(( double* )currentWkbPtr );
440  currentWkbPtr += sizeof( double );
441  }
442  else
443  {
444  z = attributeValue;
445  }
446  if ( type == POINTS )
447  {
448  //todo: handle error code -100
449  mTriangulation->addPoint( new Point3D( x, y, z ) );
450  }
451  else
452  {
453  line->insertPoint( new Point3D( x, y, z ) );
454  }
455  }
456  if ( type != POINTS )
457  {
458  mTriangulation->addLine( line, type == BREAK_LINES );
459  }
460  }
461  }
462  break;
463  }
464  default:
465  //should not happen...
466  break;
467  }
468 
469  return 0;
470 }
471