This file is indexed.

/usr/include/mia-2.4/mia/3d/datafield.hh is in libmia-2.4-dev 2.4.6-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  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
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
/* -*- mia-c++  -*-
 *
 * This file is part of MIA - a toolbox for medical image analysis 
 * Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
 *
 * MIA is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program 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 General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with MIA; if not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef __MIA_3DDATAFIELD_HH
#define __MIA_3DDATAFIELD_HH 1

#include <cstdio>
#include <vector>
#include <cmath>
#include <cassert>

#include <mia/3d/vector.hh>
#include <mia/3d/defines3d.hh>
#include <mia/3d/iterator.hh>
#include <mia/2d/datafield.hh>
#include <mia/core/msgstream.hh>
#include <mia/core/parameter.hh>
#include <mia/core/typedescr.hh>
#include <miaconfig.h>

NS_MIA_BEGIN


#define DECLARE_EXTERN_ITERATORS(TYPE)						\
	extern template class  EXPORT_3D range3d_iterator<std::vector<TYPE>::iterator>; \
	extern template class  EXPORT_3D range3d_iterator<std::vector<TYPE>::const_iterator>; \
	extern template class  EXPORT_3D range3d_iterator_with_boundary_flag<std::vector<TYPE>::iterator>; \
	extern template class  EXPORT_3D range3d_iterator_with_boundary_flag<std::vector<TYPE>::const_iterator>; \
	extern template class  EXPORT_3D range2d_iterator<std::vector<TYPE>::iterator>; \
	extern template class  EXPORT_3D range2d_iterator<std::vector<TYPE>::const_iterator>;


#ifdef __GNUC__
#pragma GCC diagnostic push
#ifndef __clang__
#pragma GCC diagnostic ignored "-Wattributes"
#endif
#endif

DECLARE_EXTERN_ITERATORS(double);
DECLARE_EXTERN_ITERATORS(float);
DECLARE_EXTERN_ITERATORS(uint32_t);
DECLARE_EXTERN_ITERATORS(int32_t);
DECLARE_EXTERN_ITERATORS(int16_t);
DECLARE_EXTERN_ITERATORS(uint16_t);
DECLARE_EXTERN_ITERATORS(int8_t);
DECLARE_EXTERN_ITERATORS(uint8_t);
DECLARE_EXTERN_ITERATORS(bool);
DECLARE_EXTERN_ITERATORS(int64_t);
DECLARE_EXTERN_ITERATORS(uint64_t);


DECLARE_EXTERN_ITERATORS(C3DFVector)
DECLARE_EXTERN_ITERATORS(C3DDVector)

#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif 

/**
   @ingroup basic 
   \brief A templated class of a 3D data field.
*/
template <class T>
class  EXPORT_3D T3DDatafield {

	typedef  ::std::vector<typename __holder_type_dispatch<T>::type> data_array;

public:
	

        /** makes a single reference of the data, after calling this, it is save to write to the data field
         */
        void make_single_ref() __attribute__((deprecated));

	/**
	   Checks whether the data hold by the data field is unique. 
	   \returns true if it is 
	 */
	bool holds_unique_data()const __attribute__((deprecated)){ 
		return true; 
	}
			

	/// a shortcut data type

	/// \cond SELFEXPLAINING 
        typedef typename data_array::iterator iterator;
        typedef typename data_array::const_iterator const_iterator;
        typedef typename data_array::const_reference const_reference;
        typedef typename data_array::reference reference;
        typedef typename data_array::const_pointer const_pointer;
        typedef typename data_array::pointer pointer;
        typedef typename data_array::value_type value_type;
        typedef typename data_array::size_type size_type;
        typedef typename data_array::difference_type difference_type;
	typedef typename atomic_data<value_type>::type atomic_type; 
	typedef range3d_iterator<iterator> range_iterator; 
	typedef range3d_iterator<const_iterator> const_range_iterator; 

	typedef range3d_iterator_with_boundary_flag<iterator> range_iterator_with_boundary_flag; 
	typedef range3d_iterator_with_boundary_flag<const_iterator> const_range_iterator_with_boundary_flag; 

	typedef C3DBounds dimsize_type;
	/// \endcond 

	/**
	   \brief This class provides access to a sub-range of the input data field

	   This class provides iterator access to a axis-aligned 3D sub-range of the 
	   corresponding data field. 
	*/

	
	
	class EXPORT_3D Range {
		friend class T3DDatafield<T>;
		friend class ConstRange;
	public:
		
		typedef T3DDatafield<T>::range_iterator iterator;
		
		iterator begin();
		
		iterator end();
		
	private:
		Range(const C3DBounds& start, const C3DBounds& end, T3DDatafield<T>& field);

		iterator m_begin;
		iterator m_end;
	}; 
	
	class EXPORT_3D ConstRange {
	public:
		friend class T3DDatafield<T>;

		typedef T3DDatafield<T>::const_range_iterator iterator;
		
		iterator begin() const;
		
		iterator end() const;

	private:
		ConstRange(const C3DBounds& start, const C3DBounds& end, const T3DDatafield<T>& field);

		ConstRange(const Range& range); 
		
		iterator m_begin;
		iterator m_end;
	}; 

	
	T3DDatafield();

        /** Constructor to create empty Datafield if given size */
        explicit T3DDatafield(const C3DBounds& _Size);

        /** Constructor to create Datafield if given size and with initialization data
            \param size the size of the 3D-field
            \param data to use for initialization
         */
        T3DDatafield(const C3DBounds& size, const T *data);


        /** Constructor to create Datafield if given size and with initialization data
            \param size the size of the 3D-field
            \param data to use for initialization
         */
        T3DDatafield(const C3DBounds& size, const data_array& data);


        /** copy - Constructor */
        T3DDatafield(const T3DDatafield<T>& org);

	/** move constructor */
	T3DDatafield(T3DDatafield<T>&& org);

        /// make sure the destructor is virtual
        virtual ~T3DDatafield();

        /**
            Gradient calculation using tri-linear interpolation
            \param p position where to evaluate the gradient
        */
	template <typename Out>
        T3DVector<Out> get_gradient(const T3DVector<float >& p) const;

        /** calculate gradient of data field at a grid point */
        template <typename Out>
	T3DVector<Out> get_gradient(size_t  x, size_t  y, size_t  z) const;

        /** calculate the gradient at a grid point given by a linear location */
        template <typename Out>
	T3DVector<Out> get_gradient(int index) const;

        /** Interpolate the value of Field at p default uses tri-linear interpolation */
        value_type get_interpol_val_at(const T3DVector<float >& p) const __attribute__((deprecated));

        /** Get the average over a given Block
         Attn: Type T must be able to hold the Sum of all Elements in Block */
        value_type get_block_avrg(const C3DBounds& Start, const C3DBounds& BlockSize) const;

        /** Assignment operator  -
            \remark it just copys a pointer to the data and increases its reference count,
            before writing it is necesary to call \a make_single_ref
        */
        T3DDatafield& operator = (const T3DDatafield& org);

	/// Moave asignment 
	T3DDatafield& operator = (T3DDatafield&& org);

        /** \returns the 3D-size of the data field */
        const C3DBounds&  get_size() const
        {
                return m_size;
        }

        /** Set alle elements of the field to T() == Zero*/
        void clear();

        /** \returns the number of elements in the datafield */
        size_type size()const
        {
                return m_data.size();
        }

	/// swap the data ofthis 3DDatafield with another one 
	void swap(T3DDatafield& other);

        /** \returns the average over the whole datafield*/
        value_type get_avg();

        /** Strip average from data
         \returns the stripped average */
        value_type strip_avg();

        /** read-only indx operator */
        const_reference operator()(size_t  x, size_t  y, size_t  z) const
	{
        	// Look if we are inside, and give reference, else give the zero
	        if (x < m_size.x && y < m_size.y && z < m_size.z) {
	                return m_data[x+ m_size.x * (y  + m_size.y * z)];
	        }
		return Zero;
	}


        /** alternate read-only indx operator */
        const_reference operator()(const C3DBounds& l)const
        {
                return (*this)(l.x,l.y,l.z);
        }

        /** Index operator witch gives write access */
        reference operator()(size_t  x, size_t  y, size_t  z)
	{
        	// Look if we are inside, and give reference, else throw exception
	        // since write access is wanted
	        assert(x < m_size.x && y < m_size.y && z < m_size.z);
		return m_data[x + m_size.x *(y + m_size.y * z)];
	}



        /** Alternate index operator witch gives write access */
        reference operator()(const C3DBounds& l)
        {
                return (*this)(l.x,l.y,l.z);
        }

        /** Get some Data along some line parallel to X axis */
        void get_data_line_x(int y, int z, std::vector<T>& buffer)const;

        /** Get some Data along some line parallel to Y axis */
        void get_data_line_y(int x, int z, std::vector<T>& buffer)const;

        /** Get some Data along some line parallel to Z axis */
        void get_data_line_z(int x, int y, std::vector<T>& buffer)const;

        /** Put some Data along some line parallel to X axis */
        void put_data_line_x(int y, int z, const std::vector<T> &buffer);

        /** Put some Data along some line parallel to Y axis */
        void put_data_line_y(int x, int z, const std::vector<T> &buffer);

        /** Put some Data along some line parallel to Z axis */
        void put_data_line_z(int x, int y, const std::vector<T> &buffer);

        /** Mask the data field with a given mask */
        template <class TMask>
        void mask(const TMask& m);

	/**
	   Read the a x-slice of the data field into a flat buffer - i.e. the 
           information about multi-dimensionality of the elements is lost. 
	   For this to work, T has to be a POD-like data type, i.e., it has no 
	   hidden elements like a virtual methods table, and, if T is a type 
	   of more then one element, all these elements have to be of the same 
	   type. Specifically, a specialization of the trait atomic_data for T 
	   must exists. 
	   \param x slice to be read 
	   \param[out] buffer Buffer where the data will be written to. It must 
	   large enough to hold size.y * size.z * number of elements 
	   
	*/
	void read_xslice_flat(size_t x, std::vector<atomic_type>& buffer) const;

	/**
	   Read the a y-slice of the data field into a flat buffer - i.e. the 
           information about multi-dimensionality of the elements is lost. 
	   For this to work, T has to be a POD-like data type, i.e., it has no 
	   hidden elements like a virtual methods table, and, if T is a type 
	   of more then one element, all these elements have to be of the same 
	   type. Specifically, a specialization of the trait atomic_data for T 
	   must exists. 
	   \param y slice to be read 
	   \param[out] buffer Buffer where the data will be written to. It must 
	   large enough to hold size.x * size.z * number of elements 
	*/
	void read_yslice_flat(size_t y, std::vector<atomic_type>& buffer) const;
	
	/**
	   Read the a z-slice of the data field into a flat buffer - i.e. the 
           information about multi-dimensionality of the elements is lost. 
	   For this to work, T has to be a POD-like data type, i.e., it has no 
	   hidden elements like a virtual methods table, and, if T is a type 
	   of more then one element, all these elements have to be of the same 
	   type. Specifically, a specialization of the trait atomic_data for T 
	   must exists. 
	   \param z slice to be read 
	   \param[out] buffer Buffer where the data will be written to. It must 
	   large enough to hold size.x * size.y * number of elements 
	*/
	void read_zslice_flat(size_t z, std::vector<atomic_type>& buffer) const;
	
	/**
	   Write a z-slice from a flat buffer to the 3D data field. For details see 
	   void read_zslice_flat(size_t z, std::vector<atomic_type>& buffer) const;
	 */
	void write_zslice_flat(size_t z, const std::vector<atomic_type>& buffer); 


	/**
	   Write a y-slice from a flat buffer to the 3D data field. For details see 
	   void read_yslice_flat(size_t y, std::vector<atomic_type>& buffer) const;
	 */
	void write_yslice_flat(size_t y, const std::vector<atomic_type>& buffer); 

	/**
	   Write a x-slice from a flat buffer to the 3D data field. For details see 
	   void read_yslice_flat(size_t x, std::vector<atomic_type>& buffer) const;
	*/
	void write_xslice_flat(size_t x, const std::vector<atomic_type>& buffer); 

	/**
	   Read a z-plane from the 3D data set.
	   \param z
	   \returns the copied data in a 2D data field 
	*/
	T2DDatafield<T> get_data_plane_xy(size_t  z)const;
	
	/**
	   Read a x-plane from the 3D data set.
	   \param x
	   \returns the copied data in a 2D data field 
	*/
	T2DDatafield<T> get_data_plane_yz(size_t  x)const;

	/**
	   Read a y-plane from the 3D data set.
	   \param y
	   \returns the copied data in a 2D data field 
	*/
	T2DDatafield<T> get_data_plane_xz(size_t  y)const;

	/**
	   write a z-plane to the 3D data set.
	   \param z
	   \param p plane data, must be of dimensions (size.x, size.y)
	*/
	void put_data_plane_xy(size_t  z, const T2DDatafield<T>& p);

	/**
	   write a x-plane to the 3D data set.
	   \param x
	   \param p plane data, must be of dimensions (size.y, size.z)
	*/
        void put_data_plane_yz(size_t  x, const T2DDatafield<T>& p);

	/**
	   write a y-plane to the 3D data set.
	   \param y
	   \param p plane data, must be of dimensions (size.x, size.z)
	*/
        void put_data_plane_xz(size_t  y, const T2DDatafield<T>& p);

        /** \returns an read only forward iterator over the whole data field */
        const_iterator begin()const
        {
                return m_data.begin();
        }
	
	/**
	   \returns an read only forward iterator over data field starting at (x,y,z)
	 */
	const_iterator begin_at(size_t x, size_t y, size_t z)const
        {
                return m_data.begin() + (z * m_size.y + y) * m_size.x + x;
        }


	/**
	   \returns the end iterator to the 3D data field 
	 */
        const_iterator end()const
        {
                return m_data.end();
        }

        /** \returns an read/write random access iterator over the whole data
            field pointing at the beginning of the data.
            The functions ensures, that the field uses a single referenced datafield */
        iterator begin()
        {
                return m_data.begin();
        }
	
	Range get_range(const C3DBounds& start, const C3DBounds& end);

	ConstRange get_range(const C3DBounds& start, const C3DBounds& end) const; 
	
/** \returns an read/write forward iterator over a subset of the data. 
            The functions ensures, that the field uses a single referenced datafield */
        range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end); 

        /** \returns the end of a read/write forward iterator over a subset of the data. */
        range_iterator end_range(const C3DBounds& begin, const C3DBounds& end); 


        /** \returns an read/write forward iterator over a subset of the data. 
            The functions ensures, that the field uses a single referenced datafield */
        const_range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end)const; 

        /** \returns the end of a read/write forward iterator over a subset of the data. */
        const_range_iterator end_range(const C3DBounds& begin, const C3DBounds& end)const; 


        /** \returns an read/write forward iterator over a subset of the data with indicator for the boundaries. */
	range_iterator_with_boundary_flag begin_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end); 

        /** \returns the end of a read/write forward iterator over a subset of the data with indicator for the boundaries. */
	range_iterator_with_boundary_flag end_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end); 


        /** \returns an read/write forward iterator over a subset of the data with indicator for the boundaries.  */
	const_range_iterator_with_boundary_flag begin_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end)const; 

        /** \returns the end of a read/write forward iterator over a subset of the data with indicator for the boundaries. */
	const_range_iterator_with_boundary_flag end_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end)const; 


	/**
	   Obtain an iterator at position (x,y,z)
	   The functions ensures, that the field uses a single referenced datafield
	   \param x
	   \param y
	   \param z
	   \returns the iterator 
	 */
	iterator begin_at(size_t x, size_t y, size_t z)
        {
		return m_data.begin() + (z * m_size.y + y) * m_size.x + x;
        }

	/** \returns an read/write random access iterator over the whole data
            field pointing at the end of the data.
            The functions ensures, that the field uses a single referenced datafield */
        iterator end()
        {
                return m_data.end();
        }

        /** a linear read only access operator */
        const_reference operator[](int i)const
        {
                return m_data[i];
        }

        /** A linear read/write access operator. The refcount of Data must be 1,
            else the program will abort with a failed assertion (if assert is enabled)
        */
        reference operator[](int i)
        {
                return m_data[i];
        }


        /** \returns the element count of one z slice */
        size_t  get_plane_size_xy()const
        {
                return m_xy;
        };

private:
	/** Size of the field */
        C3DBounds  m_size;

        /** helper: product of Size.x * Size.y */
        size_t  m_xy;

        /** Pointer to the Field of Data hold by this class */
        data_array m_data;

        /** helper: represents the zero-value */
        static const value_type Zero;
	
	static const size_t m_elements; 

};

/// a data field of float values
typedef T3DDatafield<float>  C3DFDatafield;

/// a data field of 32 bit unsigned int values
typedef T3DDatafield<uint32_t> C3DUIDatafield;

/// a data field of 32 bit signed int values
typedef T3DDatafield<int32_t>  C3DSIDatafield;

/// a data field of 32 bit unsigned int values
typedef T3DDatafield<uint16_t> C3DUSDatafield;

/// a data field of 32 bit signed int values
typedef T3DDatafield<int16_t>  C3DSSDatafield;

/// a data field of 32 bit unsigned int values
typedef T3DDatafield<uint64_t> C3DULDatafield;

/// a data field of 32 bit signed int values
typedef T3DDatafield<int64_t>  C3DLDatafield;

/// a data field of 8 bit int values
typedef T3DDatafield<uint8_t>  C3DUBDatafield;

/// a data field of 8 bit int values
typedef T3DDatafield<int8_t>  C3DSBDatafield;


	/// a data field of float values
typedef T3DDatafield<bool>  C3DBitDatafield;

/// 3D size parameter type 
typedef  CTParameter<C3DBounds> C3DBoundsParameter;

/// 3D vector parameter type 
typedef  CTParameter<C3DFVector> C3DFVectorParameter;

typedef  TTranslator<C3DFVector> C3DFVectorTranslator; 

/// @cond NEVER 
DECLARE_TYPE_DESCR(C3DBounds); 
DECLARE_TYPE_DESCR(C3DFVector);

extern template class EXPORT_3D TAttribute<C3DFVector>; 
/// @endcond 

// some implementations

template <class T>
template <typename Out>
T3DVector<Out> T3DDatafield<T>::get_gradient(size_t  x, size_t  y, size_t  z) const
{
	const int sizex = m_size.x;
	// Look if we are inside the used space
	if (x - 1 < m_size.x - 2 &&  y - 1 < m_size.y - 2 &&  z - 1 < m_size.z - 2) {

                // Lookup all neccessary Values
		const T *H  = &m_data[x + m_size.x * (y + m_size.y * z)];

		return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
				       Out((H[sizex] - H[-sizex]) * 0.5),
				       Out((H[m_xy] - H[-m_xy]) * 0.5));
	}

	return T3DVector<Out>();
}


template <class T>
template <typename Out>
T3DVector<Out> T3DDatafield<T>::get_gradient(int hardcode) const
{
	const int sizex = m_size.x;
	// Lookup all neccessary Values
	const T *H  = &m_data[hardcode];

	return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
			       Out((H[sizex] - H[-sizex]) * 0.5),
			       Out((H[m_xy] - H[-m_xy]) * 0.5));
}


/**
   Specialization to handle the wired std::vector<bool> implementation 
 */
template <>
template <typename Out>
T3DVector<Out> T3DDatafield<bool>::get_gradient(int hardcode) const
{

	// Lookup all neccessary Values
	return T3DVector<Out> (Out((m_data[hardcode + 1] - m_data[hardcode -1]) * 0.5),
			       Out((m_data[hardcode + m_size.x] - m_data[hardcode -m_size.x]) * 0.5),
			       Out((m_data[hardcode + m_xy] - m_data[hardcode -m_xy]) * 0.5));
}

template <class T>
template <typename Out>
T3DVector<Out> T3DDatafield<T>::get_gradient(const T3DVector<float >& p) const
{
        // This will become really funny
	const int sizex = m_size.x;
        // Calculate the int coordinates near the POI
        // and the distances
        size_t  x = size_t (p.x);
        float  dx = p.x - x;
        float  xm = 1 - dx;
        size_t  y = size_t (p.y);
        float  dy = p.y - y;
        float  ym = 1 - dy;
        size_t  z = size_t (p.z);
        float  dz = p.z - z;
        float  zm = 1 - dz;

	// Look if we are inside the used space
        if (x-1 < m_size.x-3 &&  y -1 < m_size.y-3 && z - 1 < m_size.z-3 ) {
                // Lookup all neccessary Values
                const T *H000  = &m_data[x + sizex * y + m_xy * z];

                const T* H_100 = &H000[-m_xy];
                const T* H_101 = &H_100[1];
                const T* H_110 = &H_100[sizex];
                const T* H_111 = &H_110[1];

                const T* H0_10 = &H000[-sizex];
                const T* H0_11 = &H0_10[1];

                const T* H00_1 = &H000[-1];
                const T* H001  = &H000[ 1];
                const T* H002  = &H000[ 2];


                const T* H010  = &H000[sizex];
                const T* H011  = &H010[ 1];
                const T* H012  = &H010[ 2];
                const T* H01_1 = &H010[-1];

                const T* H020 = &H010[sizex];
                const T* H021 = &H020[ 1];

                const T* H100 = &H000[m_xy];

                const T* H1_10 = &H100[sizex];
                const T* H1_11 = &H1_10[1];

                const T* H10_1 = &H100[-1];
                const T* H101  = &H100[ 1];
                const T* H102  = &H100[ 2];

                const T* H110  = &H100[sizex];
                const T* H111  = &H110[ 1];
                const T* H112  = &H110[ 2];
                const T* H11_1 = &H110[-1];

                const T* H120 = &H110[sizex];
                const T* H121 = &H120[ 1];

                const T* H200 = &H100[m_xy];
                const T* H201 = &H200[1];
                const T* H210 = &H200[sizex];
                const T* H211 = &H210[1];

                // use trilinear interpolation to calc the gradient
                return T3DVector<Out> (
			Out((zm * (ym * (dx * (*H002 - *H000) + xm * (*H001 - *H00_1))+
				   dy * (dx * (*H012 - *H010) + xm * (*H011 - *H01_1)))+
			     dz * (ym * (dx * (*H102 - *H100) + xm * (*H101 - *H10_1))+
				   dy * (dx * (*H112 - *H110) + xm * (*H111 - *H11_1)))) * 0.5),

			Out((zm * (ym * (xm * (*H010 - *H0_10) + dx * (*H011 - *H0_11))+
				   dy * (xm * (*H020 - *H000)  + dx * (*H021 - *H001)))+
			     dz * (ym * (xm * (*H110 - *H1_10) + dx * (*H111 - *H1_11))+
				   dy * (xm * (*H120 - *H100)  + dx * (*H121 - *H101)))) * 0.5),
			Out((zm * (ym * (xm * (*H100 - *H_100) + dx * (*H101 - *H_101))+
				   dy * (xm * (*H110 - *H_110) + dx * (*H111 - *H_111)))+
			     dz * (ym * (xm * (*H200 - *H000)  + dx * (*H201 - *H001))+
				   dy * (xm * (*H210 - *H010)  + dx * (*H211 - *H011)))) * 0.5));
        }
        return T3DVector<Out>();

}

#ifdef __GNUC__
#pragma GCC diagnostic push
#ifndef __clang__
#pragma GCC diagnostic ignored "-Wattributes"
#endif
#endif

#define DECLARE_EXTERN(TYPE)						\
	extern template class  EXPORT_3D T3DDatafield<TYPE>; 


DECLARE_EXTERN(double);
DECLARE_EXTERN(float);
DECLARE_EXTERN(uint8_t);
DECLARE_EXTERN(uint16_t);
DECLARE_EXTERN(uint32_t);
DECLARE_EXTERN(uint64_t);
DECLARE_EXTERN(int8_t);
DECLARE_EXTERN(int16_t);
DECLARE_EXTERN(int32_t);
DECLARE_EXTERN(int64_t);

DECLARE_EXTERN(C3DFVector);
DECLARE_EXTERN(C3DDVector);

extern template class EXPORT_3D CTParameter<C3DBounds>;
extern template class EXPORT_3D CTParameter<C3DFVector>;
extern template class EXPORT_3D TTranslator<C3DFVector>; 
extern template class EXPORT_3D TAttribute<C3DFVector>; 


#undef DECLARE_EXTERN

#ifdef __GNUC__
#pragma GCC diagnostic pop
#endif 

NS_MIA_END

#endif