1 /***************************************************************************
2     qgsrastercalcnode.cpp
3     ---------------------
4     begin                : October 2010
5     copyright            : (C) 2010 by Marco Hugentobler
6     email                : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************
8  *                                                                         *
9  *   This program is free software; you can redistribute it and/or modify  *
10  *   it under the terms of the GNU General Public License as published by  *
11  *   the Free Software Foundation; either version 2 of the License, or     *
12  *   (at your option) any later version.                                   *
13  *                                                                         *
14  ***************************************************************************/
15 #include "qgsrastercalcnode.h"
16 #include "qgsrasterblock.h"
17 #include "qgsrastermatrix.h"
18 
QgsRasterCalcNode(double number)19 QgsRasterCalcNode::QgsRasterCalcNode( double number )
20   : mNumber( number )
21 {
22 }
23 
QgsRasterCalcNode(QgsRasterMatrix * matrix)24 QgsRasterCalcNode::QgsRasterCalcNode( QgsRasterMatrix *matrix )
25   : mType( tMatrix )
26   , mMatrix( matrix )
27 {
28 
29 }
30 
QgsRasterCalcNode(Operator op,QgsRasterCalcNode * left,QgsRasterCalcNode * right)31 QgsRasterCalcNode::QgsRasterCalcNode( Operator op, QgsRasterCalcNode *left, QgsRasterCalcNode *right )
32   : mType( tOperator )
33   , mLeft( left )
34   , mRight( right )
35   , mOperator( op )
36 {
37 }
38 
QgsRasterCalcNode(QString functionName,QVector<QgsRasterCalcNode * > functionArgs)39 QgsRasterCalcNode::QgsRasterCalcNode( QString functionName, QVector <QgsRasterCalcNode *> functionArgs )
40   : mType( tFunction )
41   , mFunctionName( functionName )
42   , mFunctionArgs( functionArgs )
43 {
44 }
45 
QgsRasterCalcNode(const QString & rasterName)46 QgsRasterCalcNode::QgsRasterCalcNode( const QString &rasterName )
47   : mType( tRasterRef )
48   , mRasterName( rasterName )
49 {
50   if ( mRasterName.startsWith( '"' ) && mRasterName.endsWith( '"' ) )
51     mRasterName = mRasterName.mid( 1, mRasterName.size() - 2 );
52 }
53 
~QgsRasterCalcNode()54 QgsRasterCalcNode::~QgsRasterCalcNode()
55 {
56   delete mLeft;
57   delete mRight;
58 }
59 
calculate(QMap<QString,QgsRasterBlock * > & rasterData,QgsRasterMatrix & result,int row) const60 bool QgsRasterCalcNode::calculate( QMap<QString, QgsRasterBlock * > &rasterData, QgsRasterMatrix &result, int row ) const
61 {
62   //if type is raster ref: return a copy of the corresponding matrix
63 
64   //if type is operator, call the proper matrix operations
65   if ( mType == tRasterRef )
66   {
67     const QMap<QString, QgsRasterBlock *>::iterator it = rasterData.find( mRasterName );
68     if ( it == rasterData.end() )
69     {
70       QgsDebugMsg( QStringLiteral( "Error: could not find raster data for \"%1\"" ).arg( mRasterName ) );
71       return false;
72     }
73 
74     const int nRows = ( row >= 0 ? 1 : ( *it )->height() );
75     const int startRow = ( row >= 0 ? row : 0 );
76     const int endRow = startRow + nRows;
77     const int nCols = ( *it )->width();
78     const int nEntries = nCols * nRows;
79     double *data = new double[nEntries];
80 
81     //convert input raster values to double, also convert input no data to result no data
82 
83     int outRow = 0;
84     bool isNoData = false;
85     for ( int dataRow = startRow; dataRow < endRow ; ++dataRow, ++outRow )
86     {
87       for ( int dataCol = 0; dataCol < nCols; ++dataCol )
88       {
89         const double value = ( *it )->valueAndNoData( dataRow, dataCol, isNoData );
90         data[ dataCol + nCols * outRow] = isNoData ? result.nodataValue() : value;
91       }
92     }
93     result.setData( nCols, nRows, data, result.nodataValue() );
94     return true;
95   }
96   else if ( mType == tOperator )
97   {
98     QgsRasterMatrix leftMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() );
99     QgsRasterMatrix rightMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() );
100 
101     if ( !mLeft || !mLeft->calculate( rasterData, leftMatrix, row ) )
102     {
103       return false;
104     }
105     if ( mRight && !mRight->calculate( rasterData, rightMatrix, row ) )
106     {
107       return false;
108     }
109 
110     switch ( mOperator )
111     {
112       case opPLUS:
113         leftMatrix.add( rightMatrix );
114         break;
115       case opMINUS:
116         leftMatrix.subtract( rightMatrix );
117         break;
118       case opMUL:
119         leftMatrix.multiply( rightMatrix );
120         break;
121       case opDIV:
122         leftMatrix.divide( rightMatrix );
123         break;
124       case opPOW:
125         leftMatrix.power( rightMatrix );
126         break;
127       case opEQ:
128         leftMatrix.equal( rightMatrix );
129         break;
130       case opNE:
131         leftMatrix.notEqual( rightMatrix );
132         break;
133       case opGT:
134         leftMatrix.greaterThan( rightMatrix );
135         break;
136       case opLT:
137         leftMatrix.lesserThan( rightMatrix );
138         break;
139       case opGE:
140         leftMatrix.greaterEqual( rightMatrix );
141         break;
142       case opLE:
143         leftMatrix.lesserEqual( rightMatrix );
144         break;
145       case opAND:
146         leftMatrix.logicalAnd( rightMatrix );
147         break;
148       case opOR:
149         leftMatrix.logicalOr( rightMatrix );
150         break;
151       case opMIN:
152         leftMatrix.min( rightMatrix );
153         break;
154       case opMAX:
155         leftMatrix.max( rightMatrix );
156         break;
157       case opSQRT:
158         leftMatrix.squareRoot();
159         break;
160       case opSIN:
161         leftMatrix.sinus();
162         break;
163       case opCOS:
164         leftMatrix.cosinus();
165         break;
166       case opTAN:
167         leftMatrix.tangens();
168         break;
169       case opASIN:
170         leftMatrix.asinus();
171         break;
172       case opACOS:
173         leftMatrix.acosinus();
174         break;
175       case opATAN:
176         leftMatrix.atangens();
177         break;
178       case opSIGN:
179         leftMatrix.changeSign();
180         break;
181       case opLOG:
182         leftMatrix.log();
183         break;
184       case opLOG10:
185         leftMatrix.log10();
186         break;
187       case opABS:
188         leftMatrix.absoluteValue();
189         break;
190       default:
191         return false;
192     }
193     const int newNColumns = leftMatrix.nColumns();
194     const int newNRows = leftMatrix.nRows();
195     result.setData( newNColumns, newNRows, leftMatrix.takeData(), leftMatrix.nodataValue() );
196     return true;
197   }
198   else if ( mType == tNumber )
199   {
200     const size_t nEntries = static_cast<size_t>( result.nColumns() * result.nRows() );
201     double *data = new double[ nEntries ];
202     std::fill( data, data + nEntries, mNumber );
203     result.setData( result.nColumns(), 1, data, result.nodataValue() );
204 
205     return true;
206   }
207   else if ( mType == tMatrix )
208   {
209     const int nEntries = mMatrix->nColumns() * mMatrix->nRows();
210     double *data = new double[nEntries];
211     for ( int i = 0; i < nEntries; ++i )
212     {
213       data[i] = mMatrix->data()[i] == mMatrix->nodataValue() ? result.nodataValue() : mMatrix->data()[i];
214     }
215     result.setData( mMatrix->nColumns(), mMatrix->nRows(), data, result.nodataValue() );
216     return true;
217   }
218   else if ( mType == tFunction )
219   {
220     QVector <QgsRasterMatrix *> matrixContainer;
221     for ( int i = 0; i < mFunctionArgs.size(); ++i )
222     {
223       std::unique_ptr< QgsRasterMatrix > singleMatrix( new QgsRasterMatrix( result.nColumns(), result.nRows(), nullptr, result.nodataValue() ) );
224       if ( !mFunctionArgs.at( i ) || !mFunctionArgs.at( i )->calculate( rasterData, *singleMatrix, row ) )
225       {
226         return false;
227       }
228       matrixContainer.append( singleMatrix.release() );
229     }
230     evaluateFunction( matrixContainer, result );
231     return true;
232   }
233   return false;
234 }
235 
toString(bool cStyle) const236 QString QgsRasterCalcNode::toString( bool cStyle ) const
237 {
238   QString result;
239   QString left;
240   QString right;
241   if ( mLeft )
242     left = mLeft->toString( cStyle );
243   if ( mRight )
244     right = mRight->toString( cStyle );
245 
246   switch ( mType )
247   {
248     case tOperator:
249       switch ( mOperator )
250       {
251         case opPLUS:
252           result = QStringLiteral( "( %1 + %2 )" ).arg( left ).arg( right );
253           break;
254         case opMINUS:
255           result = QStringLiteral( "( %1 - %2 )" ).arg( left ).arg( right );
256           break;
257         case opSIGN:
258           result = QStringLiteral( "-%1" ).arg( left );
259           break;
260         case opMUL:
261           result = QStringLiteral( "%1 * %2" ).arg( left ).arg( right );
262           break;
263         case opDIV:
264           result = QStringLiteral( "%1 / %2" ).arg( left ).arg( right );
265           break;
266         case opPOW:
267           if ( cStyle )
268             result = QStringLiteral( "pow( %1, %2 )" ).arg( left ).arg( right );
269           else
270             result = QStringLiteral( "%1^%2" ).arg( left ).arg( right );
271           break;
272         case opEQ:
273           if ( cStyle )
274             result = QStringLiteral( "( float ) ( %1 == %2 )" ).arg( left ).arg( right );
275           else
276             result = QStringLiteral( "%1 = %2" ).arg( left ).arg( right );
277           break;
278         case opNE:
279           if ( cStyle )
280             result = QStringLiteral( "( float ) ( %1 != %2 )" ).arg( left ).arg( right );
281           else
282             result = QStringLiteral( "%1 != %2" ).arg( left ).arg( right );
283           break;
284         case opGT:
285           if ( cStyle )
286             result = QStringLiteral( "( float ) ( %1 > %2 )" ).arg( left ).arg( right );
287           else
288             result = QStringLiteral( "%1 > %2" ).arg( left ).arg( right );
289           break;
290         case opLT:
291           if ( cStyle )
292             result = QStringLiteral( "( float ) ( %1 < %2 )" ).arg( left ).arg( right );
293           else
294             result = QStringLiteral( "%1 < %2" ).arg( left ).arg( right );
295           break;
296         case opGE:
297           if ( cStyle )
298             result = QStringLiteral( "( float ) ( %1 >= %2 )" ).arg( left ).arg( right );
299           else
300             result = QStringLiteral( "%1 >= %2" ).arg( left ).arg( right );
301           break;
302         case opLE:
303           if ( cStyle )
304             result = QStringLiteral( "( float ) ( %1 <= %2 )" ).arg( left ).arg( right );
305           else
306             result = QStringLiteral( "%1 <= %2" ).arg( left ).arg( right );
307           break;
308         case opAND:
309           if ( cStyle )
310             result = QStringLiteral( "( float ) ( %1 && %2 )" ).arg( left ).arg( right );
311           else
312             result = QStringLiteral( "%1 AND %2" ).arg( left ).arg( right );
313           break;
314         case opOR:
315           if ( cStyle )
316             result = QStringLiteral( "( float ) ( %1 || %2 )" ).arg( left ).arg( right );
317           else
318             result = QStringLiteral( "%1 OR %2" ).arg( left ).arg( right );
319           break;
320         case opSQRT:
321           result = QStringLiteral( "sqrt( %1 )" ).arg( left );
322           break;
323         case opSIN:
324           result = QStringLiteral( "sin( %1 )" ).arg( left );
325           break;
326         case opCOS:
327           result = QStringLiteral( "cos( %1 )" ).arg( left );
328           break;
329         case opTAN:
330           result = QStringLiteral( "tan( %1 )" ).arg( left );
331           break;
332         case opASIN:
333           result = QStringLiteral( "asin( %1 )" ).arg( left );
334           break;
335         case opACOS:
336           result = QStringLiteral( "acos( %1 )" ).arg( left );
337           break;
338         case opATAN:
339           result = QStringLiteral( "atan( %1 )" ).arg( left );
340           break;
341         case opLOG:
342           result = QStringLiteral( "log( %1 )" ).arg( left );
343           break;
344         case opLOG10:
345           result = QStringLiteral( "log10( %1 )" ).arg( left );
346           break;
347         case opABS:
348           if ( cStyle )
349             result = QStringLiteral( "fabs( %1 )" ).arg( left );
350           else
351             // Call the floating point version
352             result = QStringLiteral( "abs( %1 )" ).arg( left );
353           break;
354         case opMIN:
355           if ( cStyle )
356             result = QStringLiteral( "min( ( float ) ( %1 ), ( float ) ( %2 ) )" ).arg( left ).arg( right );
357           else
358             result = QStringLiteral( "min( %1, %2 )" ).arg( left ).arg( right );
359           break;
360         case opMAX:
361           if ( cStyle )
362             result = QStringLiteral( "max( ( float ) ( %1 ), ( float ) ( %2 ) )" ).arg( left ).arg( right );
363           else
364             result = QStringLiteral( "max( %1, %2 )" ).arg( left ).arg( right );
365           break;
366         case opNONE:
367           break;
368       }
369       break;
370     case tRasterRef:
371       if ( cStyle )
372         result = QStringLiteral( "( float ) \"%1\"" ).arg( mRasterName );
373       else
374         result = QStringLiteral( "\"%1\"" ).arg( mRasterName );
375       break;
376     case tNumber:
377       result = QString::number( mNumber );
378       if ( cStyle )
379       {
380         result = QStringLiteral( "( float ) %1" ).arg( result );
381       }
382       break;
383     case tMatrix:
384       break;
385     case tFunction:
386       if ( mFunctionName == "if" )
387       {
388         const QString argOne = mFunctionArgs.at( 0 )->toString( cStyle );
389         const QString argTwo = mFunctionArgs.at( 1 )->toString( cStyle );
390         const QString argThree = mFunctionArgs.at( 2 )->toString( cStyle );
391         if ( cStyle )
392           result =  QStringLiteral( " ( %1 ) ? ( %2 ) : ( %3 ) " ).arg( argOne, argTwo, argThree );
393         else
394           result = QStringLiteral( "if( %1 , %2 , %3 )" ).arg( argOne, argTwo, argThree );
395       }
396       break;
397   }
398   return result;
399 }
400 
findNodes(const QgsRasterCalcNode::Type type) const401 QList<const QgsRasterCalcNode *> QgsRasterCalcNode::findNodes( const QgsRasterCalcNode::Type type ) const
402 {
403   QList<const QgsRasterCalcNode *> nodeList;
404   if ( mType == type )
405     nodeList.push_back( this );
406   if ( mLeft )
407     nodeList.append( mLeft->findNodes( type ) );
408   if ( mRight )
409     nodeList.append( mRight->findNodes( type ) );
410 
411   for ( QgsRasterCalcNode *node : mFunctionArgs )
412     nodeList.append( node->findNodes( type ) );
413 
414   return nodeList;
415 }
416 
parseRasterCalcString(const QString & str,QString & parserErrorMsg)417 QgsRasterCalcNode *QgsRasterCalcNode::parseRasterCalcString( const QString &str, QString &parserErrorMsg )
418 {
419   extern QgsRasterCalcNode *localParseRasterCalcString( const QString & str, QString & parserErrorMsg );
420   return localParseRasterCalcString( str, parserErrorMsg );
421 }
422 
referencedLayerNames()423 QStringList QgsRasterCalcNode::referencedLayerNames()
424 {
425   QStringList referencedRasters;
426 
427   QStringList rasterRef = this->cleanRasterReferences();
428   for ( const auto &i : rasterRef )
429   {
430     if ( referencedRasters.contains( i.mid( 0, i.lastIndexOf( "@" ) ) ) ) continue;
431     referencedRasters << i.mid( 0, i.lastIndexOf( "@" ) );
432   }
433 
434   return referencedRasters;
435 }
436 
cleanRasterReferences()437 QStringList QgsRasterCalcNode::cleanRasterReferences()
438 {
439   QStringList rasterReferences;
440   const QList<const QgsRasterCalcNode *> rasterRefNodes =  this->findNodes( QgsRasterCalcNode::Type::tRasterRef );
441 
442   for ( const QgsRasterCalcNode *r : rasterRefNodes )
443   {
444 
445     QString layerRef( r->toString() );
446     if ( layerRef.at( 0 ) == QLatin1String( "\"" ) && layerRef.at( layerRef.size() - 1 ) == QLatin1String( "\"" ) )
447     {
448       layerRef.remove( 0, 1 );
449       layerRef.chop( 1 );
450 
451     }
452     layerRef.remove( QChar( '\\' ), Qt::CaseInsensitive );
453     rasterReferences << layerRef;
454   }
455 
456   return rasterReferences;
457 }
458 
evaluateFunction(const QVector<QgsRasterMatrix * > & matrixVector,QgsRasterMatrix & result) const459 QgsRasterMatrix QgsRasterCalcNode::evaluateFunction( const QVector<QgsRasterMatrix *> &matrixVector, QgsRasterMatrix &result ) const
460 {
461 
462   if ( mFunctionName == "if" )
463   {
464     //scalar condition
465     if ( matrixVector.at( 0 )->isNumber() )
466     {
467       result = ( matrixVector.at( 0 )->data() ? * matrixVector.at( 1 ) : * matrixVector.at( 2 ) );
468       return result;
469     }
470     int nCols = matrixVector.at( 0 )->nColumns();
471     int nRows = matrixVector.at( 0 )->nRows();
472     int nEntries = nCols * nRows;
473     std::unique_ptr< double > dataResult( new double[nEntries] );
474     double *dataResultRawPtr =  dataResult.get();
475 
476     double *condition = matrixVector.at( 0 )->data();
477     double *firstOption = matrixVector.at( 1 )->data();
478     double *secondOption = matrixVector.at( 2 )->data();
479 
480     bool isFirstOptionNumber = matrixVector.at( 1 )->isNumber();
481     bool isSecondCOptionNumber = matrixVector.at( 2 )->isNumber();
482     double noDataValueCondition = matrixVector.at( 0 )->nodataValue();
483 
484     for ( int i = 0; i < nEntries; ++i )
485     {
486       if ( condition[i] == noDataValueCondition )
487       {
488         dataResultRawPtr[i] = result.nodataValue();
489         continue;
490       }
491       else if ( condition[i] != 0 )
492       {
493         dataResultRawPtr[i] = isFirstOptionNumber ? firstOption[0] : firstOption[i];
494         continue;
495       }
496       dataResultRawPtr[i] = isSecondCOptionNumber ? secondOption[0] : secondOption[i];
497     }
498 
499     result.setData( nCols, nRows, dataResult.release(), result.nodataValue() );
500   }
501   return result;
502 }
503