Branch data Line data Source code
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 : 0 : QgsRasterMatrix::QgsRasterMatrix( int nCols, int nRows, double *data, double nodataValue )
24 : 0 : : mColumns( nCols )
25 : 0 : , mRows( nRows )
26 : 0 : , mData( data )
27 : 0 : , mNodataValue( nodataValue )
28 : : {
29 : 0 : }
30 : :
31 : 0 : QgsRasterMatrix::QgsRasterMatrix( const QgsRasterMatrix &m )
32 : : {
33 : 0 : operator=( m );
34 : 0 : }
35 : :
36 : 0 : QgsRasterMatrix::~QgsRasterMatrix()
37 : : {
38 : 0 : delete[] mData;
39 : 0 : }
40 : :
41 : 0 : QgsRasterMatrix &QgsRasterMatrix::operator=( const QgsRasterMatrix &m )
42 : : {
43 : 0 : if ( this != &m )
44 : : {
45 : 0 : delete[] mData;
46 : 0 : mColumns = m.nColumns();
47 : 0 : mRows = m.nRows();
48 : 0 : int nEntries = mColumns * mRows;
49 : 0 : mData = new double[nEntries];
50 : 0 : memcpy( mData, m.mData, sizeof( double ) * nEntries );
51 : 0 : mNodataValue = m.nodataValue();
52 : 0 : }
53 : 0 : return *this;
54 : : }
55 : :
56 : 0 : void QgsRasterMatrix::setData( int cols, int rows, double *data, double nodataValue )
57 : : {
58 : 0 : delete[] mData;
59 : 0 : mColumns = cols;
60 : 0 : mRows = rows;
61 : 0 : mData = data;
62 : 0 : mNodataValue = nodataValue;
63 : 0 : }
64 : :
65 : 0 : double *QgsRasterMatrix::takeData()
66 : : {
67 : 0 : double *data = mData;
68 : 0 : mData = nullptr;
69 : 0 : mColumns = 0;
70 : 0 : mRows = 0;
71 : 0 : return data;
72 : : }
73 : :
74 : 0 : bool QgsRasterMatrix::add( const QgsRasterMatrix &other )
75 : : {
76 : 0 : return twoArgumentOperation( opPLUS, other );
77 : : }
78 : :
79 : 0 : bool QgsRasterMatrix::subtract( const QgsRasterMatrix &other )
80 : : {
81 : 0 : return twoArgumentOperation( opMINUS, other );
82 : : }
83 : :
84 : 0 : bool QgsRasterMatrix::multiply( const QgsRasterMatrix &other )
85 : : {
86 : 0 : return twoArgumentOperation( opMUL, other );
87 : : }
88 : :
89 : 0 : bool QgsRasterMatrix::divide( const QgsRasterMatrix &other )
90 : : {
91 : 0 : return twoArgumentOperation( opDIV, other );
92 : : }
93 : :
94 : 0 : bool QgsRasterMatrix::power( const QgsRasterMatrix &other )
95 : : {
96 : 0 : return twoArgumentOperation( opPOW, other );
97 : : }
98 : :
99 : 0 : bool QgsRasterMatrix::equal( const QgsRasterMatrix &other )
100 : : {
101 : 0 : return twoArgumentOperation( opEQ, other );
102 : : }
103 : :
104 : 0 : bool QgsRasterMatrix::notEqual( const QgsRasterMatrix &other )
105 : : {
106 : 0 : return twoArgumentOperation( opNE, other );
107 : : }
108 : :
109 : 0 : bool QgsRasterMatrix::greaterThan( const QgsRasterMatrix &other )
110 : : {
111 : 0 : return twoArgumentOperation( opGT, other );
112 : : }
113 : :
114 : 0 : bool QgsRasterMatrix::lesserThan( const QgsRasterMatrix &other )
115 : : {
116 : 0 : return twoArgumentOperation( opLT, other );
117 : : }
118 : :
119 : 0 : bool QgsRasterMatrix::greaterEqual( const QgsRasterMatrix &other )
120 : : {
121 : 0 : return twoArgumentOperation( opGE, other );
122 : : }
123 : :
124 : 0 : bool QgsRasterMatrix::lesserEqual( const QgsRasterMatrix &other )
125 : : {
126 : 0 : return twoArgumentOperation( opLE, other );
127 : : }
128 : :
129 : 0 : bool QgsRasterMatrix::logicalAnd( const QgsRasterMatrix &other )
130 : : {
131 : 0 : return twoArgumentOperation( opAND, other );
132 : : }
133 : :
134 : 0 : bool QgsRasterMatrix::logicalOr( const QgsRasterMatrix &other )
135 : : {
136 : 0 : return twoArgumentOperation( opOR, other );
137 : : }
138 : :
139 : 0 : bool QgsRasterMatrix::max( const QgsRasterMatrix &other )
140 : : {
141 : 0 : return twoArgumentOperation( opMAX, other );
142 : : }
143 : :
144 : 0 : bool QgsRasterMatrix::min( const QgsRasterMatrix &other )
145 : : {
146 : 0 : return twoArgumentOperation( opMIN, other );
147 : : }
148 : :
149 : 0 : bool QgsRasterMatrix::squareRoot()
150 : 0 : {
151 : 0 : return oneArgumentOperation( opSQRT );
152 : 0 : }
153 : :
154 : 0 : bool QgsRasterMatrix::sinus()
155 : : {
156 : 0 : return oneArgumentOperation( opSIN );
157 : : }
158 : :
159 : 0 : bool QgsRasterMatrix::asinus()
160 : : {
161 : 0 : return oneArgumentOperation( opASIN );
162 : : }
163 : :
164 : 0 : bool QgsRasterMatrix::cosinus()
165 : : {
166 : 0 : return oneArgumentOperation( opCOS );
167 : : }
168 : :
169 : 0 : bool QgsRasterMatrix::acosinus()
170 : : {
171 : 0 : return oneArgumentOperation( opACOS );
172 : : }
173 : :
174 : 0 : bool QgsRasterMatrix::tangens()
175 : : {
176 : 0 : return oneArgumentOperation( opTAN );
177 : : }
178 : :
179 : 0 : bool QgsRasterMatrix::atangens()
180 : : {
181 : 0 : return oneArgumentOperation( opATAN );
182 : : }
183 : :
184 : 0 : bool QgsRasterMatrix::changeSign()
185 : : {
186 : 0 : return oneArgumentOperation( opSIGN );
187 : : }
188 : :
189 : 0 : bool QgsRasterMatrix::log()
190 : : {
191 : 0 : return oneArgumentOperation( opLOG );
192 : : }
193 : :
194 : 0 : bool QgsRasterMatrix::log10()
195 : : {
196 : 0 : return oneArgumentOperation( opLOG10 );
197 : : }
198 : :
199 : 0 : bool QgsRasterMatrix::absoluteValue()
200 : : {
201 : 0 : return oneArgumentOperation( opABS );
202 : : }
203 : :
204 : 0 : bool QgsRasterMatrix::oneArgumentOperation( OneArgOperator op )
205 : : {
206 : 0 : if ( !mData )
207 : : {
208 : 0 : return false;
209 : : }
210 : :
211 : 0 : int nEntries = mColumns * mRows;
212 : : double value;
213 : 0 : for ( int i = 0; i < nEntries; ++i )
214 : : {
215 : 0 : value = mData[i];
216 : 0 : if ( value != mNodataValue )
217 : : {
218 : 0 : switch ( op )
219 : : {
220 : : case opSQRT:
221 : 0 : if ( value < 0 ) //no complex numbers
222 : : {
223 : 0 : mData[i] = mNodataValue;
224 : 0 : }
225 : : else
226 : : {
227 : 0 : mData[i] = std::sqrt( value );
228 : : }
229 : 0 : break;
230 : : case opSIN:
231 : 0 : mData[i] = std::sin( value );
232 : 0 : break;
233 : : case opCOS:
234 : 0 : mData[i] = std::cos( value );
235 : 0 : break;
236 : : case opTAN:
237 : 0 : mData[i] = std::tan( value );
238 : 0 : break;
239 : : case opASIN:
240 : 0 : mData[i] = std::asin( value );
241 : 0 : break;
242 : : case opACOS:
243 : 0 : mData[i] = std::acos( value );
244 : 0 : break;
245 : : case opATAN:
246 : 0 : mData[i] = std::atan( value );
247 : 0 : break;
248 : : case opSIGN:
249 : 0 : mData[i] = -value;
250 : 0 : break;
251 : : case opLOG:
252 : 0 : if ( value <= 0 )
253 : : {
254 : 0 : mData[i] = mNodataValue;
255 : 0 : }
256 : : else
257 : : {
258 : 0 : mData[i] = ::log( value );
259 : : }
260 : 0 : break;
261 : : case opLOG10:
262 : 0 : if ( value <= 0 )
263 : : {
264 : 0 : mData[i] = mNodataValue;
265 : 0 : }
266 : : else
267 : : {
268 : 0 : mData[i] = ::log10( value );
269 : : }
270 : 0 : break;
271 : : case opABS:
272 : 0 : mData[i] = ::fabs( value );
273 : 0 : break;
274 : : }
275 : 0 : }
276 : 0 : }
277 : 0 : return true;
278 : 0 : }
279 : :
280 : 0 : double QgsRasterMatrix::calculateTwoArgumentOp( TwoArgOperator op, double arg1, double arg2 ) const
281 : : {
282 : 0 : switch ( op )
283 : : {
284 : : case opPLUS:
285 : 0 : return arg1 + arg2;
286 : : case opMINUS:
287 : 0 : return arg1 - arg2;
288 : : case opMUL:
289 : 0 : return arg1 * arg2;
290 : : case opDIV:
291 : 0 : if ( arg2 == 0 )
292 : : {
293 : 0 : return mNodataValue;
294 : : }
295 : : else
296 : : {
297 : 0 : return arg1 / arg2;
298 : : }
299 : : case opPOW:
300 : 0 : if ( !testPowerValidity( arg1, arg2 ) )
301 : : {
302 : 0 : return mNodataValue;
303 : : }
304 : : else
305 : : {
306 : 0 : return std::pow( arg1, arg2 );
307 : : }
308 : : case opEQ:
309 : 0 : return ( arg1 == arg2 ? 1.0 : 0.0 );
310 : : case opNE:
311 : 0 : return ( arg1 == arg2 ? 0.0 : 1.0 );
312 : : case opGT:
313 : 0 : return ( arg1 > arg2 ? 1.0 : 0.0 );
314 : : case opLT:
315 : 0 : return ( arg1 < arg2 ? 1.0 : 0.0 );
316 : : case opGE:
317 : 0 : return ( arg1 >= arg2 ? 1.0 : 0.0 );
318 : : case opLE:
319 : 0 : return ( arg1 <= arg2 ? 1.0 : 0.0 );
320 : : case opAND:
321 : 0 : return ( arg1 && arg2 ? 1.0 : 0.0 );
322 : : case opOR:
323 : 0 : return ( arg1 || arg2 ? 1.0 : 0.0 );
324 : : case opMAX:
325 : 0 : return std::max( arg1, arg2 );
326 : : case opMIN:
327 : 0 : return std::min( arg1, arg2 );
328 : : }
329 : 0 : return mNodataValue;
330 : 0 : }
331 : :
332 : 0 : bool QgsRasterMatrix::twoArgumentOperation( TwoArgOperator op, const QgsRasterMatrix &other )
333 : : {
334 : 0 : if ( isNumber() && other.isNumber() ) //operation on two 1x1 matrices
335 : : {
336 : : //operations with nodata values always generate nodata
337 : 0 : if ( mData[0] == mNodataValue || other.number() == other.nodataValue() )
338 : : {
339 : 0 : mData[0] = mNodataValue;
340 : 0 : }
341 : : else
342 : : {
343 : 0 : mData[0] = calculateTwoArgumentOp( op, mData[0], other.number() );
344 : : }
345 : 0 : return true;
346 : : }
347 : :
348 : : //two matrices
349 : 0 : if ( !isNumber() && !other.isNumber() )
350 : : {
351 : 0 : double *matrix = other.mData;
352 : 0 : int nEntries = mColumns * mRows;
353 : : double value1, value2;
354 : :
355 : 0 : for ( int i = 0; i < nEntries; ++i )
356 : : {
357 : 0 : value1 = mData[i];
358 : 0 : value2 = matrix[i];
359 : 0 : if ( value1 == mNodataValue || value2 == other.mNodataValue )
360 : : {
361 : 0 : mData[i] = mNodataValue;
362 : 0 : }
363 : : else
364 : : {
365 : 0 : mData[i] = calculateTwoArgumentOp( op, value1, value2 );
366 : : }
367 : 0 : }
368 : 0 : return true;
369 : : }
370 : :
371 : : //this matrix is a single number and the other one a real matrix
372 : 0 : if ( isNumber() )
373 : : {
374 : 0 : double *matrix = other.mData;
375 : 0 : int nEntries = other.nColumns() * other.nRows();
376 : 0 : double value = mData[0];
377 : 0 : delete[] mData;
378 : 0 : mData = new double[nEntries];
379 : 0 : mColumns = other.nColumns();
380 : 0 : mRows = other.nRows();
381 : 0 : mNodataValue = other.nodataValue();
382 : :
383 : 0 : if ( value == mNodataValue )
384 : : {
385 : 0 : for ( int i = 0; i < nEntries; ++i )
386 : : {
387 : 0 : mData[i] = mNodataValue;
388 : 0 : }
389 : 0 : return true;
390 : : }
391 : :
392 : 0 : for ( int i = 0; i < nEntries; ++i )
393 : : {
394 : 0 : if ( matrix[i] == other.mNodataValue )
395 : : {
396 : 0 : mData[i] = mNodataValue;
397 : 0 : continue;
398 : : }
399 : :
400 : 0 : mData[i] = calculateTwoArgumentOp( op, value, matrix[i] );
401 : 0 : }
402 : 0 : return true;
403 : : }
404 : : else //this matrix is a real matrix and the other a number
405 : : {
406 : 0 : double value = other.number();
407 : 0 : int nEntries = mColumns * mRows;
408 : :
409 : 0 : if ( other.number() == other.mNodataValue )
410 : : {
411 : 0 : for ( int i = 0; i < nEntries; ++i )
412 : : {
413 : 0 : mData[i] = mNodataValue;
414 : 0 : }
415 : 0 : return true;
416 : : }
417 : :
418 : 0 : for ( int i = 0; i < nEntries; ++i )
419 : : {
420 : 0 : if ( mData[i] == mNodataValue )
421 : : {
422 : 0 : continue;
423 : : }
424 : :
425 : 0 : mData[i] = calculateTwoArgumentOp( op, mData[i], value );
426 : 0 : }
427 : 0 : return true;
428 : : }
429 : 0 : }
430 : :
431 : 0 : bool QgsRasterMatrix::testPowerValidity( double base, double power ) const
432 : : {
433 : 0 : return !( ( base == 0 && power < 0 ) || ( base < 0 && ( power - std::floor( power ) ) > 0 ) );
434 : : }
|
