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