Quantum GIS API Documentation  1.8
src/analysis/interpolation/qgsinterpolator.cpp
Go to the documentation of this file.
00001 /***************************************************************************
00002                               qgsinterpolator.cpp
00003                               -------------------
00004   begin                : Marco 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 "qgsinterpolator.h"
00019 #include "qgsvectordataprovider.h"
00020 #include "qgsvectorlayer.h"
00021 #include "qgsgeometry.h"
00022 
00023 QgsInterpolator::QgsInterpolator( const QList<LayerData>& layerData ): mDataIsCached( false ), mLayerData( layerData )
00024 {
00025 
00026 }
00027 
00028 QgsInterpolator::QgsInterpolator()
00029 {
00030 
00031 }
00032 
00033 QgsInterpolator::~QgsInterpolator()
00034 {
00035 
00036 }
00037 
00038 int QgsInterpolator::cacheBaseData()
00039 {
00040   if ( mLayerData.size() < 1 )
00041   {
00042     return 0;
00043   }
00044 
00045   //reserve initial memory for 100000 vertices
00046   mCachedBaseData.clear();
00047   mCachedBaseData.reserve( 100000 );
00048 
00049   QList<LayerData>::iterator v_it = mLayerData.begin();
00050 
00051   for ( ; v_it != mLayerData.end(); ++v_it )
00052   {
00053     if ( v_it->vectorLayer == 0 )
00054     {
00055       continue;
00056     }
00057 
00058     QgsVectorDataProvider* provider = v_it->vectorLayer->dataProvider();
00059     if ( !provider )
00060     {
00061       return 2;
00062     }
00063 
00064     QgsAttributeList attList;
00065     if ( !v_it->zCoordInterpolation )
00066     {
00067       attList.push_back( v_it->interpolationAttribute );
00068     }
00069 
00070     provider->select( attList );
00071 
00072     QgsFeature theFeature;
00073     double attributeValue = 0.0;
00074     bool attributeConversionOk = false;
00075 
00076     while ( provider->nextFeature( theFeature ) )
00077     {
00078       if ( !v_it->zCoordInterpolation )
00079       {
00080         QgsAttributeMap attMap = theFeature.attributeMap();
00081         QgsAttributeMap::const_iterator att_it = attMap.find( v_it->interpolationAttribute );
00082         if ( att_it == attMap.end() ) //attribute not found, something must be wrong (e.g. NULL value)
00083         {
00084           continue;
00085         }
00086         attributeValue = att_it.value().toDouble( &attributeConversionOk );
00087         if ( !attributeConversionOk || qIsNaN( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
00088         {
00089           continue;
00090         }
00091       }
00092 
00093       if ( addVerticesToCache( theFeature.geometry(), v_it->zCoordInterpolation, attributeValue ) != 0 )
00094       {
00095         return 3;
00096       }
00097     }
00098   }
00099 
00100   return 0;
00101 }
00102 
00103 int QgsInterpolator::addVerticesToCache( QgsGeometry* geom, bool zCoord, double attributeValue )
00104 {
00105   if ( !geom )
00106   {
00107     return 1;
00108   }
00109 
00110   bool hasZValue = false;
00111   unsigned char* currentWkbPtr = geom->asWkb();
00112   vertexData theVertex; //the current vertex
00113 
00114   QGis::WkbType wkbType = geom->wkbType();
00115   switch ( wkbType )
00116   {
00117     case QGis::WKBPoint25D:
00118       hasZValue = true;
00119     case QGis::WKBPoint:
00120     {
00121       currentWkbPtr += ( 1 + sizeof( int ) );
00122       theVertex.x = *(( double * )( currentWkbPtr ) );
00123       currentWkbPtr += sizeof( double );
00124       theVertex.y = *(( double * )( currentWkbPtr ) );
00125       if ( zCoord && hasZValue )
00126       {
00127         currentWkbPtr += sizeof( double );
00128         theVertex.z = *(( double * )( currentWkbPtr ) );
00129       }
00130       else
00131       {
00132         theVertex.z = attributeValue;
00133       }
00134       mCachedBaseData.push_back( theVertex );
00135       break;
00136     }
00137     case QGis::WKBLineString25D:
00138       hasZValue = true;
00139     case QGis::WKBLineString:
00140     {
00141       currentWkbPtr += ( 1 + sizeof( int ) );
00142       int* npoints = ( int* )currentWkbPtr;
00143       currentWkbPtr += sizeof( int );
00144       for ( int index = 0; index < *npoints; ++index )
00145       {
00146         theVertex.x = *(( double * )( currentWkbPtr ) );
00147         currentWkbPtr += sizeof( double );
00148         theVertex.y = *(( double * )( currentWkbPtr ) );
00149         currentWkbPtr += sizeof( double );
00150         if ( zCoord && hasZValue ) //skip z-coordinate for 25D geometries
00151         {
00152           theVertex.z = *(( double * )( currentWkbPtr ) );
00153           currentWkbPtr += sizeof( double );
00154         }
00155         else
00156         {
00157           theVertex.z = attributeValue;
00158         }
00159         mCachedBaseData.push_back( theVertex );
00160       }
00161       break;
00162     }
00163 #if 0
00164     case QGis::WKBPolygon25D:
00165       hasZValue = true;
00166     case QGis::WKBPolygon:
00167     {
00168       int* nrings = ( int* )( mGeometry + 5 );
00169       int* npoints;
00170       unsigned char* ptr = mGeometry + 9;
00171       for ( int index = 0; index < *nrings; ++index )
00172       {
00173         npoints = ( int* )ptr;
00174         ptr += sizeof( int );
00175         for ( int index2 = 0; index2 < *npoints; ++index2 )
00176         {
00177           tempx = ( double* )ptr;
00178           ptr += sizeof( double );
00179           tempy = ( double* )ptr;
00180           if ( point.sqrDist( *tempx, *tempy ) < actdist )
00181           {
00182             x = *tempx;
00183             y = *tempy;
00184             actdist = point.sqrDist( *tempx, *tempy );
00185             vertexnr = vertexcounter;
00186             //assign the rubber band indices
00187             if ( index2 == 0 )
00188             {
00189               beforeVertex = vertexcounter + ( *npoints - 2 );
00190               afterVertex = vertexcounter + 1;
00191             }
00192             else if ( index2 == ( *npoints - 1 ) )
00193             {
00194               beforeVertex = vertexcounter - 1;
00195               afterVertex = vertexcounter - ( *npoints - 2 );
00196             }
00197             else
00198             {
00199               beforeVertex = vertexcounter - 1;
00200               afterVertex = vertexcounter + 1;
00201             }
00202           }
00203           ptr += sizeof( double );
00204           if ( hasZValue ) //skip z-coordinate for 25D geometries
00205           {
00206             ptr += sizeof( double );
00207           }
00208           ++vertexcounter;
00209         }
00210       }
00211       break;
00212     }
00213     case QGis::WKBMultiPoint25D:
00214       hasZValue = true;
00215     case QGis::WKBMultiPoint:
00216     {
00217       unsigned char* ptr = mGeometry + 5;
00218       int* npoints = ( int* )ptr;
00219       ptr += sizeof( int );
00220       for ( int index = 0; index < *npoints; ++index )
00221       {
00222         ptr += ( 1 + sizeof( int ) ); //skip endian and point type
00223         tempx = ( double* )ptr;
00224         tempy = ( double* )( ptr + sizeof( double ) );
00225         if ( point.sqrDist( *tempx, *tempy ) < actdist )
00226         {
00227           x = *tempx;
00228           y = *tempy;
00229           actdist = point.sqrDist( *tempx, *tempy );
00230           vertexnr = index;
00231         }
00232         ptr += ( 2 * sizeof( double ) );
00233         if ( hasZValue ) //skip z-coordinate for 25D geometries
00234         {
00235           ptr += sizeof( double );
00236         }
00237       }
00238       break;
00239     }
00240     case QGis::WKBMultiLineString25D:
00241       hasZValue = true;
00242     case QGis::WKBMultiLineString:
00243     {
00244       unsigned char* ptr = mGeometry + 5;
00245       int* nlines = ( int* )ptr;
00246       int* npoints = 0;
00247       ptr += sizeof( int );
00248       for ( int index = 0; index < *nlines; ++index )
00249       {
00250         ptr += ( sizeof( int ) + 1 );
00251         npoints = ( int* )ptr;
00252         ptr += sizeof( int );
00253         for ( int index2 = 0; index2 < *npoints; ++index2 )
00254         {
00255           tempx = ( double* )ptr;
00256           ptr += sizeof( double );
00257           tempy = ( double* )ptr;
00258           ptr += sizeof( double );
00259           if ( point.sqrDist( *tempx, *tempy ) < actdist )
00260           {
00261             x = *tempx;
00262             y = *tempy;
00263             actdist = point.sqrDist( *tempx, *tempy );
00264             vertexnr = vertexcounter;
00265 
00266             if ( index2 == 0 )//assign the rubber band indices
00267             {
00268               beforeVertex = -1;
00269             }
00270             else
00271             {
00272               beforeVertex = vertexnr - 1;
00273             }
00274             if ( index2 == ( *npoints ) - 1 )
00275             {
00276               afterVertex = -1;
00277             }
00278             else
00279             {
00280               afterVertex = vertexnr + 1;
00281             }
00282           }
00283           if ( hasZValue ) //skip z-coordinate for 25D geometries
00284           {
00285             ptr += sizeof( double );
00286           }
00287           ++vertexcounter;
00288         }
00289       }
00290       break;
00291     }
00292     case QGis::WKBMultiPolygon25D:
00293       hasZValue = true;
00294     case QGis::WKBMultiPolygon:
00295     {
00296       unsigned char* ptr = mGeometry + 5;
00297       int* npolys = ( int* )ptr;
00298       int* nrings;
00299       int* npoints;
00300       ptr += sizeof( int );
00301       for ( int index = 0; index < *npolys; ++index )
00302       {
00303         ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type
00304         nrings = ( int* )ptr;
00305         ptr += sizeof( int );
00306         for ( int index2 = 0; index2 < *nrings; ++index2 )
00307         {
00308           npoints = ( int* )ptr;
00309           ptr += sizeof( int );
00310           for ( int index3 = 0; index3 < *npoints; ++index3 )
00311           {
00312             tempx = ( double* )ptr;
00313             ptr += sizeof( double );
00314             tempy = ( double* )ptr;
00315             if ( point.sqrDist( *tempx, *tempy ) < actdist )
00316             {
00317               x = *tempx;
00318               y = *tempy;
00319               actdist = point.sqrDist( *tempx, *tempy );
00320               vertexnr = vertexcounter;
00321 
00322               //assign the rubber band indices
00323               if ( index3 == 0 )
00324               {
00325                 beforeVertex = vertexcounter + ( *npoints - 2 );
00326                 afterVertex = vertexcounter + 1;
00327               }
00328               else if ( index3 == ( *npoints - 1 ) )
00329               {
00330                 beforeVertex = vertexcounter - 1;
00331                 afterVertex = vertexcounter - ( *npoints - 2 );
00332               }
00333               else
00334               {
00335                 beforeVertex = vertexcounter - 1;
00336                 afterVertex = vertexcounter + 1;
00337               }
00338             }
00339             ptr += sizeof( double );
00340             if ( hasZValue ) //skip z-coordinate for 25D geometries
00341             {
00342               ptr += sizeof( double );
00343             }
00344             ++vertexcounter;
00345           }
00346         }
00347       }
00348       break;
00349     }
00350 #endif //0
00351     default:
00352       break;
00353   }
00354   mDataIsCached = true;
00355   return 0;
00356 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines