Quantum GIS API Documentation
1.8
|
00001 /*************************************************************************** 00002 qgstininterpolator.cpp 00003 ---------------------- 00004 begin : March 10, 2008 00005 copyright : (C) 2008 by Marco Hugentobler 00006 email : marco dot hugentobler at karto dot baug dot ethz dot ch 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * * 00011 * This program is free software; you can redistribute it and/or modify * 00012 * it under the terms of the GNU General Public License as published by * 00013 * the Free Software Foundation; either version 2 of the License, or * 00014 * (at your option) any later version. * 00015 * * 00016 ***************************************************************************/ 00017 00018 #include "qgstininterpolator.h" 00019 #include "CloughTocherInterpolator.h" 00020 #include "DualEdgeTriangulation.h" 00021 #include "NormVecDecorator.h" 00022 #include "LinTriangleInterpolator.h" 00023 #include "Point3D.h" 00024 #include "qgsfeature.h" 00025 #include "qgsgeometry.h" 00026 #include "qgssinglesymbolrenderer.h" 00027 #include "qgsvectorlayer.h" 00028 #include <QProgressDialog> 00029 00030 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData>& inputData, TIN_INTERPOLATION interpolation, bool showProgressDialog ) 00031 : QgsInterpolator( inputData ) 00032 , mTriangulation( 0 ) 00033 , mTriangleInterpolator( 0 ) 00034 , mIsInitialized( false ) 00035 , mShowProgressDialog( showProgressDialog ) 00036 , mExportTriangulationToFile( false ) 00037 , mInterpolation( interpolation ) 00038 { 00039 } 00040 00041 QgsTINInterpolator::~QgsTINInterpolator() 00042 { 00043 delete mTriangulation; 00044 delete mTriangleInterpolator; 00045 } 00046 00047 int QgsTINInterpolator::interpolatePoint( double x, double y, double& result ) 00048 { 00049 if ( !mIsInitialized ) 00050 { 00051 initialize(); 00052 } 00053 00054 if ( !mTriangleInterpolator ) 00055 { 00056 return 1; 00057 } 00058 00059 Point3D r; 00060 if ( !mTriangleInterpolator->calcPoint( x, y, &r ) ) 00061 { 00062 return 2; 00063 } 00064 result = r.getZ(); 00065 return 0; 00066 } 00067 00068 void QgsTINInterpolator::initialize() 00069 { 00070 DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( 100000, 0 ); 00071 if ( mInterpolation == CloughTocher ) 00072 { 00073 NormVecDecorator* dec = new NormVecDecorator(); 00074 dec->addTriangulation( theDualEdgeTriangulation ); 00075 mTriangulation = dec; 00076 } 00077 else 00078 { 00079 mTriangulation = theDualEdgeTriangulation; 00080 } 00081 00082 //get number of features if we use a progress bar 00083 int nFeatures = 0; 00084 int nProcessedFeatures = 0; 00085 if ( mShowProgressDialog ) 00086 { 00087 QList<LayerData>::iterator layerDataIt = mLayerData.begin(); 00088 for ( ; layerDataIt != mLayerData.end(); ++layerDataIt ) 00089 { 00090 if ( layerDataIt->vectorLayer ) 00091 { 00092 nFeatures += layerDataIt->vectorLayer->featureCount(); 00093 } 00094 } 00095 } 00096 00097 QProgressDialog* theProgressDialog = 0; 00098 if ( mShowProgressDialog ) 00099 { 00100 theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, 0 ); 00101 theProgressDialog->setWindowModality( Qt::WindowModal ); 00102 } 00103 00104 00105 QgsFeature f; 00106 QList<LayerData>::iterator layerDataIt = mLayerData.begin(); 00107 for ( ; layerDataIt != mLayerData.end(); ++layerDataIt ) 00108 { 00109 if ( layerDataIt->vectorLayer ) 00110 { 00111 QgsAttributeList attList; 00112 if ( !layerDataIt->zCoordInterpolation ) 00113 { 00114 attList.push_back( layerDataIt->interpolationAttribute ); 00115 } 00116 layerDataIt->vectorLayer->select( attList ); 00117 while ( layerDataIt->vectorLayer->nextFeature( f ) ) 00118 { 00119 if ( mShowProgressDialog ) 00120 { 00121 if ( theProgressDialog->wasCanceled() ) 00122 { 00123 break; 00124 } 00125 theProgressDialog->setValue( nProcessedFeatures ); 00126 } 00127 insertData( &f, layerDataIt->zCoordInterpolation, layerDataIt->interpolationAttribute, layerDataIt->mInputType ); 00128 ++nProcessedFeatures; 00129 } 00130 } 00131 } 00132 00133 delete theProgressDialog; 00134 00135 if ( mInterpolation == CloughTocher ) 00136 { 00137 CloughTocherInterpolator* ctInterpolator = new CloughTocherInterpolator(); 00138 NormVecDecorator* dec = dynamic_cast<NormVecDecorator*>( mTriangulation ); 00139 if ( dec ) 00140 { 00141 QProgressDialog* progressDialog = 0; 00142 if ( mShowProgressDialog ) //show a progress dialog because it can take a long time... 00143 { 00144 progressDialog = new QProgressDialog(); 00145 progressDialog->setLabelText( QObject::tr( "Estimating normal derivatives..." ) ); 00146 } 00147 dec->estimateFirstDerivatives( progressDialog ); 00148 delete progressDialog; 00149 ctInterpolator->setTriangulation( dec ); 00150 dec->setTriangleInterpolator( ctInterpolator ); 00151 mTriangleInterpolator = ctInterpolator; 00152 } 00153 } 00154 else //linear 00155 { 00156 mTriangleInterpolator = new LinTriangleInterpolator( theDualEdgeTriangulation ); 00157 } 00158 mIsInitialized = true; 00159 00160 //debug 00161 if ( mExportTriangulationToFile ) 00162 { 00163 theDualEdgeTriangulation->saveAsShapefile( mTriangulationFilePath ); 00164 } 00165 } 00166 00167 int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputType type ) 00168 { 00169 if ( !f ) 00170 { 00171 return 1; 00172 } 00173 00174 QgsGeometry* g = f->geometry(); 00175 { 00176 if ( !g ) 00177 { 00178 return 2; 00179 } 00180 } 00181 00182 //check attribute value 00183 double attributeValue = 0; 00184 bool attributeConversionOk = false; 00185 if ( !zCoord ) 00186 { 00187 QgsAttributeMap attMap = f->attributeMap(); 00188 QgsAttributeMap::const_iterator att_it = attMap.find( attr ); 00189 if ( att_it == attMap.end() ) //attribute not found, something must be wrong (e.g. NULL value) 00190 { 00191 return 3; 00192 } 00193 attributeValue = att_it.value().toDouble( &attributeConversionOk ); 00194 if ( !attributeConversionOk || qIsNaN( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation 00195 { 00196 return 4; 00197 } 00198 } 00199 00200 //parse WKB. It is ugly, but we cannot use the methods with QgsPoint because they don't contain z-values for 25D types 00201 bool hasZValue = false; 00202 double x, y, z; 00203 unsigned char* currentWkbPtr = g->asWkb(); 00204 //maybe a structure or break line 00205 Line3D* line = 0; 00206 00207 QGis::WkbType wkbType = g->wkbType(); 00208 switch ( wkbType ) 00209 { 00210 case QGis::WKBPoint25D: 00211 hasZValue = true; 00212 case QGis::WKBPoint: 00213 { 00214 currentWkbPtr += ( 1 + sizeof( int ) ); 00215 x = *(( double * )( currentWkbPtr ) ); 00216 currentWkbPtr += sizeof( double ); 00217 y = *(( double * )( currentWkbPtr ) ); 00218 if ( zCoord && hasZValue ) 00219 { 00220 currentWkbPtr += sizeof( double ); 00221 z = *(( double * )( currentWkbPtr ) ); 00222 } 00223 else 00224 { 00225 z = attributeValue; 00226 } 00227 Point3D* thePoint = new Point3D( x, y, z ); 00228 if ( mTriangulation->addPoint( thePoint ) == -100 ) 00229 { 00230 return -1; 00231 } 00232 break; 00233 } 00234 case QGis::WKBMultiPoint25D: 00235 hasZValue = true; 00236 case QGis::WKBMultiPoint: 00237 { 00238 currentWkbPtr += ( 1 + sizeof( int ) ); 00239 int* npoints = ( int* )currentWkbPtr; 00240 currentWkbPtr += sizeof( int ); 00241 for ( int index = 0; index < *npoints; ++index ) 00242 { 00243 currentWkbPtr += ( 1 + sizeof( int ) ); //skip endian and point type 00244 x = *(( double* )currentWkbPtr ); 00245 currentWkbPtr += sizeof( double ); 00246 y = *(( double* )currentWkbPtr ); 00247 currentWkbPtr += sizeof( double ); 00248 if ( hasZValue ) //skip z-coordinate for 25D geometries 00249 { 00250 z = *(( double* )currentWkbPtr ); 00251 currentWkbPtr += sizeof( double ); 00252 } 00253 else 00254 { 00255 z = attributeValue; 00256 } 00257 } 00258 break; 00259 } 00260 case QGis::WKBLineString25D: 00261 hasZValue = true; 00262 case QGis::WKBLineString: 00263 { 00264 if ( type != POINTS ) 00265 { 00266 line = new Line3D(); 00267 } 00268 currentWkbPtr += ( 1 + sizeof( int ) ); 00269 int* npoints = ( int* )currentWkbPtr; 00270 currentWkbPtr += sizeof( int ); 00271 for ( int index = 0; index < *npoints; ++index ) 00272 { 00273 x = *(( double * )( currentWkbPtr ) ); 00274 currentWkbPtr += sizeof( double ); 00275 y = *(( double * )( currentWkbPtr ) ); 00276 currentWkbPtr += sizeof( double ); 00277 if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries 00278 { 00279 z = *(( double * )( currentWkbPtr ) ); 00280 } 00281 else 00282 { 00283 z = attributeValue; 00284 } 00285 if ( hasZValue ) 00286 { 00287 currentWkbPtr += sizeof( double ); 00288 } 00289 00290 if ( type == POINTS ) 00291 { 00292 //todo: handle error code -100 00293 mTriangulation->addPoint( new Point3D( x, y, z ) ); 00294 } 00295 else 00296 { 00297 line->insertPoint( new Point3D( x, y, z ) ); 00298 } 00299 } 00300 00301 if ( type != POINTS ) 00302 { 00303 mTriangulation->addLine( line, type == BREAK_LINES ); 00304 } 00305 break; 00306 } 00307 case QGis::WKBMultiLineString25D: 00308 hasZValue = true; 00309 case QGis::WKBMultiLineString: 00310 { 00311 currentWkbPtr += ( 1 + sizeof( int ) ); 00312 int* nlines = ( int* )currentWkbPtr; 00313 int* npoints = 0; 00314 currentWkbPtr += sizeof( int ); 00315 for ( int index = 0; index < *nlines; ++index ) 00316 { 00317 if ( type != POINTS ) 00318 { 00319 line = new Line3D(); 00320 } 00321 currentWkbPtr += ( sizeof( int ) + 1 ); 00322 npoints = ( int* )currentWkbPtr; 00323 currentWkbPtr += sizeof( int ); 00324 for ( int index2 = 0; index2 < *npoints; ++index2 ) 00325 { 00326 x = *(( double* )currentWkbPtr ); 00327 currentWkbPtr += sizeof( double ); 00328 y = *(( double* )currentWkbPtr ); 00329 currentWkbPtr += sizeof( double ); 00330 00331 if ( hasZValue ) //skip z-coordinate for 25D geometries 00332 { 00333 z = *(( double* ) currentWkbPtr ); 00334 currentWkbPtr += sizeof( double ); 00335 } 00336 else 00337 { 00338 z = attributeValue; 00339 } 00340 00341 if ( type == POINTS ) 00342 { 00343 //todo: handle error code -100 00344 mTriangulation->addPoint( new Point3D( x, y, z ) ); 00345 } 00346 else 00347 { 00348 line->insertPoint( new Point3D( x, y, z ) ); 00349 } 00350 } 00351 if ( type != POINTS ) 00352 { 00353 mTriangulation->addLine( line, type == BREAK_LINES ); 00354 } 00355 } 00356 break; 00357 } 00358 case QGis::WKBPolygon25D: 00359 hasZValue = true; 00360 case QGis::WKBPolygon: 00361 { 00362 currentWkbPtr += ( 1 + sizeof( int ) ); 00363 int* nrings = ( int* )currentWkbPtr; 00364 currentWkbPtr += sizeof( int ); 00365 int* npoints; 00366 for ( int index = 0; index < *nrings; ++index ) 00367 { 00368 if ( type != POINTS ) 00369 { 00370 line = new Line3D(); 00371 } 00372 00373 npoints = ( int* )currentWkbPtr; 00374 currentWkbPtr += sizeof( int ); 00375 for ( int index2 = 0; index2 < *npoints; ++index2 ) 00376 { 00377 x = *(( double* )currentWkbPtr ); 00378 currentWkbPtr += sizeof( double ); 00379 y = *(( double* )currentWkbPtr ); 00380 currentWkbPtr += sizeof( double ); 00381 if ( hasZValue ) //skip z-coordinate for 25D geometries 00382 { 00383 z = *(( double* )currentWkbPtr );; 00384 currentWkbPtr += sizeof( double ); 00385 } 00386 else 00387 { 00388 z = attributeValue; 00389 } 00390 if ( type == POINTS ) 00391 { 00392 //todo: handle error code -100 00393 mTriangulation->addPoint( new Point3D( x, y, z ) ); 00394 } 00395 else 00396 { 00397 line->insertPoint( new Point3D( x, y, z ) ); 00398 } 00399 } 00400 00401 if ( type != POINTS ) 00402 { 00403 mTriangulation->addLine( line, type == BREAK_LINES ); 00404 } 00405 } 00406 break; 00407 } 00408 00409 case QGis::WKBMultiPolygon25D: 00410 hasZValue = true; 00411 case QGis::WKBMultiPolygon: 00412 { 00413 currentWkbPtr += ( 1 + sizeof( int ) ); 00414 int* npolys = ( int* )currentWkbPtr; 00415 int* nrings; 00416 int* npoints; 00417 currentWkbPtr += sizeof( int ); 00418 for ( int index = 0; index < *npolys; ++index ) 00419 { 00420 currentWkbPtr += ( 1 + sizeof( int ) ); //skip endian and polygon type 00421 nrings = ( int* )currentWkbPtr; 00422 currentWkbPtr += sizeof( int ); 00423 for ( int index2 = 0; index2 < *nrings; ++index2 ) 00424 { 00425 if ( type != POINTS ) 00426 { 00427 line = new Line3D(); 00428 } 00429 npoints = ( int* )currentWkbPtr; 00430 currentWkbPtr += sizeof( int ); 00431 for ( int index3 = 0; index3 < *npoints; ++index3 ) 00432 { 00433 x = *(( double* )currentWkbPtr ); 00434 currentWkbPtr += sizeof( double ); 00435 y = *(( double* )currentWkbPtr ); 00436 currentWkbPtr += sizeof( double ); 00437 if ( hasZValue ) //skip z-coordinate for 25D geometries 00438 { 00439 z = *(( double* )currentWkbPtr ); 00440 currentWkbPtr += sizeof( double ); 00441 } 00442 else 00443 { 00444 z = attributeValue; 00445 } 00446 if ( type == POINTS ) 00447 { 00448 //todo: handle error code -100 00449 mTriangulation->addPoint( new Point3D( x, y, z ) ); 00450 } 00451 else 00452 { 00453 line->insertPoint( new Point3D( x, y, z ) ); 00454 } 00455 } 00456 if ( type != POINTS ) 00457 { 00458 mTriangulation->addLine( line, type == BREAK_LINES ); 00459 } 00460 } 00461 } 00462 break; 00463 } 00464 default: 00465 //should not happen... 00466 break; 00467 } 00468 00469 return 0; 00470 } 00471