Quantum GIS API Documentation  1.8
src/analysis/interpolation/qgstininterpolator.cpp
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines