Quantum GIS API Documentation
1.8
|
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 }