QGIS API Documentation  2.8.2-Wien
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups 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 
185 bool QgsRasterMatrix::oneArgumentOperation( OneArgOperator op )
186 {
187  if ( !mData )
188  {
189  return false;
190  }
191 
192  int nEntries = mColumns * mRows;
193  double value;
194  for ( int i = 0; i < nEntries; ++i )
195  {
196  value = mData[i];
197  if ( value != mNodataValue )
198  {
199  switch ( op )
200  {
201  case opSQRT:
202  if ( value < 0 ) //no complex numbers
203  {
204  mData[i] = static_cast<float>( mNodataValue );
205  }
206  else
207  {
208  mData[i] = static_cast<float>( sqrt( value ) );
209  }
210  break;
211  case opSIN:
212  mData[i] = static_cast<float>( sin( value ) );
213  break;
214  case opCOS:
215  mData[i] = static_cast<float>( cos( value ) );
216  break;
217  case opTAN:
218  mData[i] = static_cast<float>( tan( value ) );
219  break;
220  case opASIN:
221  mData[i] = static_cast<float>( asin( value ) );
222  break;
223  case opACOS:
224  mData[i] = static_cast<float>( acos( value ) );
225  break;
226  case opATAN:
227  mData[i] = static_cast<float>( atan( value ) );
228  break;
229  case opSIGN:
230  mData[i] = static_cast<float>( -value );
231  }
232  }
233  }
234  return true;
235 }
236 
237 bool QgsRasterMatrix::twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix& other )
238 {
239  if ( isNumber() && other.isNumber() ) //operation on two 1x1 matrices
240  {
241  //operations with nodata values always generate nodata
242  if ( mData[0] == mNodataValue || other.number() == other.nodataValue() )
243  {
244  mData[0] = static_cast<float>( mNodataValue );
245  return true;
246  }
247  switch ( op )
248  {
249  case opPLUS:
250  mData[0] = static_cast<float>( number() + other.number() );
251  break;
252  case opMINUS:
253  mData[0] = static_cast<float>( number() - other.number() );
254  break;
255  case opMUL:
256  mData[0] = static_cast<float>( number() * other.number() );
257  break;
258  case opDIV:
259  if ( other.number() == 0 )
260  {
261  mData[0] = static_cast<float>( mNodataValue );
262  }
263  else
264  {
265  mData[0] = static_cast<float>( number() / other.number() );
266  }
267  break;
268  case opPOW:
269  if ( !testPowerValidity( mData[0], ( float ) other.number() ) )
270  {
271  mData[0] = static_cast<float>( mNodataValue );
272  }
273  else
274  {
275  mData[0] = pow( mData[0], ( float ) other.number() );
276  }
277  break;
278  case opEQ:
279  mData[0] = mData[0] == other.number() ? 1.0f : 0.0f;
280  break;
281  case opNE:
282  mData[0] = mData[0] == other.number() ? 0.0f : 1.0f;
283  break;
284  case opGT:
285  mData[0] = mData[0] > other.number() ? 1.0f : 0.0f;
286  break;
287  case opLT:
288  mData[0] = mData[0] < other.number() ? 1.0f : 0.0f;
289  break;
290  case opGE:
291  mData[0] = mData[0] >= other.number() ? 1.0f : 0.0f;
292  break;
293  case opLE:
294  mData[0] = mData[0] <= other.number() ? 1.0f : 0.0f;
295  break;
296  case opAND:
297  mData[0] = mData[0] && other.number() ? 1.0f : 0.0f;
298  break;
299  case opOR:
300  mData[0] = mData[0] || other.number() ? 1.0f : 0.0f;
301  break;
302  }
303  return true;
304  }
305 
306  //two matrices
307  if ( !isNumber() && !other.isNumber() )
308  {
309  float* matrix = other.mData;
310  int nEntries = mColumns * mRows;
311  double value1, value2;
312 
313  for ( int i = 0; i < nEntries; ++i )
314  {
315  value1 = mData[i]; value2 = matrix[i];
316  if ( value1 == mNodataValue || value2 == other.mNodataValue )
317  {
318  mData[i] = static_cast<float>( mNodataValue );
319  }
320  else
321  {
322  switch ( op )
323  {
324  case opPLUS:
325  mData[i] = static_cast<float>( value1 + value2 );
326  break;
327  case opMINUS:
328  mData[i] = static_cast<float>( value1 - value2 );
329  break;
330  case opMUL:
331  mData[i] = static_cast<float>( value1 * value2 );
332  break;
333  case opDIV:
334  if ( value2 == 0 )
335  {
336  mData[i] = static_cast<float>( mNodataValue );
337  }
338  else
339  {
340  mData[i] = static_cast<float>( value1 / value2 );
341  }
342  break;
343  case opPOW:
344  if ( !testPowerValidity( value1, value2 ) )
345  {
346  mData[i] = static_cast<float>( mNodataValue );
347  }
348  else
349  {
350  mData[i] = static_cast<float>( pow( value1, value2 ) );
351  }
352  break;
353  case opEQ:
354  mData[i] = value1 == value2 ? 1.0f : 0.0f;
355  break;
356  case opNE:
357  mData[i] = value1 == value2 ? 0.0f : 1.0f;
358  break;
359  case opGT:
360  mData[i] = value1 > value2 ? 1.0f : 0.0f;
361  break;
362  case opLT:
363  mData[i] = value1 < value2 ? 1.0f : 0.0f;
364  break;
365  case opGE:
366  mData[i] = value1 >= value2 ? 1.0f : 0.0f;
367  break;
368  case opLE:
369  mData[i] = value1 <= value2 ? 1.0f : 0.0f;
370  break;
371  case opAND:
372  mData[i] = value1 && value2 ? 1.0f : 0.0f;
373  break;
374  case opOR:
375  mData[i] = value1 || value2 ? 1.0f : 0.0f;
376  break;
377  }
378  }
379  }
380  return true;
381  }
382 
383  //this matrix is a single number and the other one a real matrix
384  if ( isNumber() )
385  {
386  float* matrix = other.mData;
387  int nEntries = other.nColumns() * other.nRows();
388  double value = mData[0];
389  delete[] mData;
390  mData = new float[nEntries]; mColumns = other.nColumns(); mRows = other.nRows();
391  mNodataValue = other.nodataValue();
392 
393  if ( value == mNodataValue )
394  {
395  for ( int i = 0; i < nEntries; ++i )
396  {
397  mData[i] = static_cast<float>( mNodataValue );
398  }
399  return true;
400  }
401 
402  for ( int i = 0; i < nEntries; ++i )
403  {
404  if ( matrix[i] == other.mNodataValue )
405  {
406  mData[i] = static_cast<float>( mNodataValue );
407  continue;
408  }
409 
410  switch ( op )
411  {
412  case opPLUS:
413  mData[i] = static_cast<float>( value + matrix[i] );
414  break;
415  case opMINUS:
416  mData[i] = static_cast<float>( value - matrix[i] );
417  break;
418  case opMUL:
419  mData[i] = static_cast<float>( value * matrix[i] );
420  break;
421  case opDIV:
422  if ( matrix[i] == 0 )
423  {
424  mData[i] = static_cast<float>( mNodataValue );
425  }
426  else
427  {
428  mData[i] = static_cast<float>( value / matrix[i] );
429  }
430  break;
431  case opPOW:
432  if ( !testPowerValidity( value, matrix[i] ) )
433  {
434  mData[i] = static_cast<float>( mNodataValue );
435  }
436  else
437  {
438  mData[i] = pow( static_cast<float>( value ), matrix[i] );
439  }
440  break;
441  case opEQ:
442  mData[i] = value == matrix[i] ? 1.0f : 0.0f;
443  break;
444  case opNE:
445  mData[i] = value == matrix[i] ? 0.0f : 1.0f;
446  break;
447  case opGT:
448  mData[i] = value > matrix[i] ? 1.0f : 0.0f;
449  break;
450  case opLT:
451  mData[i] = value < matrix[i] ? 1.0f : 0.0f;
452  break;
453  case opGE:
454  mData[i] = value >= matrix[i] ? 1.0f : 0.0f;
455  break;
456  case opLE:
457  mData[i] = value <= matrix[i] ? 1.0f : 0.0f;
458  break;
459  case opAND:
460  mData[i] = value && matrix[i] ? 1.0f : 0.0f;
461  break;
462  case opOR:
463  mData[i] = value || matrix[i] ? 1.0f : 0.0f;
464  break;
465  }
466  }
467  return true;
468  }
469  else //this matrix is a real matrix and the other a number
470  {
471  double value = other.number();
472  int nEntries = mColumns * mRows;
473 
474  if ( other.number() == other.mNodataValue )
475  {
476  for ( int i = 0; i < nEntries; ++i )
477  {
478  mData[i] = static_cast<float>( mNodataValue );
479  }
480  return true;
481  }
482 
483  for ( int i = 0; i < nEntries; ++i )
484  {
485  if ( mData[i] == mNodataValue )
486  {
487  continue;
488  }
489 
490  switch ( op )
491  {
492  case opPLUS:
493  mData[i] = static_cast<float>( mData[i] + value );
494  break;
495  case opMINUS:
496  mData[i] = static_cast<float>( mData[i] - value );
497  break;
498  case opMUL:
499  mData[i] = static_cast<float>( mData[i] * value );
500  break;
501  case opDIV:
502  if ( value == 0 )
503  {
504  mData[i] = static_cast<float>( mNodataValue );
505  }
506  else
507  {
508  mData[i] = static_cast<float>( mData[i] / value );
509  }
510  break;
511  case opPOW:
512  if ( !testPowerValidity( mData[i], value ) )
513  {
514  mData[i] = static_cast<float>( mNodataValue );
515  }
516  else
517  {
518  mData[i] = pow( mData[i], ( float ) value );
519  }
520  break;
521  case opEQ:
522  mData[i] = mData[i] == value ? 1.0f : 0.0f;
523  break;
524  case opNE:
525  mData[i] = mData[i] == value ? 0.0f : 1.0f;
526  break;
527  case opGT:
528  mData[i] = mData[i] > value ? 1.0f : 0.0f;
529  break;
530  case opLT:
531  mData[i] = mData[i] < value ? 1.0f : 0.0f;
532  break;
533  case opGE:
534  mData[i] = mData[i] >= value ? 1.0f : 0.0f;
535  break;
536  case opLE:
537  mData[i] = mData[i] <= value ? 1.0f : 0.0f;
538  break;
539  case opAND:
540  mData[i] = mData[i] && value ? 1.0f : 0.0f;
541  break;
542  case opOR:
543  mData[i] = mData[i] || value ? 1.0f : 0.0f;
544  break;
545  }
546  }
547  return true;
548  }
549 }
550 
551 bool QgsRasterMatrix::testPowerValidity( double base, double power )
552 {
553  if (( base == 0 && power < 0 ) || ( base < 0 && ( power - floor( power ) ) > 0 ) )
554  {
555  return false;
556  }
557  return true;
558 }