QGIS API Documentation  2.9.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 <string.h>
20 
21 #include <cmath>
22 
24  : mColumns( 0 )
25  , mRows( 0 )
26  , mData( 0 )
27  , mNodataValue( -1 )
28 {
29 }
30 
31 QgsRasterMatrix::QgsRasterMatrix( int nCols, int nRows, float* data, double nodataValue )
32  : mColumns( nCols )
33  , mRows( nRows )
34  , mData( data )
35  , mNodataValue( nodataValue )
36 {
37 }
38 
40  : mColumns( 0 )
41  , mRows( 0 )
42  , mData( 0 )
43 {
44  operator=( m );
45 }
46 
48 {
49  delete[] mData;
50 }
51 
53 {
54  delete[] mData;
55  mColumns = m.nColumns();
56  mRows = m.nRows();
57  int nEntries = mColumns * mRows;
58  mData = new float[nEntries];
59  memcpy( mData, m.mData, sizeof( float ) * nEntries );
60  mNodataValue = m.nodataValue();
61  return *this;
62 }
63 
64 void QgsRasterMatrix::setData( int cols, int rows, float* data, double nodataValue )
65 {
66  delete[] mData;
67  mColumns = cols;
68  mRows = rows;
69  mData = data;
70  mNodataValue = nodataValue;
71 }
72 
74 {
75  float* data = mData;
76  mData = 0; mColumns = 0; mRows = 0;
77  return data;
78 }
79 
81 {
82  return twoArgumentOperation( opPLUS, other );
83 }
84 
86 {
87  return twoArgumentOperation( opMINUS, other );
88 }
89 
91 {
92  return twoArgumentOperation( opMUL, other );
93 }
94 
96 {
97  return twoArgumentOperation( opDIV, other );
98 }
99 
101 {
102  return twoArgumentOperation( opPOW, other );
103 }
104 
106 {
107  return twoArgumentOperation( opEQ, other );
108 }
109 
111 {
112  return twoArgumentOperation( opNE, other );
113 }
114 
116 {
117  return twoArgumentOperation( opGT, other );
118 }
119 
121 {
122  return twoArgumentOperation( opLT, other );
123 }
124 
126 {
127  return twoArgumentOperation( opGE, other );
128 }
129 
131 {
132  return twoArgumentOperation( opLE, other );
133 }
134 
136 {
137  return twoArgumentOperation( opAND, other );
138 }
139 
141 {
142  return twoArgumentOperation( opOR, other );
143 }
144 
146 {
147  return oneArgumentOperation( opSQRT );
148 }
149 
151 {
152  return oneArgumentOperation( opSIN );
153 }
154 
156 {
157  return oneArgumentOperation( opASIN );
158 }
159 
161 {
162  return oneArgumentOperation( opCOS );
163 }
164 
166 {
167  return oneArgumentOperation( opACOS );
168 }
169 
171 {
172  return oneArgumentOperation( opTAN );
173 }
174 
176 {
177  return oneArgumentOperation( opATAN );
178 }
179 
181 {
182  return oneArgumentOperation( opSIGN );
183 }
184 
186 {
187  return oneArgumentOperation( opLOG );
188 }
189 
191 {
192  return oneArgumentOperation( opLOG10 );
193 }
194 
195 bool QgsRasterMatrix::oneArgumentOperation( OneArgOperator op )
196 {
197  if ( !mData )
198  {
199  return false;
200  }
201 
202  int nEntries = mColumns * mRows;
203  double value;
204  for ( int i = 0; i < nEntries; ++i )
205  {
206  value = mData[i];
207  if ( value != mNodataValue )
208  {
209  switch ( op )
210  {
211  case opSQRT:
212  if ( value < 0 ) //no complex numbers
213  {
214  mData[i] = static_cast<float>( mNodataValue );
215  }
216  else
217  {
218  mData[i] = static_cast<float>( sqrt( value ) );
219  }
220  break;
221  case opSIN:
222  mData[i] = static_cast<float>( sin( value ) );
223  break;
224  case opCOS:
225  mData[i] = static_cast<float>( cos( value ) );
226  break;
227  case opTAN:
228  mData[i] = static_cast<float>( tan( value ) );
229  break;
230  case opASIN:
231  mData[i] = static_cast<float>( asin( value ) );
232  break;
233  case opACOS:
234  mData[i] = static_cast<float>( acos( value ) );
235  break;
236  case opATAN:
237  mData[i] = static_cast<float>( atan( value ) );
238  break;
239  case opSIGN:
240  mData[i] = static_cast<float>( -value );
241  break;
242  case opLOG:
243  mData[i] = static_cast<float>( ::log( value ) );
244  break;
245  case opLOG10:
246  mData[i] = static_cast<float>( ::log10( value ) );
247  break;
248  }
249  }
250  }
251  return true;
252 }
253 
254 bool QgsRasterMatrix::twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix& other )
255 {
256  if ( isNumber() && other.isNumber() ) //operation on two 1x1 matrices
257  {
258  //operations with nodata values always generate nodata
259  if ( mData[0] == mNodataValue || other.number() == other.nodataValue() )
260  {
261  mData[0] = static_cast<float>( mNodataValue );
262  return true;
263  }
264  switch ( op )
265  {
266  case opPLUS:
267  mData[0] = static_cast<float>( number() + other.number() );
268  break;
269  case opMINUS:
270  mData[0] = static_cast<float>( number() - other.number() );
271  break;
272  case opMUL:
273  mData[0] = static_cast<float>( number() * other.number() );
274  break;
275  case opDIV:
276  if ( other.number() == 0 )
277  {
278  mData[0] = static_cast<float>( mNodataValue );
279  }
280  else
281  {
282  mData[0] = static_cast<float>( number() / other.number() );
283  }
284  break;
285  case opPOW:
286  if ( !testPowerValidity( mData[0], ( float ) other.number() ) )
287  {
288  mData[0] = static_cast<float>( mNodataValue );
289  }
290  else
291  {
292  mData[0] = pow( mData[0], ( float ) other.number() );
293  }
294  break;
295  case opEQ:
296  mData[0] = mData[0] == other.number() ? 1.0f : 0.0f;
297  break;
298  case opNE:
299  mData[0] = mData[0] == other.number() ? 0.0f : 1.0f;
300  break;
301  case opGT:
302  mData[0] = mData[0] > other.number() ? 1.0f : 0.0f;
303  break;
304  case opLT:
305  mData[0] = mData[0] < other.number() ? 1.0f : 0.0f;
306  break;
307  case opGE:
308  mData[0] = mData[0] >= other.number() ? 1.0f : 0.0f;
309  break;
310  case opLE:
311  mData[0] = mData[0] <= other.number() ? 1.0f : 0.0f;
312  break;
313  case opAND:
314  mData[0] = mData[0] && other.number() ? 1.0f : 0.0f;
315  break;
316  case opOR:
317  mData[0] = mData[0] || other.number() ? 1.0f : 0.0f;
318  break;
319  }
320  return true;
321  }
322 
323  //two matrices
324  if ( !isNumber() && !other.isNumber() )
325  {
326  float* matrix = other.mData;
327  int nEntries = mColumns * mRows;
328  double value1, value2;
329 
330  for ( int i = 0; i < nEntries; ++i )
331  {
332  value1 = mData[i]; value2 = matrix[i];
333  if ( value1 == mNodataValue || value2 == other.mNodataValue )
334  {
335  mData[i] = static_cast<float>( mNodataValue );
336  }
337  else
338  {
339  switch ( op )
340  {
341  case opPLUS:
342  mData[i] = static_cast<float>( value1 + value2 );
343  break;
344  case opMINUS:
345  mData[i] = static_cast<float>( value1 - value2 );
346  break;
347  case opMUL:
348  mData[i] = static_cast<float>( value1 * value2 );
349  break;
350  case opDIV:
351  if ( value2 == 0 )
352  {
353  mData[i] = static_cast<float>( mNodataValue );
354  }
355  else
356  {
357  mData[i] = static_cast<float>( value1 / value2 );
358  }
359  break;
360  case opPOW:
361  if ( !testPowerValidity( value1, value2 ) )
362  {
363  mData[i] = static_cast<float>( mNodataValue );
364  }
365  else
366  {
367  mData[i] = static_cast<float>( pow( value1, value2 ) );
368  }
369  break;
370  case opEQ:
371  mData[i] = value1 == value2 ? 1.0f : 0.0f;
372  break;
373  case opNE:
374  mData[i] = value1 == value2 ? 0.0f : 1.0f;
375  break;
376  case opGT:
377  mData[i] = value1 > value2 ? 1.0f : 0.0f;
378  break;
379  case opLT:
380  mData[i] = value1 < value2 ? 1.0f : 0.0f;
381  break;
382  case opGE:
383  mData[i] = value1 >= value2 ? 1.0f : 0.0f;
384  break;
385  case opLE:
386  mData[i] = value1 <= value2 ? 1.0f : 0.0f;
387  break;
388  case opAND:
389  mData[i] = value1 && value2 ? 1.0f : 0.0f;
390  break;
391  case opOR:
392  mData[i] = value1 || value2 ? 1.0f : 0.0f;
393  break;
394  }
395  }
396  }
397  return true;
398  }
399 
400  //this matrix is a single number and the other one a real matrix
401  if ( isNumber() )
402  {
403  float* matrix = other.mData;
404  int nEntries = other.nColumns() * other.nRows();
405  double value = mData[0];
406  delete[] mData;
407  mData = new float[nEntries]; mColumns = other.nColumns(); mRows = other.nRows();
408  mNodataValue = other.nodataValue();
409 
410  if ( value == mNodataValue )
411  {
412  for ( int i = 0; i < nEntries; ++i )
413  {
414  mData[i] = static_cast<float>( mNodataValue );
415  }
416  return true;
417  }
418 
419  for ( int i = 0; i < nEntries; ++i )
420  {
421  if ( matrix[i] == other.mNodataValue )
422  {
423  mData[i] = static_cast<float>( mNodataValue );
424  continue;
425  }
426 
427  switch ( op )
428  {
429  case opPLUS:
430  mData[i] = static_cast<float>( value + matrix[i] );
431  break;
432  case opMINUS:
433  mData[i] = static_cast<float>( value - matrix[i] );
434  break;
435  case opMUL:
436  mData[i] = static_cast<float>( value * matrix[i] );
437  break;
438  case opDIV:
439  if ( matrix[i] == 0 )
440  {
441  mData[i] = static_cast<float>( mNodataValue );
442  }
443  else
444  {
445  mData[i] = static_cast<float>( value / matrix[i] );
446  }
447  break;
448  case opPOW:
449  if ( !testPowerValidity( value, matrix[i] ) )
450  {
451  mData[i] = static_cast<float>( mNodataValue );
452  }
453  else
454  {
455  mData[i] = pow( static_cast<float>( value ), matrix[i] );
456  }
457  break;
458  case opEQ:
459  mData[i] = value == matrix[i] ? 1.0f : 0.0f;
460  break;
461  case opNE:
462  mData[i] = value == matrix[i] ? 0.0f : 1.0f;
463  break;
464  case opGT:
465  mData[i] = value > matrix[i] ? 1.0f : 0.0f;
466  break;
467  case opLT:
468  mData[i] = value < matrix[i] ? 1.0f : 0.0f;
469  break;
470  case opGE:
471  mData[i] = value >= matrix[i] ? 1.0f : 0.0f;
472  break;
473  case opLE:
474  mData[i] = value <= matrix[i] ? 1.0f : 0.0f;
475  break;
476  case opAND:
477  mData[i] = value && matrix[i] ? 1.0f : 0.0f;
478  break;
479  case opOR:
480  mData[i] = value || matrix[i] ? 1.0f : 0.0f;
481  break;
482  }
483  }
484  return true;
485  }
486  else //this matrix is a real matrix and the other a number
487  {
488  double value = other.number();
489  int nEntries = mColumns * mRows;
490 
491  if ( other.number() == other.mNodataValue )
492  {
493  for ( int i = 0; i < nEntries; ++i )
494  {
495  mData[i] = static_cast<float>( mNodataValue );
496  }
497  return true;
498  }
499 
500  for ( int i = 0; i < nEntries; ++i )
501  {
502  if ( mData[i] == mNodataValue )
503  {
504  continue;
505  }
506 
507  switch ( op )
508  {
509  case opPLUS:
510  mData[i] = static_cast<float>( mData[i] + value );
511  break;
512  case opMINUS:
513  mData[i] = static_cast<float>( mData[i] - value );
514  break;
515  case opMUL:
516  mData[i] = static_cast<float>( mData[i] * value );
517  break;
518  case opDIV:
519  if ( value == 0 )
520  {
521  mData[i] = static_cast<float>( mNodataValue );
522  }
523  else
524  {
525  mData[i] = static_cast<float>( mData[i] / value );
526  }
527  break;
528  case opPOW:
529  if ( !testPowerValidity( mData[i], value ) )
530  {
531  mData[i] = static_cast<float>( mNodataValue );
532  }
533  else
534  {
535  mData[i] = pow( mData[i], ( float ) value );
536  }
537  break;
538  case opEQ:
539  mData[i] = mData[i] == value ? 1.0f : 0.0f;
540  break;
541  case opNE:
542  mData[i] = mData[i] == value ? 0.0f : 1.0f;
543  break;
544  case opGT:
545  mData[i] = mData[i] > value ? 1.0f : 0.0f;
546  break;
547  case opLT:
548  mData[i] = mData[i] < value ? 1.0f : 0.0f;
549  break;
550  case opGE:
551  mData[i] = mData[i] >= value ? 1.0f : 0.0f;
552  break;
553  case opLE:
554  mData[i] = mData[i] <= value ? 1.0f : 0.0f;
555  break;
556  case opAND:
557  mData[i] = mData[i] && value ? 1.0f : 0.0f;
558  break;
559  case opOR:
560  mData[i] = mData[i] || value ? 1.0f : 0.0f;
561  break;
562  }
563  }
564  return true;
565  }
566 }
567 
568 bool QgsRasterMatrix::testPowerValidity( double base, double power )
569 {
570  if (( base == 0 && power < 0 ) || ( base < 0 && ( power - floor( power ) ) > 0 ) )
571  {
572  return false;
573  }
574  return true;
575 }
QgsRasterMatrix & operator=(const QgsRasterMatrix &m)
void setData(int cols, int rows, float *data, double nodataValue)
float * takeData()
Returns data and ownership.
bool add(const QgsRasterMatrix &other)
Adds another matrix to this one.
bool power(const QgsRasterMatrix &other)
QgsRasterMatrix()
Takes ownership of data array.
bool greaterThan(const QgsRasterMatrix &other)
int nColumns() const
bool notEqual(const QgsRasterMatrix &other)
bool equal(const QgsRasterMatrix &other)
double number() const
bool logicalOr(const QgsRasterMatrix &other)
bool subtract(const QgsRasterMatrix &other)
Subtracts another matrix from this one.
int nRows() const
bool multiply(const QgsRasterMatrix &other)
bool isNumber() const
Returns true if matrix is 1x1 (=scalar number)
bool lesserEqual(const QgsRasterMatrix &other)
bool divide(const QgsRasterMatrix &other)
double ANALYSIS_EXPORT power(double a, int b)
power function for integer coefficients
double nodataValue() const
bool greaterEqual(const QgsRasterMatrix &other)
bool lesserThan(const QgsRasterMatrix &other)
float * data()
Returns data array (but not ownership)
bool logicalAnd(const QgsRasterMatrix &other)