00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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 #ifndef Q_OS_MACX
00031 #include <cmath>
00032 #else
00033 #include <math.h>
00034 #endif
00035 #ifdef WIN32
00036 #include <float.h>
00037 #define isnan(f) _isnan(f)
00038 #endif
00039
00040 QgsTINInterpolator::QgsTINInterpolator( const QList<LayerData>& inputData, TIN_INTERPOLATION interpolation, bool showProgressDialog ): QgsInterpolator( inputData ), mTriangulation( 0 ), \
00041 mTriangleInterpolator( 0 ), mIsInitialized( false ), mShowProgressDialog( showProgressDialog ), mExportTriangulationToFile( false ), mInterpolation( interpolation )
00042 {
00043 }
00044
00045 QgsTINInterpolator::~QgsTINInterpolator()
00046 {
00047 delete mTriangulation;
00048 delete mTriangleInterpolator;
00049 }
00050
00051 int QgsTINInterpolator::interpolatePoint( double x, double y, double& result )
00052 {
00053 if ( !mIsInitialized )
00054 {
00055 initialize();
00056 }
00057
00058 if ( !mTriangleInterpolator )
00059 {
00060 return 1;
00061 }
00062
00063 Point3D r;
00064 if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
00065 {
00066 return 2;
00067 }
00068 result = r.getZ();
00069 return 0;
00070 }
00071
00072 void QgsTINInterpolator::initialize()
00073 {
00074 DualEdgeTriangulation* theDualEdgeTriangulation = new DualEdgeTriangulation( 100000, 0 );
00075 if ( mInterpolation == CloughTocher )
00076 {
00077 NormVecDecorator* dec = new NormVecDecorator();
00078 dec->addTriangulation( theDualEdgeTriangulation );
00079 mTriangulation = dec;
00080 }
00081 else
00082 {
00083 mTriangulation = theDualEdgeTriangulation;
00084 }
00085
00086
00087 int nFeatures = 0;
00088 int nProcessedFeatures = 0;
00089 if ( mShowProgressDialog )
00090 {
00091 QList<LayerData>::iterator layerDataIt = mLayerData.begin();
00092 for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
00093 {
00094 if ( layerDataIt->vectorLayer )
00095 {
00096 nFeatures += layerDataIt->vectorLayer->featureCount();
00097 }
00098 }
00099 }
00100
00101 QProgressDialog* theProgressDialog = 0;
00102 if ( mShowProgressDialog )
00103 {
00104 theProgressDialog = new QProgressDialog( QObject::tr( "Building triangulation..." ), QObject::tr( "Abort" ), 0, nFeatures, 0 );
00105 theProgressDialog->setWindowModality( Qt::WindowModal );
00106 }
00107
00108
00109 QgsFeature f;
00110 QList<LayerData>::iterator layerDataIt = mLayerData.begin();
00111 for ( ; layerDataIt != mLayerData.end(); ++layerDataIt )
00112 {
00113 if ( layerDataIt->vectorLayer )
00114 {
00115 QgsAttributeList attList;
00116 if ( !layerDataIt->zCoordInterpolation )
00117 {
00118 attList.push_back( layerDataIt->interpolationAttribute );
00119 }
00120 layerDataIt->vectorLayer->select( attList );
00121 while ( layerDataIt->vectorLayer->nextFeature( f ) )
00122 {
00123 if ( mShowProgressDialog )
00124 {
00125 if ( theProgressDialog->wasCanceled() )
00126 {
00127 break;
00128 }
00129 theProgressDialog->setValue( nProcessedFeatures );
00130 }
00131 insertData( &f, layerDataIt->zCoordInterpolation, layerDataIt->interpolationAttribute, layerDataIt->mInputType );
00132 ++nProcessedFeatures;
00133 }
00134 }
00135 }
00136
00137 delete theProgressDialog;
00138
00139 if ( mInterpolation == CloughTocher )
00140 {
00141 CloughTocherInterpolator* ctInterpolator = new CloughTocherInterpolator();
00142 NormVecDecorator* dec = dynamic_cast<NormVecDecorator*>( mTriangulation );
00143 if ( dec )
00144 {
00145 QProgressDialog* progressDialog = 0;
00146 if ( mShowProgressDialog )
00147 {
00148 progressDialog = new QProgressDialog();
00149 progressDialog->setLabelText( QObject::tr( "Estimating normal derivatives..." ) );
00150 }
00151 dec->estimateFirstDerivatives( progressDialog );
00152 delete progressDialog;
00153 ctInterpolator->setTriangulation( dec );
00154 dec->setTriangleInterpolator( ctInterpolator );
00155 mTriangleInterpolator = ctInterpolator;
00156 }
00157 }
00158 else
00159 {
00160 mTriangleInterpolator = new LinTriangleInterpolator( theDualEdgeTriangulation );
00161 }
00162 mIsInitialized = true;
00163
00164
00165 if ( mExportTriangulationToFile )
00166 {
00167 theDualEdgeTriangulation->saveAsShapefile( mTriangulationFilePath );
00168 }
00169 }
00170
00171 int QgsTINInterpolator::insertData( QgsFeature* f, bool zCoord, int attr, InputType type )
00172 {
00173 if ( !f )
00174 {
00175 return 1;
00176 }
00177
00178 QgsGeometry* g = f->geometry();
00179 {
00180 if ( !g )
00181 {
00182 return 2;
00183 }
00184 }
00185
00186
00187 double attributeValue = 0;
00188 bool attributeConversionOk = false;
00189 if ( !zCoord )
00190 {
00191 QgsAttributeMap attMap = f->attributeMap();
00192 QgsAttributeMap::const_iterator att_it = attMap.find( attr );
00193 if ( att_it == attMap.end() )
00194 {
00195 return 3;
00196 }
00197 attributeValue = att_it.value().toDouble( &attributeConversionOk );
00198 if ( !attributeConversionOk || isnan( attributeValue ) )
00199 {
00200 return 4;
00201 }
00202 }
00203
00204
00205 bool hasZValue = false;
00206 double x, y, z;
00207 unsigned char* currentWkbPtr = g->asWkb();
00208
00209 Line3D* line = 0;
00210
00211 QGis::WkbType wkbType = g->wkbType();
00212 switch ( wkbType )
00213 {
00214 case QGis::WKBPoint25D:
00215 hasZValue = true;
00216 case QGis::WKBPoint:
00217 {
00218 currentWkbPtr += ( 1 + sizeof( int ) );
00219 x = *(( double * )( currentWkbPtr ) );
00220 currentWkbPtr += sizeof( double );
00221 y = *(( double * )( currentWkbPtr ) );
00222 if ( zCoord && hasZValue )
00223 {
00224 currentWkbPtr += sizeof( double );
00225 z = *(( double * )( currentWkbPtr ) );
00226 }
00227 else
00228 {
00229 z = attributeValue;
00230 }
00231 Point3D* thePoint = new Point3D( x, y, z );
00232 if ( mTriangulation->addPoint( thePoint ) == -100 )
00233 {
00234 return -1;
00235 }
00236 break;
00237 }
00238 case QGis::WKBMultiPoint25D:
00239 hasZValue = true;
00240 case QGis::WKBMultiPoint:
00241 {
00242 currentWkbPtr += ( 1 + sizeof( int ) );
00243 int* npoints = ( int* )currentWkbPtr;
00244 currentWkbPtr += sizeof( int );
00245 for ( int index = 0; index < *npoints; ++index )
00246 {
00247 currentWkbPtr += ( 1 + sizeof( int ) );
00248 x = *(( double* )currentWkbPtr );
00249 currentWkbPtr += sizeof( double );
00250 y = *(( double* )currentWkbPtr );
00251 currentWkbPtr += sizeof( double );
00252 if ( hasZValue )
00253 {
00254 z = *(( double* )currentWkbPtr );
00255 currentWkbPtr += sizeof( double );
00256 }
00257 else
00258 {
00259 z = attributeValue;
00260 }
00261 }
00262 break;
00263 }
00264 case QGis::WKBLineString25D:
00265 hasZValue = true;
00266 case QGis::WKBLineString:
00267 {
00268 if ( type != POINTS )
00269 {
00270 line = new Line3D();
00271 }
00272 currentWkbPtr += ( 1 + sizeof( int ) );
00273 int* npoints = ( int* )currentWkbPtr;
00274 currentWkbPtr += sizeof( int );
00275 for ( int index = 0; index < *npoints; ++index )
00276 {
00277 x = *(( double * )( currentWkbPtr ) );
00278 currentWkbPtr += sizeof( double );
00279 y = *(( double * )( currentWkbPtr ) );
00280 currentWkbPtr += sizeof( double );
00281 if ( zCoord && hasZValue )
00282 {
00283 z = *(( double * )( currentWkbPtr ) );
00284 }
00285 else
00286 {
00287 z = attributeValue;
00288 }
00289 if ( hasZValue )
00290 {
00291 currentWkbPtr += sizeof( double );
00292 }
00293
00294 if ( type == POINTS )
00295 {
00296
00297 mTriangulation->addPoint( new Point3D( x, y, z ) );
00298 }
00299 else
00300 {
00301 line->insertPoint( new Point3D( x, y, z ) );
00302 }
00303 }
00304
00305 if ( type != POINTS )
00306 {
00307 mTriangulation->addLine( line, type == BREAK_LINES );
00308 }
00309 break;
00310 }
00311 case QGis::WKBMultiLineString25D:
00312 hasZValue = true;
00313 case QGis::WKBMultiLineString:
00314 {
00315 currentWkbPtr += ( 1 + sizeof( int ) );
00316 int* nlines = ( int* )currentWkbPtr;
00317 int* npoints = 0;
00318 currentWkbPtr += sizeof( int );
00319 for ( int index = 0; index < *nlines; ++index )
00320 {
00321 if ( type != POINTS )
00322 {
00323 line = new Line3D();
00324 }
00325 currentWkbPtr += ( sizeof( int ) + 1 );
00326 npoints = ( int* )currentWkbPtr;
00327 currentWkbPtr += sizeof( int );
00328 for ( int index2 = 0; index2 < *npoints; ++index2 )
00329 {
00330 x = *(( double* )currentWkbPtr );
00331 currentWkbPtr += sizeof( double );
00332 y = *(( double* )currentWkbPtr );
00333 currentWkbPtr += sizeof( double );
00334
00335 if ( hasZValue )
00336 {
00337 z = *(( double* ) currentWkbPtr );
00338 currentWkbPtr += sizeof( double );
00339 }
00340 else
00341 {
00342 z = attributeValue;
00343 }
00344
00345 if ( type == POINTS )
00346 {
00347
00348 mTriangulation->addPoint( new Point3D( x, y, z ) );
00349 }
00350 else
00351 {
00352 line->insertPoint( new Point3D( x, y, z ) );
00353 }
00354 }
00355 if ( type != POINTS )
00356 {
00357 mTriangulation->addLine( line, type == BREAK_LINES );
00358 }
00359 }
00360 break;
00361 }
00362 case QGis::WKBPolygon25D:
00363 hasZValue = true;
00364 case QGis::WKBPolygon:
00365 {
00366 currentWkbPtr += ( 1 + sizeof( int ) );
00367 int* nrings = ( int* )currentWkbPtr;
00368 currentWkbPtr += sizeof( int );
00369 int* npoints;
00370 for ( int index = 0; index < *nrings; ++index )
00371 {
00372 if ( type != POINTS )
00373 {
00374 line = new Line3D();
00375 }
00376
00377 npoints = ( int* )currentWkbPtr;
00378 currentWkbPtr += sizeof( int );
00379 for ( int index2 = 0; index2 < *npoints; ++index2 )
00380 {
00381 x = *(( double* )currentWkbPtr );
00382 currentWkbPtr += sizeof( double );
00383 y = *(( double* )currentWkbPtr );
00384 currentWkbPtr += sizeof( double );
00385 if ( hasZValue )
00386 {
00387 z = *(( double* )currentWkbPtr );;
00388 currentWkbPtr += sizeof( double );
00389 }
00390 else
00391 {
00392 z = attributeValue;
00393 }
00394 if ( type == POINTS )
00395 {
00396
00397 mTriangulation->addPoint( new Point3D( x, y, z ) );
00398 }
00399 else
00400 {
00401 line->insertPoint( new Point3D( x, y, z ) );
00402 }
00403 }
00404
00405 if ( type != POINTS )
00406 {
00407 mTriangulation->addLine( line, type == BREAK_LINES );
00408 }
00409 }
00410 break;
00411 }
00412
00413 case QGis::WKBMultiPolygon25D:
00414 hasZValue = true;
00415 case QGis::WKBMultiPolygon:
00416 {
00417 currentWkbPtr += ( 1 + sizeof( int ) );
00418 int* npolys = ( int* )currentWkbPtr;
00419 int* nrings;
00420 int* npoints;
00421 currentWkbPtr += sizeof( int );
00422 for ( int index = 0; index < *npolys; ++index )
00423 {
00424 currentWkbPtr += ( 1 + sizeof( int ) );
00425 nrings = ( int* )currentWkbPtr;
00426 currentWkbPtr += sizeof( int );
00427 for ( int index2 = 0; index2 < *nrings; ++index2 )
00428 {
00429 if ( type != POINTS )
00430 {
00431 line = new Line3D();
00432 }
00433 npoints = ( int* )currentWkbPtr;
00434 currentWkbPtr += sizeof( int );
00435 for ( int index3 = 0; index3 < *npoints; ++index3 )
00436 {
00437 x = *(( double* )currentWkbPtr );
00438 currentWkbPtr += sizeof( double );
00439 y = *(( double* )currentWkbPtr );
00440 currentWkbPtr += sizeof( double );
00441 if ( hasZValue )
00442 {
00443 z = *(( double* )currentWkbPtr );
00444 currentWkbPtr += sizeof( double );
00445 }
00446 else
00447 {
00448 z = attributeValue;
00449 }
00450 if ( type == POINTS )
00451 {
00452
00453 mTriangulation->addPoint( new Point3D( x, y, z ) );
00454 }
00455 else
00456 {
00457 line->insertPoint( new Point3D( x, y, z ) );
00458 }
00459 }
00460 if ( type != POINTS )
00461 {
00462 mTriangulation->addLine( line, type == BREAK_LINES );
00463 }
00464 }
00465 }
00466 break;
00467 }
00468 default:
00469
00470 break;
00471 }
00472
00473 return 0;
00474 }
00475