Quantum GIS API Documentation
1.7.4
|
00001 /*************************************************************************** 00002 qgsrastermatrix.cpp 00003 ------------------- 00004 begin : 2010-10-23 00005 copyright : (C) 20010 by Marco Hugentobler 00006 email : marco dot hugentobler at sourcepole 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 "qgsrastermatrix.h" 00019 #include <string.h> 00020 00021 #include <cmath> 00022 00023 QgsRasterMatrix::QgsRasterMatrix(): mColumns( 0 ), mRows( 0 ), mData( 0 ) 00024 { 00025 } 00026 00027 QgsRasterMatrix::QgsRasterMatrix( int nCols, int nRows, float* data, double nodataValue ): 00028 mColumns( nCols ), mRows( nRows ), mData( data ), mNodataValue( nodataValue ) 00029 { 00030 } 00031 00032 QgsRasterMatrix::QgsRasterMatrix( const QgsRasterMatrix& m ): mColumns( 0 ), mRows( 0 ), mData( 0 ) 00033 { 00034 operator=( m ); 00035 } 00036 00037 QgsRasterMatrix::~QgsRasterMatrix() 00038 { 00039 delete[] mData; 00040 } 00041 00042 QgsRasterMatrix& QgsRasterMatrix::operator=( const QgsRasterMatrix & m ) 00043 { 00044 delete[] mData; 00045 mColumns = m.nColumns(); 00046 mRows = m.nRows(); 00047 int nEntries = mColumns * mRows; 00048 mData = new float[nEntries]; 00049 memcpy( mData, m.mData, sizeof( float ) * nEntries ); 00050 mNodataValue = m.nodataValue(); 00051 return *this; 00052 } 00053 00054 void QgsRasterMatrix::setData( int cols, int rows, float* data, double nodataValue ) 00055 { 00056 delete[] mData; 00057 mColumns = cols; 00058 mRows = rows; 00059 mData = data; 00060 mNodataValue = nodataValue; 00061 } 00062 00063 float* QgsRasterMatrix::takeData() 00064 { 00065 float* data = mData; 00066 mData = 0; mColumns = 0; mRows = 0; 00067 return data; 00068 } 00069 00070 bool QgsRasterMatrix::add( const QgsRasterMatrix& other ) 00071 { 00072 return twoArgumentOperation( opPLUS, other ); 00073 } 00074 00075 bool QgsRasterMatrix::subtract( const QgsRasterMatrix& other ) 00076 { 00077 return twoArgumentOperation( opMINUS, other ); 00078 } 00079 00080 bool QgsRasterMatrix::multiply( const QgsRasterMatrix& other ) 00081 { 00082 return twoArgumentOperation( opMUL, other ); 00083 } 00084 00085 bool QgsRasterMatrix::divide( const QgsRasterMatrix& other ) 00086 { 00087 return twoArgumentOperation( opDIV, other ); 00088 } 00089 00090 bool QgsRasterMatrix::power( const QgsRasterMatrix& other ) 00091 { 00092 return twoArgumentOperation( opPOW, other ); 00093 } 00094 00095 bool QgsRasterMatrix::equal( const QgsRasterMatrix& other ) 00096 { 00097 return twoArgumentOperation( opEQ, other ); 00098 } 00099 00100 bool QgsRasterMatrix::notEqual( const QgsRasterMatrix& other ) 00101 { 00102 return twoArgumentOperation( opNE, other ); 00103 } 00104 00105 bool QgsRasterMatrix::greaterThan( const QgsRasterMatrix& other ) 00106 { 00107 return twoArgumentOperation( opGT, other ); 00108 } 00109 00110 bool QgsRasterMatrix::lesserThan( const QgsRasterMatrix& other ) 00111 { 00112 return twoArgumentOperation( opLT, other ); 00113 } 00114 00115 bool QgsRasterMatrix::greaterEqual( const QgsRasterMatrix& other ) 00116 { 00117 return twoArgumentOperation( opGE, other ); 00118 } 00119 00120 bool QgsRasterMatrix::lesserEqual( const QgsRasterMatrix& other ) 00121 { 00122 return twoArgumentOperation( opLE, other ); 00123 } 00124 00125 bool QgsRasterMatrix::logicalAnd( const QgsRasterMatrix& other ) 00126 { 00127 return twoArgumentOperation( opAND, other ); 00128 } 00129 00130 bool QgsRasterMatrix::logicalOr( const QgsRasterMatrix& other ) 00131 { 00132 return twoArgumentOperation( opOR, other ); 00133 } 00134 00135 bool QgsRasterMatrix::squareRoot() 00136 { 00137 return oneArgumentOperation( opSQRT ); 00138 } 00139 00140 bool QgsRasterMatrix::sinus() 00141 { 00142 return oneArgumentOperation( opSIN ); 00143 } 00144 00145 bool QgsRasterMatrix::asinus() 00146 { 00147 return oneArgumentOperation( opASIN ); 00148 } 00149 00150 bool QgsRasterMatrix::cosinus() 00151 { 00152 return oneArgumentOperation( opCOS ); 00153 } 00154 00155 bool QgsRasterMatrix::acosinus() 00156 { 00157 return oneArgumentOperation( opACOS ); 00158 } 00159 00160 bool QgsRasterMatrix::tangens() 00161 { 00162 return oneArgumentOperation( opTAN ); 00163 } 00164 00165 bool QgsRasterMatrix::atangens() 00166 { 00167 return oneArgumentOperation( opATAN ); 00168 } 00169 00170 bool QgsRasterMatrix::changeSign() 00171 { 00172 return oneArgumentOperation( opSIGN ); 00173 } 00174 00175 bool QgsRasterMatrix::oneArgumentOperation( OneArgOperator op ) 00176 { 00177 if ( !mData ) 00178 { 00179 return false; 00180 } 00181 00182 int nEntries = mColumns * mRows; 00183 double value; 00184 for ( int i = 0; i < nEntries; ++i ) 00185 { 00186 value = mData[i]; 00187 if ( value != mNodataValue ) 00188 { 00189 switch ( op ) 00190 { 00191 case opSQRT: 00192 if ( value < 0 ) //no complex numbers 00193 { 00194 mData[i] = static_cast<float>( mNodataValue ); 00195 } 00196 else 00197 { 00198 mData[i] = static_cast<float>( sqrt( value ) ); 00199 } 00200 break; 00201 case opSIN: 00202 mData[i] = static_cast<float>( sin( value ) ); 00203 break; 00204 case opCOS: 00205 mData[i] = static_cast<float>( cos( value ) ); 00206 break; 00207 case opTAN: 00208 mData[i] = static_cast<float>( tan( value ) ); 00209 break; 00210 case opASIN: 00211 mData[i] = static_cast<float>( asin( value ) ); 00212 break; 00213 case opACOS: 00214 mData[i] = static_cast<float>( acos( value ) ); 00215 break; 00216 case opATAN: 00217 mData[i] = static_cast<float>( atan( value ) ); 00218 break; 00219 case opSIGN: 00220 mData[i] = -value; 00221 } 00222 } 00223 } 00224 return true; 00225 } 00226 00227 bool QgsRasterMatrix::twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix& other ) 00228 { 00229 if ( isNumber() && other.isNumber() ) //operation on two 1x1 matrices 00230 { 00231 //operations with nodata values always generate nodata 00232 if ( mData[0] == mNodataValue || other.number() == other.nodataValue() ) 00233 { 00234 mData[0] = static_cast<float>( mNodataValue ); 00235 return true; 00236 } 00237 switch ( op ) 00238 { 00239 case opPLUS: 00240 mData[0] = static_cast<float>( number() + other.number() ); 00241 break; 00242 case opMINUS: 00243 mData[0] = static_cast<float>( number() - other.number() ); 00244 break; 00245 case opMUL: 00246 mData[0] = static_cast<float>( number() * other.number() ); 00247 break; 00248 case opDIV: 00249 if ( other.number() == 0 ) 00250 { 00251 mData[0] = static_cast<float>( mNodataValue ); 00252 } 00253 else 00254 { 00255 mData[0] = static_cast<float>( number() / other.number() ); 00256 } 00257 break; 00258 case opPOW: 00259 if ( !testPowerValidity( mData[0], ( float ) other.number() ) ) 00260 { 00261 mData[0] = static_cast<float>( mNodataValue ); 00262 } 00263 else 00264 { 00265 mData[0] = pow( mData[0], ( float ) other.number() ); 00266 } 00267 break; 00268 case opEQ: 00269 mData[0] = mData[0] == other.number() ? 1.0f : 0.0f; 00270 break; 00271 case opNE: 00272 mData[0] = mData[0] == other.number() ? 0.0f : 1.0f; 00273 break; 00274 case opGT: 00275 mData[0] = mData[0] > other.number() ? 1.0f : 0.0f; 00276 break; 00277 case opLT: 00278 mData[0] = mData[0] < other.number() ? 1.0f : 0.0f; 00279 break; 00280 case opGE: 00281 mData[0] = mData[0] >= other.number() ? 1.0f : 0.0f; 00282 break; 00283 case opLE: 00284 mData[0] = mData[0] <= other.number() ? 1.0f : 0.0f; 00285 break; 00286 case opAND: 00287 mData[0] = mData[0] && other.number() ? 1.0f : 0.0f; 00288 break; 00289 case opOR: 00290 mData[0] = mData[0] || other.number() ? 1.0f : 0.0f; 00291 break; 00292 } 00293 return true; 00294 } 00295 00296 //two matrices 00297 if ( !isNumber() && !other.isNumber() ) 00298 { 00299 float* matrix = other.mData; 00300 int nEntries = mColumns * mRows; 00301 double value1, value2; 00302 00303 for ( int i = 0; i < nEntries; ++i ) 00304 { 00305 value1 = mData[i]; value2 = matrix[i]; 00306 if ( value1 == mNodataValue || value2 == other.mNodataValue ) 00307 { 00308 mData[i] = static_cast<float>( mNodataValue ); 00309 } 00310 else 00311 { 00312 switch ( op ) 00313 { 00314 case opPLUS: 00315 mData[i] = static_cast<float>( value1 + value2 ); 00316 break; 00317 case opMINUS: 00318 mData[i] = static_cast<float>( value1 - value2 ); 00319 break; 00320 case opMUL: 00321 mData[i] = static_cast<float>( value1 * value2 ); 00322 break; 00323 case opDIV: 00324 if ( value2 == 0 ) 00325 { 00326 mData[i] = static_cast<float>( mNodataValue ); 00327 } 00328 else 00329 { 00330 mData[i] = static_cast<float>( value1 / value2 ); 00331 } 00332 break; 00333 case opPOW: 00334 if ( !testPowerValidity( value1, value2 ) ) 00335 { 00336 mData[i] = static_cast<float>( mNodataValue ); 00337 } 00338 else 00339 { 00340 mData[i] = static_cast<float>( pow( value1, value2 ) ); 00341 } 00342 break; 00343 case opEQ: 00344 mData[i] = value1 == value2 ? 1.0f : 0.0f; 00345 break; 00346 case opNE: 00347 mData[i] = value1 == value2 ? 0.0f : 1.0f; 00348 break; 00349 case opGT: 00350 mData[i] = value1 > value2 ? 1.0f : 0.0f; 00351 break; 00352 case opLT: 00353 mData[i] = value1 < value2 ? 1.0f : 0.0f; 00354 break; 00355 case opGE: 00356 mData[i] = value1 >= value2 ? 1.0f : 0.0f; 00357 break; 00358 case opLE: 00359 mData[i] = value1 <= value2 ? 1.0f : 0.0f; 00360 break; 00361 case opAND: 00362 mData[i] = value1 && value2 ? 1.0f : 0.0f; 00363 break; 00364 case opOR: 00365 mData[i] = value1 || value2 ? 1.0f : 0.0f; 00366 break; 00367 } 00368 } 00369 } 00370 return true; 00371 } 00372 00373 //this matrix is a single number and the other one a real matrix 00374 if ( isNumber() ) 00375 { 00376 float* matrix = other.mData; 00377 int nEntries = other.nColumns() * other.nRows(); 00378 double value = mData[0]; 00379 delete[] mData; 00380 mData = new float[nEntries]; mColumns = other.nColumns(); mRows = other.nRows(); 00381 mNodataValue = other.nodataValue(); 00382 00383 if ( value == mNodataValue ) 00384 { 00385 for ( int i = 0; i < nEntries; ++i ) 00386 { 00387 mData[i] = static_cast<float>( mNodataValue ); 00388 } 00389 return true; 00390 } 00391 00392 for ( int i = 0; i < nEntries; ++i ) 00393 { 00394 if ( matrix[i] == other.mNodataValue ) 00395 { 00396 mData[i] = static_cast<float>( mNodataValue ); 00397 continue; 00398 } 00399 00400 switch ( op ) 00401 { 00402 case opPLUS: 00403 mData[i] = static_cast<float>( value + matrix[i] ); 00404 break; 00405 case opMINUS: 00406 mData[i] = static_cast<float>( value - matrix[i] ); 00407 break; 00408 case opMUL: 00409 mData[i] = static_cast<float>( value * matrix[i] ); 00410 break; 00411 case opDIV: 00412 if ( matrix[i] == 0 ) 00413 { 00414 mData[i] = static_cast<float>( mNodataValue ); 00415 } 00416 else 00417 { 00418 mData[i] = static_cast<float>( value / matrix[i] ); 00419 } 00420 break; 00421 case opPOW: 00422 if ( !testPowerValidity( value, matrix[i] ) ) 00423 { 00424 mData[i] = static_cast<float>( mNodataValue ); 00425 } 00426 else 00427 { 00428 mData[i] = pow( static_cast<float>( value ), matrix[i] ); 00429 } 00430 break; 00431 case opEQ: 00432 mData[i] = value == matrix[i] ? 1.0f : 0.0f; 00433 break; 00434 case opNE: 00435 mData[i] = value == matrix[i] ? 0.0f : 1.0f; 00436 break; 00437 case opGT: 00438 mData[i] = value > matrix[i] ? 1.0f : 0.0f; 00439 break; 00440 case opLT: 00441 mData[i] = value < matrix[i] ? 1.0f : 0.0f; 00442 break; 00443 case opGE: 00444 mData[i] = value >= matrix[i] ? 1.0f : 0.0f; 00445 break; 00446 case opLE: 00447 mData[i] = value <= matrix[i] ? 1.0f : 0.0f; 00448 break; 00449 case opAND: 00450 mData[i] = value && matrix[i] ? 1.0f : 0.0f; 00451 break; 00452 case opOR: 00453 mData[i] = value || matrix[i] ? 1.0f : 0.0f; 00454 break; 00455 } 00456 } 00457 return true; 00458 } 00459 else //this matrix is a real matrix and the other a number 00460 { 00461 double value = other.number(); 00462 int nEntries = mColumns * mRows; 00463 00464 if ( other.number() == other.mNodataValue ) 00465 { 00466 for ( int i = 0; i < nEntries; ++i ) 00467 { 00468 mData[i] = static_cast<float>( mNodataValue ); 00469 } 00470 return true; 00471 } 00472 00473 for ( int i = 0; i < nEntries; ++i ) 00474 { 00475 if ( mData[i] == mNodataValue ) 00476 { 00477 continue; 00478 } 00479 00480 switch ( op ) 00481 { 00482 case opPLUS: 00483 mData[i] = static_cast<float>( mData[i] + value ); 00484 break; 00485 case opMINUS: 00486 mData[i] = static_cast<float>( mData[i] - value ); 00487 break; 00488 case opMUL: 00489 mData[i] = static_cast<float>( mData[i] * value ); 00490 break; 00491 case opDIV: 00492 if ( value == 0 ) 00493 { 00494 mData[i] = static_cast<float>( mNodataValue ); 00495 } 00496 else 00497 { 00498 mData[i] = static_cast<float>( mData[i] / value ); 00499 } 00500 break; 00501 case opPOW: 00502 if ( !testPowerValidity( mData[i], value ) ) 00503 { 00504 mData[i] = static_cast<float>( mNodataValue ); 00505 } 00506 else 00507 { 00508 mData[i] = pow( mData[i], ( float ) value ); 00509 } 00510 break; 00511 case opEQ: 00512 mData[i] = mData[i] == value ? 1.0f : 0.0f; 00513 break; 00514 case opNE: 00515 mData[i] = mData[i] == value ? 0.0f : 1.0f; 00516 break; 00517 case opGT: 00518 mData[i] = mData[i] > value ? 1.0f : 0.0f; 00519 break; 00520 case opLT: 00521 mData[i] = mData[i] < value ? 1.0f : 0.0f; 00522 break; 00523 case opGE: 00524 mData[i] = mData[i] >= value ? 1.0f : 0.0f; 00525 break; 00526 case opLE: 00527 mData[i] = mData[i] <= value ? 1.0f : 0.0f; 00528 break; 00529 case opAND: 00530 mData[i] = mData[i] && value ? 1.0f : 0.0f; 00531 break; 00532 case opOR: 00533 mData[i] = mData[i] || value ? 1.0f : 0.0f; 00534 break; 00535 } 00536 } 00537 return true; 00538 } 00539 } 00540 00541 bool QgsRasterMatrix::testPowerValidity( double base, double power ) 00542 { 00543 if (( base == 0 && power < 0 ) || ( power < 0 && ( power - floor( power ) ) > 0 ) ) 00544 { 00545 return false; 00546 } 00547 return true; 00548 }