QGIS API Documentation  3.13.0-Master (b73bd58cfb)
qgsrastermatrix.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrastermatrix.cpp
3  -------------------
4  begin : 2010-10-23
5  copyright : (C) 20010 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7 ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgsrastermatrix.h"
19 #include <cstring>
20 #include <cmath>
21 #include <algorithm>
22 
23 QgsRasterMatrix::QgsRasterMatrix( int nCols, int nRows, double *data, double nodataValue )
24  : mColumns( nCols )
25  , mRows( nRows )
26  , mData( data )
27  , mNodataValue( nodataValue )
28 {
29 }
30 
32 {
33  operator=( m );
34 }
35 
37 {
38  delete[] mData;
39 }
40 
42 {
43  delete[] mData;
44  mColumns = m.nColumns();
45  mRows = m.nRows();
46  int nEntries = mColumns * mRows;
47  mData = new double[nEntries];
48  memcpy( mData, m.mData, sizeof( double ) * nEntries );
49  mNodataValue = m.nodataValue();
50  return *this;
51 }
52 
53 void QgsRasterMatrix::setData( int cols, int rows, double *data, double nodataValue )
54 {
55  delete[] mData;
56  mColumns = cols;
57  mRows = rows;
58  mData = data;
59  mNodataValue = nodataValue;
60 }
61 
63 {
64  double *data = mData;
65  mData = nullptr;
66  mColumns = 0;
67  mRows = 0;
68  return data;
69 }
70 
72 {
73  return twoArgumentOperation( opPLUS, other );
74 }
75 
77 {
78  return twoArgumentOperation( opMINUS, other );
79 }
80 
82 {
83  return twoArgumentOperation( opMUL, other );
84 }
85 
87 {
88  return twoArgumentOperation( opDIV, other );
89 }
90 
92 {
93  return twoArgumentOperation( opPOW, other );
94 }
95 
97 {
98  return twoArgumentOperation( opEQ, other );
99 }
100 
102 {
103  return twoArgumentOperation( opNE, other );
104 }
105 
107 {
108  return twoArgumentOperation( opGT, other );
109 }
110 
112 {
113  return twoArgumentOperation( opLT, other );
114 }
115 
117 {
118  return twoArgumentOperation( opGE, other );
119 }
120 
122 {
123  return twoArgumentOperation( opLE, other );
124 }
125 
127 {
128  return twoArgumentOperation( opAND, other );
129 }
130 
132 {
133  return twoArgumentOperation( opOR, other );
134 }
135 
137 {
138  return twoArgumentOperation( opMAX, other );
139 }
140 
142 {
143  return twoArgumentOperation( opMIN, other );
144 }
145 
147 {
148  return oneArgumentOperation( opSQRT );
149 }
150 
152 {
153  return oneArgumentOperation( opSIN );
154 }
155 
157 {
158  return oneArgumentOperation( opASIN );
159 }
160 
162 {
163  return oneArgumentOperation( opCOS );
164 }
165 
167 {
168  return oneArgumentOperation( opACOS );
169 }
170 
172 {
173  return oneArgumentOperation( opTAN );
174 }
175 
177 {
178  return oneArgumentOperation( opATAN );
179 }
180 
182 {
183  return oneArgumentOperation( opSIGN );
184 }
185 
187 {
188  return oneArgumentOperation( opLOG );
189 }
190 
192 {
193  return oneArgumentOperation( opLOG10 );
194 }
195 
197 {
198  return oneArgumentOperation( opABS );
199 }
200 
201 bool QgsRasterMatrix::oneArgumentOperation( OneArgOperator op )
202 {
203  if ( !mData )
204  {
205  return false;
206  }
207 
208  int nEntries = mColumns * mRows;
209  double value;
210  for ( int i = 0; i < nEntries; ++i )
211  {
212  value = mData[i];
213  if ( value != mNodataValue )
214  {
215  switch ( op )
216  {
217  case opSQRT:
218  if ( value < 0 ) //no complex numbers
219  {
220  mData[i] = mNodataValue;
221  }
222  else
223  {
224  mData[i] = std::sqrt( value );
225  }
226  break;
227  case opSIN:
228  mData[i] = std::sin( value );
229  break;
230  case opCOS:
231  mData[i] = std::cos( value );
232  break;
233  case opTAN:
234  mData[i] = std::tan( value );
235  break;
236  case opASIN:
237  mData[i] = std::asin( value );
238  break;
239  case opACOS:
240  mData[i] = std::acos( value );
241  break;
242  case opATAN:
243  mData[i] = std::atan( value );
244  break;
245  case opSIGN:
246  mData[i] = -value;
247  break;
248  case opLOG:
249  if ( value <= 0 )
250  {
251  mData[i] = mNodataValue;
252  }
253  else
254  {
255  mData[i] = ::log( value );
256  }
257  break;
258  case opLOG10:
259  if ( value <= 0 )
260  {
261  mData[i] = mNodataValue;
262  }
263  else
264  {
265  mData[i] = ::log10( value );
266  }
267  break;
268  case opABS:
269  mData[i] = ::fabs( value );
270  break;
271  }
272  }
273  }
274  return true;
275 }
276 
277 double QgsRasterMatrix::calculateTwoArgumentOp( TwoArgOperator op, double arg1, double arg2 ) const
278 {
279  switch ( op )
280  {
281  case opPLUS:
282  return arg1 + arg2;
283  case opMINUS:
284  return arg1 - arg2;
285  case opMUL:
286  return arg1 * arg2;
287  case opDIV:
288  if ( arg2 == 0 )
289  {
290  return mNodataValue;
291  }
292  else
293  {
294  return arg1 / arg2;
295  }
296  case opPOW:
297  if ( !testPowerValidity( arg1, arg2 ) )
298  {
299  return mNodataValue;
300  }
301  else
302  {
303  return std::pow( arg1, arg2 );
304  }
305  case opEQ:
306  return ( arg1 == arg2 ? 1.0 : 0.0 );
307  case opNE:
308  return ( arg1 == arg2 ? 0.0 : 1.0 );
309  case opGT:
310  return ( arg1 > arg2 ? 1.0 : 0.0 );
311  case opLT:
312  return ( arg1 < arg2 ? 1.0 : 0.0 );
313  case opGE:
314  return ( arg1 >= arg2 ? 1.0 : 0.0 );
315  case opLE:
316  return ( arg1 <= arg2 ? 1.0 : 0.0 );
317  case opAND:
318  return ( arg1 && arg2 ? 1.0 : 0.0 );
319  case opOR:
320  return ( arg1 || arg2 ? 1.0 : 0.0 );
321  case opMAX:
322  return std::max( arg1, arg2 );
323  case opMIN:
324  return std::min( arg1, arg2 );
325  }
326  return mNodataValue;
327 }
328 
329 bool QgsRasterMatrix::twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix &other )
330 {
331  if ( isNumber() && other.isNumber() ) //operation on two 1x1 matrices
332  {
333  //operations with nodata values always generate nodata
334  if ( mData[0] == mNodataValue || other.number() == other.nodataValue() )
335  {
336  mData[0] = mNodataValue;
337  }
338  else
339  {
340  mData[0] = calculateTwoArgumentOp( op, mData[0], other.number() );
341  }
342  return true;
343  }
344 
345  //two matrices
346  if ( !isNumber() && !other.isNumber() )
347  {
348  double *matrix = other.mData;
349  int nEntries = mColumns * mRows;
350  double value1, value2;
351 
352  for ( int i = 0; i < nEntries; ++i )
353  {
354  value1 = mData[i];
355  value2 = matrix[i];
356  if ( value1 == mNodataValue || value2 == other.mNodataValue )
357  {
358  mData[i] = mNodataValue;
359  }
360  else
361  {
362  mData[i] = calculateTwoArgumentOp( op, value1, value2 );
363  }
364  }
365  return true;
366  }
367 
368  //this matrix is a single number and the other one a real matrix
369  if ( isNumber() )
370  {
371  double *matrix = other.mData;
372  int nEntries = other.nColumns() * other.nRows();
373  double value = mData[0];
374  delete[] mData;
375  mData = new double[nEntries];
376  mColumns = other.nColumns();
377  mRows = other.nRows();
378  mNodataValue = other.nodataValue();
379 
380  if ( value == mNodataValue )
381  {
382  for ( int i = 0; i < nEntries; ++i )
383  {
384  mData[i] = mNodataValue;
385  }
386  return true;
387  }
388 
389  for ( int i = 0; i < nEntries; ++i )
390  {
391  if ( matrix[i] == other.mNodataValue )
392  {
393  mData[i] = mNodataValue;
394  continue;
395  }
396 
397  mData[i] = calculateTwoArgumentOp( op, value, matrix[i] );
398  }
399  return true;
400  }
401  else //this matrix is a real matrix and the other a number
402  {
403  double value = other.number();
404  int nEntries = mColumns * mRows;
405 
406  if ( other.number() == other.mNodataValue )
407  {
408  for ( int i = 0; i < nEntries; ++i )
409  {
410  mData[i] = mNodataValue;
411  }
412  return true;
413  }
414 
415  for ( int i = 0; i < nEntries; ++i )
416  {
417  if ( mData[i] == mNodataValue )
418  {
419  continue;
420  }
421 
422  mData[i] = calculateTwoArgumentOp( op, mData[i], value );
423  }
424  return true;
425  }
426 }
427 
428 bool QgsRasterMatrix::testPowerValidity( double base, double power ) const
429 {
430  return !( ( base == 0 && power < 0 ) || ( base < 0 && ( power - std::floor( power ) ) > 0 ) );
431 }
double * takeData()
Returns data and ownership.
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
QgsRasterMatrix & operator=(const QgsRasterMatrix &m)
bool add(const QgsRasterMatrix &other)
Adds another matrix to this one.
double nodataValue() const
bool power(const QgsRasterMatrix &other)
bool absoluteValue()
Calculates the absolute value.
bool greaterThan(const QgsRasterMatrix &other)
bool notEqual(const QgsRasterMatrix &other)
bool equal(const QgsRasterMatrix &other)
int nRows() const
bool logicalOr(const QgsRasterMatrix &other)
void setData(int cols, int rows, double *data, double nodataValue)
QgsRasterMatrix()=default
Takes ownership of data array.
bool subtract(const QgsRasterMatrix &other)
Subtracts another matrix from this one.
bool multiply(const QgsRasterMatrix &other)
int nColumns() const
bool lesserEqual(const QgsRasterMatrix &other)
bool divide(const QgsRasterMatrix &other)
bool greaterEqual(const QgsRasterMatrix &other)
bool min(const QgsRasterMatrix &other)
Calculates the minimum value between two matrices.
bool max(const QgsRasterMatrix &other)
Calculates the maximum value between two matrices.
double number() const
double * data()
Returns data array (but not ownership)
bool lesserThan(const QgsRasterMatrix &other)
bool logicalAnd(const QgsRasterMatrix &other)