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