WSL/SLF GitLab Repository

Array3D.h 27 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef ARRAY3D_H
#define ARRAY3D_H

#include <meteoio/IOUtils.h>
#include <meteoio/IOExceptions.h>

#include <vector>
#include <limits>
#include <iostream>

namespace mio {

template <class T> class Array3D;
template <class T> class Array3DProxy2;

/**
 * @class Array3DProxy
 * @brief The template class Array3DProxy is a helper class for the template class Array3D
36
 *        with the purpose of adding the [][] operator to Array3D
37
38
39
40
41
42
 *
 * @author Thomas Egger
 */
template <class T> class Array3DProxy {
 	public:
		friend class Array3D<T>;
43
44
		Array3DProxy2<T> operator[](const unsigned int& i_any) {
			return Array3DProxy2<T>(array3D, anx, i_any);
45
46
47
		}

 	private:
48
 		Array3DProxy(Array3D<T>& i_array3D, const unsigned int& i_anx) : array3D(i_array3D), anx(i_anx){}
49
50
51
52
53
54
55
56
57
58
59
60
61
62
		Array3D<T>& array3D;
		const unsigned int anx;
};

/**
 * @class Array3DProxy2
 * @brief The template class Array3DProxy2 is a helper class for the template class Array3D
 *        with the purpose of adding the [][][] operator to Array3D
 *
 * @author Thomas Egger
 */
template <class T> class Array3DProxy2 {
 	public:
		friend class Array3DProxy<T>;
63
64
		T& operator[](const unsigned int& i_anz) {
			return array3D(anx, any, i_anz);
65
66
67
		}

	private:
68
69
 		Array3DProxy2(Array3D<T>& i_array3D, const unsigned int& i_anx,
				    const unsigned int& i_any) : array3D(i_array3D), anx(i_anx), any(i_any){}
70
71
72
		Array3D<T>& array3D;
		const unsigned int anx;
		const unsigned int any;
73
};
74
75
76

/**
 * @class Array3D
77
 * @brief The template class Array3D is a 3D Array (Tensor) able to hold any type of object as datatype.
78
 * It relies on the Array3DProxy2 class to provide the [][][] operator (slower than the (i,j,k) call).
79
 * If the compilation flag NOSAFECHECKS is used, bounds check is turned off (leading to increased performances).
80
 * @ingroup data_str
81
82
83
84
85
86
87
88
89
 * @date  2009-07-19
 * @author Thomas Egger
 */
template<class T> class Array3D {
	public:
		Array3D();

		/**
		* A constructor that can be used to create an Array3D object that is contained in the
90
91
92
93
94
95
96
97
98
		* one passed as i_array3D argument. The resulting Array3D object is a by value copy of
		* a subvolume of the volume spanned by the i_array3D
		* @param i_array3D array containing to extract the values from
		* @param i_nx lower left corner cell X index
		* @param i_ny lower left corner cell Y index
		* @param i_nz lower left corner cell Z index
		* @param i_ncols number of columns of the new array
		* @param i_nrows number of rows of the new array
		* @param i_ndepth number of depths of the new array
99
		*/
100
101
102
		Array3D(const Array3D<T>& i_array3D,
		        const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz,
		        const unsigned int& i_ncols, const unsigned int& i_nrows, const unsigned int& i_ndepth);
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

		/**
		* A constructor that creates an array of a given size
		* @param anx number of columns of the new array
		* @param any number of rows of the new array
		* @param anz number of rows of the new array
		*/
		Array3D(const unsigned int& anx, const unsigned int& any, const unsigned int& anz);

		/**
		* A constructor that creates an array filled with constant values
		* @param anx number of columns of the new array
		* @param any number of rows of the new array
		* @param anz number of depths of the new array
		* @param init initial value to fill the array with
		*/
		Array3D(const unsigned int& anx, const unsigned int& any, const unsigned int& anz, const T& init);

		/**
		* A method that can be used to create an Array3D object that is contained in the
123
124
125
126
127
128
129
130
131
		* one passed as i_array3D argument. The resulting Array3D object is a by value copy of
		* a subvolume of the volume spanned by the i_array3D
		* @param i_array3D array containing to extract the values from
		* @param i_nx lower left corner cell X index
		* @param i_ny lower left corner cell Y index
		* @param i_nz lower left corner cell Z index
		* @param i_ncols number of columns of the new array
		* @param i_nrows number of rows of the new array
		* @param i_ndepth number of depths of the new array
132
		*/
133
134
135
		void subset(const Array3D<T>& i_array3D,
		            const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz,
		            const unsigned int& i_ncols, const unsigned int& i_nrows, const unsigned int& i_ndepth);
136

137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
		/**
		* @brief A method that can be used to insert a subplane into an existing Array2D object
		* that is passed as i_array2D argument. This is exactly the opposite of the subset method
		* an can be used to rebuild an array from subsets.
		* @param i_array3D array containing to extract the values from
		* @param i_nx lower left corner cell X index
		* @param i_ny lower left corner cell Y index
		* @param i_nz lower left corner cell Z index
		* @param i_ncols number of columns of the new array
		* @param i_nrows number of rows of the new array
		* @param i_ndepth number of depths of the new array
		*/
		void fill(const Array3D<T>& i_array3D,
		          const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz,
		          const unsigned int& i_ncols, const unsigned int& i_nrows, const unsigned int& i_ndepth);

153
154
		void fill(const Array3D<T>& i_array3D, const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz);

155
156
		/**
		* @brief set how to process nodata values (ie: as nodata or as normal numbers)
157
158
		* @param i_keep_nodata true means that NODATA is interpreted as NODATA, false means that it is a normal number
		By default, arrays keep nodata.
159
		*/
160
161
162
163
164
165
166
		void setKeepNodata(const bool i_keep_nodata);

		/**
		* @brief get how to process nodata values (ie: as nodata or as normal numbers)
		* @return true means that NODATA is interpreted as NODATA, false means that it is a normal number
		*/
		bool getKeepNodata();
167

168
169
170
		void resize(const unsigned int& anx, const unsigned int& any, const unsigned int& anz);
		void resize(const unsigned int& anx, const unsigned int& any, const unsigned int& anz, const T& init);
		void size(unsigned int& anx, unsigned int& any, unsigned int& anz) const;
171
172
173
		unsigned int getNx() const;
		unsigned int getNy() const;
		unsigned int getNz() const;
174
		void clear();
175
		bool isEmpty() const;
176
177
178
179
180

		/**
		* @brief returns the minimum value contained in the grid
		* @return minimum value
		*/
181
		T getMin() const;
182
183
184
185
		/**
		* @brief returns the maximum value contained in the grid
		* @return maximum value
		*/
186
		T getMax() const;
187
188
189
190
		/**
		* @brief returns the mean value contained in the grid
		* @return mean value
		*/
191
		T getMean() const;
192
		/**
193
		* @brief returns the number of points contained in the grid.
194
195
		* If setNodataHandling(IOUtils::RAW_NODATA), then the number of points is the size of the grid.
		* If setNodataHandling(IOUtils::PARSE_NODATA), then it is the number of non-nodata values in the grid
196
		* @return count
197
		*/
198
		size_t getCount() const;
199
200
201
202
		/**
		* @brief returns the grid of the absolute value of values contained in the grid
		* @return grid of abs(grid)
		*/
203
204
		const Array3D<T> getAbs() const;
		void abs();
205

206
207
208
		const std::string toString() const;
		template<class P> friend std::iostream& operator<<(std::iostream& os, const Array3D<P>& array);
		template<class P> friend std::iostream& operator>>(std::iostream& is, Array3D<P>& array);
209

210
211
212
		bool checkEpsilonEquality(const Array3D<double>& rhs, const double& epsilon) const;
		static bool checkEpsilonEquality(const Array3D<double>& rhs1, const Array3D<double>& rhs2, const double& epsilon);

213
214
		T& operator ()(const unsigned int& i);
		const T operator ()(const unsigned int& i) const;
215
216
217
218
219
		T& operator ()(const unsigned int& x, const unsigned int& y, const unsigned int& z);
		const T operator ()(const unsigned int& x, const unsigned int& y, const unsigned int& z) const;
		Array3DProxy<T> operator[](const unsigned int& i);

		Array3D<T>& operator =(const Array3D<T>&);
220
		Array3D<T>& operator =(const T& value);
221

222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
		Array3D<T>& operator+=(const T& rhs);
		const Array3D<T> operator+(const T& rhs);
		Array3D<T>& operator+=(const Array3D<T>& rhs);
		const Array3D<T> operator+(const Array3D<T>& rhs);

		Array3D<T>& operator-=(const T& rhs);
		const Array3D<T> operator-(const T& rhs);
		Array3D<T>& operator-=(const Array3D<T>& rhs);
		const Array3D<T> operator-(const Array3D<T>& rhs);

		Array3D<T>& operator*=(const T& rhs);
		const Array3D<T> operator*(const T& rhs);
		Array3D<T>& operator*=(const Array3D<T>& rhs);
		const Array3D<T> operator*(const Array3D<T>& rhs);

		Array3D<T>& operator/=(const T& rhs);
		const Array3D<T> operator/(const T& rhs);
		Array3D<T>& operator/=(const Array3D<T>& rhs);
		const Array3D<T> operator/(const Array3D<T>& rhs);

242
243
244
		bool operator==(const Array3D<T>&) const; ///<Operator that tests for equality
		bool operator!=(const Array3D<T>&) const; ///<Operator that tests for inequality

245
246
247
248
249
250
	protected:
		std::vector<T> vecData; ///< The actual objects are stored in a one-dimensional vector
		unsigned int nx;
		unsigned int ny;
		unsigned int nz;
		unsigned int nxny; //nx times ny
251
		bool keep_nodata;
252
253
};

254
template<class T> inline T& Array3D<T>::operator()(const unsigned int& i) {
255
256
257
#ifndef NOSAFECHECKS
	return vecData.at(i);
#else
258
	return vecData[i];
259
#endif
260
261
}

262
template<class T> inline const T Array3D<T>::operator()(const unsigned int& i) const {
263
264
265
#ifndef NOSAFECHECKS
	return vecData.at(i);
#else
266
	return vecData[i];
267
#endif
268
269
}

270
template<class T> inline T& Array3D<T>::operator()(const unsigned int& x, const unsigned int& y, const unsigned int& z) {
271
#ifndef NOSAFECHECKS
272
273
274
	if ((x >= nx) || (y >= ny) || (z >= nz))  {
		std::stringstream ss;
		ss << "Trying to access array(" << x << "," << y << "," << z << ")";
275
		ss << " while array is (" << nx << "," << ny << "," << nz << ")";
276
		throw IndexOutOfBoundsException(ss.str(), AT);
277
278
279
280
281
282
283
	}
#endif

	//ROW-MAJOR alignment of the vector: fully C-compatible memory layout
	return vecData[x + y*nx + z*nxny];
}

284
template<class T> inline const T Array3D<T>::operator()(const unsigned int& x, const unsigned int& y, const unsigned int& z) const {
285
#ifndef NOSAFECHECKS
286
287
288
	if ((x >= nx) || (y >= ny) || (z >= nz))  {
		std::stringstream ss;
		ss << "Trying to access array(" << x << "," << y << "," << z << ")";
289
		ss << " while array is (" << nx << "," << ny << "," << nz << ")";
290
		throw IndexOutOfBoundsException(ss.str(), AT);
291
292
293
294
295
296
	}
#endif
	return vecData[x + y*nx + z*nxny];
}

template<class T> Array3DProxy<T> Array3D<T>::operator[](const unsigned int& i) {
297
	return Array3DProxy<T>(*this, i);
298
299
}

300
301
template<class T> Array3D<T>::Array3D() : vecData(), nx(0), ny(0), nz(0), nxny(0), keep_nodata(true)
{
302
303
}

304
305
306
template<class T> Array3D<T>::Array3D(const Array3D<T>& i_array3D,
                                      const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz,
                                      const unsigned int& i_ncols, const unsigned int& i_nrows, const unsigned int& i_ndepth)
307
                               : vecData(i_ncols*i_nrows*i_ndepth), nx(i_ncols), ny(i_nrows), nz(i_ndepth), nxny(i_ncols*i_nrows), keep_nodata(true)
308
{
309
	subset(i_array3D, i_nx, i_ny, i_nz, i_ncols, i_nrows, i_ndepth);
310
311
}

312
313
314
template<class T> void Array3D<T>::subset(const Array3D<T>& i_array3D,
                                     const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz,
                                     const unsigned int& i_ncols, const unsigned int& i_nrows, const unsigned int& i_ndepth)
315
{
316

317
318
319
320
321
322
323
	if (((i_nx+i_ncols) > i_array3D.nx) || ((i_ny+i_nrows) > i_array3D.ny) || ((i_nz+i_ndepth) > i_array3D.nz)) {
		std::stringstream ss;
		ss << "Trying to cut an array of size (" << nx << "," << ny << "," << nz << ") ";
		ss << "to size (" << i_ncols << "," << i_nrows << "," << i_ndepth << ") ";
		ss << "starting at (" << i_nx << "," << i_ny << "," << i_nz << ")";
		throw IndexOutOfBoundsException(ss.str(), AT);
	}
324

325
	if ((i_ncols == 0) || (i_nrows == 0) || (i_ndepth == 0)) //the space has to make sense
326
		throw IndexOutOfBoundsException("Trying to cut an array into a null sized array!", AT);
327

328
	resize(i_ncols, i_nrows, i_ndepth); //create new Array3D object
329
	//Copy by value subspace
330
	for (unsigned int ii=0; ii<nz; ii++) {
331
332
333
		for (unsigned int jj=0; jj<ny; jj++) {
			for (unsigned int kk=0; kk<nx; kk++) {
				//Running through the vector in order of memory alignment
334
				operator()(kk,jj,ii) = i_array3D(i_nx+kk, i_ny+jj, i_nz+ii);
335
336
337
338
339
			}
		}
	}
}

340
341
342
343
344
345
346
template<class T> void Array3D<T>::fill(const Array3D<T>& i_array3D, const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz)
{
	unsigned int i_ncols, i_nrows, i_ndepth;
	i_array3D.size(i_ncols, i_nrows, i_ndepth);
	fill(i_array3D, i_nx, i_ny, i_nz, i_ncols, i_nrows, i_ndepth);
}

347
348
349
350
351
template<class T> void Array3D<T>::fill(const Array3D<T>& i_array3D,
                                     const unsigned int& i_nx, const unsigned int& i_ny, const unsigned int& i_nz,
                                     const unsigned int& i_ncols, const unsigned int& i_nrows, const unsigned int& i_ndepth)
{

352
353
354
355
356
357
358
	if (((i_nx+i_ncols) > i_array3D.nx) || ((i_ny+i_nrows) > i_array3D.ny) || ((i_nz+i_ndepth) > i_array3D.nz)) {
		std::stringstream ss;
		ss << "Filling an array of size (" << nx << "," << ny << "," << nz << ") ";
		ss << "with an array of size (" << i_ncols << "," << i_nrows << "," << i_ndepth << ") ";
		ss << "starting at (" << i_nx << "," << i_ny << "," << i_nz << ")";
		throw IndexOutOfBoundsException(ss.str(), AT);
	}
359
360

	if ((i_ncols == 0) || (i_nrows == 0) || (i_ndepth == 0)) //the space has to make sense
361
		throw IndexOutOfBoundsException("Filling an array with a null sized array!", AT);
362
363
364
365
366
367
368
369
370
371
372
373
374
375

	//Copy by value subspace
	for (unsigned int ii=i_nz; ii<(i_nz+i_ndepth); ii++) {
		for (unsigned int jj=i_ny; jj<(i_ny+i_nrows); jj++) {
			for (unsigned int kk=i_nx; kk<(i_nx+i_ncols); kk++) {
				const unsigned int ix = kk-i_nx;
				const unsigned int iy = jj-i_ny;
				const unsigned int iz = ii-i_nz;
				operator()(kk,jj,ii) = i_array3D(ix, iy, iz);
			}
		}
	}
}

376
377
378
379
template<class T> Array3D<T>::Array3D(const unsigned int& anx, const unsigned int& any, const unsigned int& anz)
                             : vecData(anx*any*anz), nx(anx), ny(any), nz(anz), nxny(anx*any), keep_nodata(true)
{
	//resize(anx, any, anz);
380
381
}

382
383
384
385
template<class T> Array3D<T>::Array3D(const unsigned int& anx, const unsigned int& any, const unsigned int& anz, const T& init)
                             : vecData(anx*any*anz, init), nx(anx), ny(any), nz(anz), nxny(anx*any), keep_nodata(true)
{
	//resize(anx, any, anz, init);
386
387
}

388
389
390
391
392
393
template<class T> void Array3D<T>::setKeepNodata(const bool i_keep_nodata) {
	keep_nodata = i_keep_nodata;
}

template<class T> bool Array3D<T>::getKeepNodata() {
	return keep_nodata;
394
395
}

396
template<class T> void Array3D<T>::resize(const unsigned int& anx, const unsigned int& any, const unsigned int& anz) {
397
398
399
400
401
402
	clear();  //we won't be able to "rescue" old values, so we reset the whole vector
	vecData.resize(anx*any*anz);
	nx = anx;
	ny = any;
	nz = anz;
	nxny = nx*ny;
403
404
}

405
template<class T> void Array3D<T>::resize(const unsigned int& anx, const unsigned int& any, const unsigned int& anz, const T& init) {
406
	clear();  //we won't be able to "rescue" old values, so we reset the whole vector
407
408
409
410
411
	vecData.resize(anx*any*anz, init);
	nx = anx;
	ny = any;
	nz = anz;
	nxny = nx*ny;
412
413
}

414
415
416
417
template<class T> void Array3D<T>::size(unsigned int& anx, unsigned int& any, unsigned int& anz) const {
	anx=nx;
	any=ny;
	anz=nz;
418
419
}

420
421
422
423
424
425
426
427
428
429
430
431
template<class T> unsigned int Array3D<T>::getNx() const {
	return nx;
}

template<class T> unsigned int Array3D<T>::getNy() const {
	return ny;
}

template<class T> unsigned int Array3D<T>::getNz() const {
	return nz;
}

432
433
434
435
436
template<class T> void Array3D<T>::clear() {
	vecData.clear();
	nx = ny = nz = nxny = 0;
}

437
438
439
440
template<class T> bool Array3D<T>::isEmpty() const {
	return (nx==0 && ny==0 && nz==0);
}

441
template<class T> const std::string Array3D<T>::toString() const {
442
	std::ostringstream os;
443
	os << "<array3d>\n";
444
	for (unsigned int kk=0; kk<nz; kk++) {
445
		os << "depth[" << kk << "]\n";
446
447
448
		for(unsigned int ii=0; ii<nx; ii++) {
			for (unsigned int jj=0; jj<ny; jj++) {
				os << operator()(ii,jj,kk) << " ";
449
450
451
452
453
			}
			os << "\n";
		}
	}
	os << "</array3d>\n";
454
455
456
457
458
459
460
461
462
	return os.str();
}

template<class P> std::iostream& operator<<(std::iostream& os, const Array3D<P>& array) {
	os.write(reinterpret_cast<const char*>(&array.keep_nodata), sizeof(array.keep_nodata));
	os.write(reinterpret_cast<const char*>(&array.nx), sizeof(array.nx));
	os.write(reinterpret_cast<const char*>(&array.ny), sizeof(array.ny));
	os.write(reinterpret_cast<const char*>(&array.nz), sizeof(array.nz));
	os.write(reinterpret_cast<const char*>(&array.vecData[0]), array.nx*array.ny*array.nz*sizeof(P));
463
464
465
	return os;
}

466
467
468
469
470
471
472
template<class P> std::iostream& operator>>(std::iostream& is, Array3D<P>& array) {
	is.read(reinterpret_cast<char*>(&array.keep_nodata), sizeof(array.keep_nodata));
	is.read(reinterpret_cast<char*>(&array.nx), sizeof(array.nx));
	is.read(reinterpret_cast<char*>(&array.ny), sizeof(array.ny));
	is.read(reinterpret_cast<char*>(&array.nz), sizeof(array.nz));
	array.vecData.resize(array.nx*array.ny*array.nz);
	is.read(reinterpret_cast<char*>(&array.vecData[0]), array.nx*array.ny*array.nz*sizeof(P)); //30 times faster than assign() or copy()
473
	return is;
474
475
}

476
template<class T> T Array3D<T>::getMin() const {
477
478

	T min = std::numeric_limits<T>::max();
479
	const unsigned int nxyz = ny*nx*nz;
480

481
	if(keep_nodata==false) {
482
		for (unsigned int jj=0; jj<nxyz; jj++) {
483
			const T val = vecData[jj];
484
			if(val<min) min=val;
485
486
		}
		return min;
487
	} else {
488
		for (unsigned int jj=0; jj<nxyz; jj++) {
489
			const T val = vecData[jj];
490
			if(val!=IOUtils::nodata && val<min) min=val;
491
492
493
494
495
496
		}
		if(min!=std::numeric_limits<T>::max()) return min;
		else return (T)IOUtils::nodata;
	}
}

497
template<class T> T Array3D<T>::getMax() const {
498
499

	T max = -std::numeric_limits<T>::max();
500
	const unsigned int nxyz = ny*nx*nz;
501

502
	if(keep_nodata==false) {
503
		for (unsigned int jj=0; jj<nxyz; jj++) {
504
			const T val = vecData[jj];
505
			if(val>max) max=val;
506
507
		}
		return max;
508
	} else {
509
		for (unsigned int jj=0; jj<nxyz; jj++) {
510
			const T val = vecData[jj];
511
			if(val!=IOUtils::nodata && val>max) max=val;
512
513
514
515
516
517
		}
		if(max!=-std::numeric_limits<T>::max()) return max;
		else return (T)IOUtils::nodata;
	}
}

518
template<class T> T Array3D<T>::getMean() const {
519
520

	T mean = 0;
521
	const unsigned int nxyz = nx*ny*nz;
522

523
	if(keep_nodata==false) {
524
		for (unsigned int jj=0; jj<nxyz; jj++) {
525
			const T val = vecData[jj];
526
			mean += val;
527
		}
528
		if(nxyz>0) return mean/(T)(nxyz);
529
		else return (T)0;
530
	} else {
531
		unsigned int count = 0;
532
		for (unsigned int jj=0; jj<nxyz; jj++) {
533
			const T val = vecData[jj];
534
535
536
			if(val!=IOUtils::nodata) {
				mean += val;
				count++;
537
538
539
540
541
542
543
			}
		}
		if(count>0) return mean/(T)(count);
		else return (T)IOUtils::nodata;
	}
}

544
template<class T> size_t Array3D<T>::getCount() const
545
{
546
547
	const unsigned int nxyz = nx*ny*nz;

548
	if(keep_nodata==false) {
549
		return (size_t)nxyz;
550
	} else {
551
		size_t count = 0;
552
		for (unsigned int ii=0; ii<nxyz; ii++) {
553
			if(vecData[ii]!=IOUtils::nodata) count++;
554
		}
555
		return count;
556
557
558
	}
}

559
template<class T> void Array3D<T>::abs() {
560
561
	if(std::numeric_limits<T>::is_signed) {
		const unsigned int nxyz = nx*ny*nz;
562
		if(keep_nodata==false) {
563
564
			for (unsigned int jj=0; jj<nxyz; jj++) {
				T& val = vecData[jj];
565
566
				if(val<0) val=-val;
			}
567
		} else {
568
569
			for (unsigned int jj=0; jj<nxyz; jj++) {
				T& val = vecData[jj];
570
571
572
573
574
575
				if(val<0 && val!=IOUtils::nodata) val=-val;
			}
		}
	}
}

576
template<class T> const Array3D<T> Array3D<T>::getAbs() const {
577
	Array3D<T> result = *this; //make a copy
578
	result.abs(); //already implemented
579
580
581
582

	return result;
}

583
//arithmetic operators
584
585
586
587
588
589
590
591
592
593
594
595
596
597
template<class T> bool Array3D<T>::checkEpsilonEquality(const Array3D<double>& rhs, const double& epsilon) const {
	if(nx!=rhs.nx || ny!=rhs.ny || nz!=rhs.nz) return false;

	const unsigned int nxyz = nx*ny*nz;
	for (unsigned int jj=0; jj<nxyz; jj++)
		if(IOUtils::checkEpsilonEquality(vecData[jj], rhs.vecData[jj], epsilon)==false) return false;

	return true;
}

template<class T> bool Array3D<T>::checkEpsilonEquality(const Array3D<double>& rhs1, const Array3D<double>& rhs2, const double& epsilon) {
	return rhs1.checkEpsilonEquality(rhs2, epsilon);
}

598
599
template<class T> Array3D<T>& Array3D<T>::operator=(const Array3D<T>& source) {
	if(this != &source) {
600
		keep_nodata = source.keep_nodata;
601
602
603
604
605
606
607
608
609
		nx = source.nx;
		ny = source.ny;
		nz = source.nz;
		nxny = source.nxny;
		vecData = source.vecData;
	}
	return *this;
}

610
611
612
613
614
template<class T> Array3D<T>& Array3D<T>::operator=(const T& value) {
	std::fill(vecData.begin(), vecData.end(), value);
	return *this;
}

615
616
617
template<class T> Array3D<T>& Array3D<T>::operator+=(const Array3D<T>& rhs)
{
	//They have to have equal size
618
619
620
621
622
623
	if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
		std::stringstream ss;
		ss << "Trying to add two Array3D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << "," << nz << ") + (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
		throw IOException(ss.str(), AT);
	}
624
	//Add to every single member of the Array3D<T>
625
626
627
	const unsigned int nxyz = nx*ny*nz;
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxyz; jj++) {
628
			vecData[jj] += rhs(jj);
629
630
631
		}
	} else {
		for (unsigned int jj=0; jj<nxyz; jj++) {
632
633
			if(vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				vecData[jj] = IOUtils::nodata;
634
			else
635
				vecData[jj] += rhs(jj);
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
		}
	}

	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator+(const Array3D<T>& rhs)
{
	Array3D<T> result = *this; //make a copy
	result += rhs; //already implemented

	return result;
}

template<class T> Array3D<T>& Array3D<T>::operator+=(const T& rhs)
{
	//Add to every single member of the Array3D<T>
653
654
655
	const unsigned int nxyz = nx*ny*nz;
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxyz; jj++) {
656
			vecData[jj] += rhs;
657
658
659
		}
	} else {
		for (unsigned int jj=0; jj<nxyz; jj++) {
660
661
			if(vecData[jj]!=IOUtils::nodata)
				vecData[jj] += rhs;
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
		}
	}

	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator+(const T& rhs)
{
	Array3D<T> result = *this;
	result += rhs; //already implemented

	return result;
}

template<class T> Array3D<T>& Array3D<T>::operator-=(const Array3D<T>& rhs)
{
	//They have to have equal size
679
680
681
682
683
684
	if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
		std::stringstream ss;
		ss << "Trying to substract two Array3D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << "," << nz << ") - (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
		throw IOException(ss.str(), AT);
	}
685
	//Substract to every single member of the Array3D<T>
686
687
688
	const unsigned int nxyz = nx*ny*nz;
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxyz; jj++) {
689
			vecData[jj] -= rhs(jj);
690
691
692
		}
	} else {
		for (unsigned int jj=0; jj<nxyz; jj++) {
693
694
			if(vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				vecData[jj] = IOUtils::nodata;
695
			else
696
				vecData[jj] -= rhs(jj);
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
		}
	}

	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator-(const Array3D<T>& rhs)
{
	Array3D<T> result = *this; //make a copy
	result -= rhs; //already implemented

	return result;
}

template<class T> Array3D<T>& Array3D<T>::operator-=(const T& rhs)
{
713
	*this += -rhs;
714
715
716
717
718
719
	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator-(const T& rhs)
{
	Array3D<T> result = *this;
720
	result += -rhs; //already implemented
721
722
723
724
725
726
727

	return result;
}

template<class T> Array3D<T>& Array3D<T>::operator*=(const Array3D<T>& rhs)
{
	//They have to have equal size
728
729
730
731
732
733
	if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
		std::stringstream ss;
		ss << "Trying to multiply two Array3D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << "," << nz << ") * (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
		throw IOException(ss.str(), AT);
	}
734
	//Multiply every single member of the Array3D<T>
735
736
737
	const unsigned int nxyz = nx*ny*nz;
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxyz; jj++) {
738
			vecData[jj] *= rhs(jj);
739
740
741
		}
	} else {
		for (unsigned int jj=0; jj<nxyz; jj++) {
742
743
			if(vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				vecData[jj] = IOUtils::nodata;
744
			else
745
				vecData[jj] *= rhs(jj);
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
		}
	}

	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator*(const Array3D<T>& rhs)
{
	Array3D<T> result = *this; //make a copy
	result *= rhs; //already implemented

	return result;
}

template<class T> Array3D<T>& Array3D<T>::operator*=(const T& rhs)
{
	//Multiply every single member of the Array3D<T>
763
764
765
	const unsigned int nxyz = nx*ny*nz;
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxyz; jj++) {
766
			vecData[jj] *= rhs;
767
768
769
		}
	} else {
		for (unsigned int jj=0; jj<nxyz; jj++) {
770
771
			if(vecData[jj]!=IOUtils::nodata)
				vecData[jj] *= rhs;
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
		}
	}

	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator*(const T& rhs)
{
	Array3D<T> result = *this;
	result *= rhs; //already implemented

	return result;
}

template<class T> Array3D<T>& Array3D<T>::operator/=(const Array3D<T>& rhs)
{
	//They have to have equal size
789
790
791
792
793
794
	if ((rhs.nx != nx) || (rhs.ny != ny) || (rhs.nz != nz)) {
		std::stringstream ss;
		ss << "Trying to divide two Array3D objects with different dimensions: ";
		ss << "(" << nx << "," << ny << "," << nz << ") / (" << rhs.nx << "," << rhs.ny << "," << rhs.nz << ")";
		throw IOException(ss.str(), AT);
	}
795
	//Divide every single member of the Array3D<T>
796
797
798
	const unsigned int nxyz = nx*ny*nz;
	if(keep_nodata==false) {
		for (unsigned int jj=0; jj<nxyz; jj++) {
799
			vecData[jj] /= rhs(jj);
800
801
802
		}
	} else {
		for (unsigned int jj=0; jj<nxyz; jj++) {
803
804
			if(vecData[jj]==IOUtils::nodata || rhs(jj)==IOUtils::nodata)
				vecData[jj] = IOUtils::nodata;
805
			else
806
				vecData[jj] /= rhs(jj);
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
		}
	}

	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator/(const Array3D<T>& rhs)
{
	Array3D<T> result = *this; //make a copy
	result /= rhs; //already implemented

	return result;
}

template<class T> Array3D<T>& Array3D<T>::operator/=(const T& rhs)
{
823
	*this *= (1./rhs);
824
825
826
827
828
829
	return *this;
}

template<class T> const Array3D<T> Array3D<T>::operator/(const T& rhs)
{
	Array3D<T> result = *this;
830
	result *= (1./rhs); //already implemented
831
832
833
834

	return result;
}

835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
template<class T> bool Array3D<T>::operator==(const Array3D<T>& in) const {
	unsigned int in_nx=in.getNx(), in_ny=in.getNy(), in_nz=in.getNz();

	if(nx!=in_nx || ny!=in_ny || nz!=in_nz)
		return false;

	const unsigned int nxyz = nx*ny*nz;
	for(unsigned int jj=0; jj<nxyz; jj++)
		if( !IOUtils::checkEpsilonEquality( vecData[jj] , in.vecData[jj], 1e-6) ) return false;

	return true;
}

template<class T> bool Array3D<T>::operator!=(const Array3D<T>& in) const {
	return !(*this==in);
}

852
853
854
} //end namespace mio

#endif