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