3d/datafield.hh
Go to the documentation of this file.
1 /* -*- mia-c++ -*-
2  *
3  * This file is part of MIA - a toolbox for medical image analysis
4  * Copyright (c) Leipzig, Madrid 1999-2016 Gert Wollny
5  *
6  * MIA is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with MIA; if not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef __MIA_3DDATAFIELD_HH
22 #define __MIA_3DDATAFIELD_HH 1
23 
24 #include <cstdio>
25 #include <vector>
26 #include <cmath>
27 #include <cassert>
28 
29 #include <mia/3d/vector.hh>
30 #include <mia/3d/defines3d.hh>
31 #include <mia/3d/iterator.hh>
32 #include <mia/2d/datafield.hh>
33 #include <mia/core/msgstream.hh>
34 #include <mia/core/parameter.hh>
35 #include <mia/core/typedescr.hh>
36 #include <miaconfig.h>
37 
39 
40 
41 #define DECLARE_EXTERN_ITERATORS(TYPE) \
42  extern template class EXPORT_3D range3d_iterator<std::vector<TYPE>::iterator>; \
43  extern template class EXPORT_3D range3d_iterator<std::vector<TYPE>::const_iterator>; \
44  extern template class EXPORT_3D range3d_iterator_with_boundary_flag<std::vector<TYPE>::iterator>; \
45  extern template class EXPORT_3D range3d_iterator_with_boundary_flag<std::vector<TYPE>::const_iterator>; \
46  extern template class EXPORT_3D range2d_iterator<std::vector<TYPE>::iterator>; \
47  extern template class EXPORT_3D range2d_iterator<std::vector<TYPE>::const_iterator>;
48 
49 
50 #ifdef __GNUC__
51 #pragma GCC diagnostic push
52 #ifndef __clang__
53 #pragma GCC diagnostic ignored "-Wattributes"
54 #endif
55 #endif
56 
59 DECLARE_EXTERN_ITERATORS(unsigned int);
62 DECLARE_EXTERN_ITERATORS(unsigned short);
63 DECLARE_EXTERN_ITERATORS(unsigned char );
64 DECLARE_EXTERN_ITERATORS(signed char);
66 
67 #ifdef LONG_64BIT
68 DECLARE_EXTERN_ITERATORS(signed long);
69 DECLARE_EXTERN_ITERATORS(unsigned long);
70 #endif
71 
74 
75 #ifdef __GNUC__
76 #pragma GCC diagnostic pop
77 #endif
78 
83 template <class T>
85 
86  typedef ::std::vector<typename __holder_type_dispatch<T>::type> data_array;
87 
88  typedef std::shared_ptr<data_array> ref_data_type;
89 
90 public:
91 
92 
95  void make_single_ref();
96 
101  bool holds_unique_data()const {
102  return m_data.unique();
103  }
104 
105 
107 
109  typedef typename data_array::iterator iterator;
110  typedef typename data_array::const_iterator const_iterator;
111  typedef typename data_array::const_reference const_reference;
112  typedef typename data_array::reference reference;
113  typedef typename data_array::const_pointer const_pointer;
114  typedef typename data_array::pointer pointer;
115  typedef typename data_array::value_type value_type;
116  typedef typename data_array::size_type size_type;
117  typedef typename data_array::difference_type difference_type;
118  typedef typename atomic_data<value_type>::type atomic_type;
119  typedef range3d_iterator<iterator> range_iterator;
120  typedef range3d_iterator<const_iterator> const_range_iterator;
121 
122  typedef range3d_iterator_with_boundary_flag<iterator> range_iterator_with_boundary_flag;
123  typedef range3d_iterator_with_boundary_flag<const_iterator> const_range_iterator_with_boundary_flag;
124 
125  typedef C3DBounds dimsize_type;
127 
137  class EXPORT_3D Range {
138  friend class T3DDatafield<T>;
139  friend class ConstRange;
140  public:
141 
143 
144  iterator begin();
145 
146  iterator end();
147 
148  private:
149  Range(const C3DBounds& start, const C3DBounds& end, T3DDatafield<T>& field);
150 
151  iterator m_begin;
152  iterator m_end;
153  };
154 
156  public:
157  friend class T3DDatafield<T>;
158 
160 
161  iterator begin() const;
162 
163  iterator end() const;
164 
165  private:
166  ConstRange(const C3DBounds& start, const C3DBounds& end, const T3DDatafield<T>& field);
167 
168  ConstRange(const Range& range);
169 
170  iterator m_begin;
171  iterator m_end;
172  };
173 
174 
175  T3DDatafield();
176 
178  T3DDatafield(const C3DBounds& _Size);
179 
184  T3DDatafield(const C3DBounds& size, const T *data);
185 
186 
191  T3DDatafield(const C3DBounds& size, const data_array& data);
192 
193 
195  T3DDatafield(const T3DDatafield<T>& org);
196 
198  virtual ~T3DDatafield();
199 
204  template <typename Out>
206 
208  template <typename Out>
209  T3DVector<Out> get_gradient(size_t x, size_t y, size_t z) const;
210 
212  template <typename Out>
213  T3DVector<Out> get_gradient(int index) const;
214 
216  value_type get_interpol_val_at(const T3DVector<float >& p) const;
217 
218  /* some rough interpolation using barycentric coordinates, needs less addition and
219  multiplications then tri-linear interp. but is usally of low quality
220  \remark this function may vanish
221  value_type get_barycent_interpol_val_at(const T3DVector<float >& p) const;
222  */
223 
225  value_type get_trilin_interpol_val_at(const T3DVector<float >& p) const;
226 
229  value_type get_block_avrg(const C3DBounds& Start, const C3DBounds& BlockSize) const;
230 
236 
238  const C3DBounds& get_size() const
239  {
240  return m_size;
241  }
242 
244  void clear();
245 
247  size_type size()const
248  {
249  return m_data->size();
250  }
251 
253  void swap(T3DDatafield& other);
254 
256  value_type get_avg();
257 
260  value_type strip_avg();
261 
262 
264  value_type operator()(const T3DVector<float >& pos)const;
265 
267  const_reference operator()(size_t x, size_t y, size_t z) const
268  {
269  // Look if we are inside, and give reference, else give the zero
270  if (x < m_size.x && y < m_size.y && z < m_size.z) {
271  const data_array& data = *m_data;
272  return data[x+ m_size.x * (y + m_size.y * z)];
273  }
274  return Zero;
275  }
276 
277 
279  const_reference operator()(const C3DBounds& l)const
280  {
281  return (*this)(l.x,l.y,l.z);
282  }
283 
285  reference operator()(size_t x, size_t y, size_t z)
286  {
287  // Look if we are inside, and give reference, else throw exception
288  // since write access is wanted
289  assert(x < m_size.x && y < m_size.y && z < m_size.z);
290  return (*m_data)[x + m_size.x *(y + m_size.y * z)];
291  }
292 
293 
294 
296  reference operator()(const C3DBounds& l)
297  {
298  return (*this)(l.x,l.y,l.z);
299  }
300 
302  void get_data_line_x(int y, int z, std::vector<T>& buffer)const;
303 
305  void get_data_line_y(int x, int z, std::vector<T>& buffer)const;
306 
308  void get_data_line_z(int x, int y, std::vector<T>& buffer)const;
309 
311  void put_data_line_x(int y, int z, const std::vector<T> &buffer);
312 
314  void put_data_line_y(int x, int z, const std::vector<T> &buffer);
315 
317  void put_data_line_z(int x, int y, const std::vector<T> &buffer);
318 
320  template <class TMask>
321  void mask(const TMask& m);
322 
336  void read_xslice_flat(size_t x, std::vector<atomic_type>& buffer) const;
337 
350  void read_yslice_flat(size_t y, std::vector<atomic_type>& buffer) const;
351 
364  void read_zslice_flat(size_t z, std::vector<atomic_type>& buffer) const;
365 
370  void write_zslice_flat(size_t z, const std::vector<atomic_type>& buffer);
371 
372 
377  void write_yslice_flat(size_t y, const std::vector<atomic_type>& buffer);
378 
383  void write_xslice_flat(size_t x, const std::vector<atomic_type>& buffer);
384 
390  T2DDatafield<T> get_data_plane_xy(size_t z)const;
391 
397  T2DDatafield<T> get_data_plane_yz(size_t x)const;
398 
404  T2DDatafield<T> get_data_plane_xz(size_t y)const;
405 
411  void put_data_plane_xy(size_t z, const T2DDatafield<T>& p);
412 
418  void put_data_plane_yz(size_t x, const T2DDatafield<T>& p);
419 
425  void put_data_plane_xz(size_t y, const T2DDatafield<T>& p);
426 
428  const_iterator begin()const
429  {
430  return m_data->begin();
431  }
432 
436  const_iterator begin_at(size_t x, size_t y, size_t z)const
437  {
438  return m_data->begin() + (z * m_size.y + y) * m_size.x + x;
439  }
440 
441 
445  const_iterator end()const
446  {
447  return m_data->end();
448  }
449 
453  iterator begin()
454  {
455  make_single_ref();
456  return m_data->begin();
457  }
458 
459  Range get_range(const C3DBounds& start, const C3DBounds& end);
460 
461  ConstRange get_range(const C3DBounds& start, const C3DBounds& end) const;
462 
465  range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end);
466 
468  range_iterator end_range(const C3DBounds& begin, const C3DBounds& end);
469 
470 
473  const_range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end)const;
474 
476  const_range_iterator end_range(const C3DBounds& begin, const C3DBounds& end)const;
477 
478 
480  range_iterator_with_boundary_flag begin_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end);
481 
483  range_iterator_with_boundary_flag end_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end);
484 
485 
487  const_range_iterator_with_boundary_flag begin_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end)const;
488 
490  const_range_iterator_with_boundary_flag end_range_with_boundary_flags(const C3DBounds& begin, const C3DBounds& end)const;
491 
492 
501  iterator begin_at(size_t x, size_t y, size_t z)
502  {
503  make_single_ref();
504  return m_data->begin() + (z * m_size.y + y) * m_size.x + x;
505  }
506 
510  iterator end()
511  {
512  make_single_ref();
513  return m_data->end();
514  }
515 
517  const_reference operator[](int i)const
518  {
519  return (*m_data)[i];
520  }
521 
525  reference operator[](int i)
526  {
527  assert(m_data.unique());
528  return (*m_data)[i];
529  }
530 
531 
533  size_t get_plane_size_xy()const
534  {
535  return m_xy;
536  };
537 
538 private:
540  C3DBounds m_size;
541 
543  size_t m_xy;
544 
546  ref_data_type m_data;
547 
549  static const value_type Zero;
550 
551  static const size_t m_elements;
552 
553 };
554 
557 
560 
563 
564 
567 
570 
573 
576 
579 
582 
584 
586 DECLARE_TYPE_DESCR(C3DBounds);
587 DECLARE_TYPE_DESCR(C3DFVector);
588 
589 extern template class EXPORT_3D TAttribute<C3DFVector>;
591 
592 // some implementations
593 
594 template <class T>
595 template <typename Out>
596 T3DVector<Out> T3DDatafield<T>::get_gradient(size_t x, size_t y, size_t z) const
597 {
598  const std::vector<T>& data = *m_data;
599  const int sizex = m_size.x;
600  // Look if we are inside the used space
601  if (x - 1 < m_size.x - 2 && y - 1 < m_size.y - 2 && z - 1 < m_size.z - 2) {
602 
603  // Lookup all neccessary Values
604  const T *H = &data[x + m_size.x * (y + m_size.y * z)];
605 
606  return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
607  Out((H[sizex] - H[-sizex]) * 0.5),
608  Out((H[m_xy] - H[-m_xy]) * 0.5));
609  }
610 
611  return T3DVector<Out>();
612 }
613 
614 
615 template <class T>
616 template <typename Out>
618 {
619  const int sizex = m_size.x;
620  // Lookup all neccessary Values
621  const T *H = &(*m_data)[hardcode];
622 
623  return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
624  Out((H[sizex] - H[-sizex]) * 0.5),
625  Out((H[m_xy] - H[-m_xy]) * 0.5));
626 }
627 
628 
632 template <>
633 template <typename Out>
635 {
636 
637  // Lookup all neccessary Values
638  return T3DVector<Out> (Out(((*m_data)[hardcode + 1] - (*m_data)[hardcode -1]) * 0.5),
639  Out(((*m_data)[hardcode + m_size.x] - (*m_data)[hardcode -m_size.x]) * 0.5),
640  Out(((*m_data)[hardcode + m_xy] - (*m_data)[hardcode -m_xy]) * 0.5));
641 }
642 
643 template <class T>
644 template <typename Out>
646 {
647  // This will become really funny
648  const int sizex = m_size.x;
649  const std::vector<T>& data = *m_data;
650  // Calculate the int coordinates near the POI
651  // and the distances
652  size_t x = size_t (p.x);
653  float dx = p.x - x;
654  float xm = 1 - dx;
655  size_t y = size_t (p.y);
656  float dy = p.y - y;
657  float ym = 1 - dy;
658  size_t z = size_t (p.z);
659  float dz = p.z - z;
660  float zm = 1 - dz;
661 
662  // Look if we are inside the used space
663  if (x-1 < m_size.x-3 && y -1 < m_size.y-3 && z - 1 < m_size.z-3 ) {
664  // Lookup all neccessary Values
665  const T *H000 = &data[x + sizex * y + m_xy * z];
666 
667  const T* H_100 = &H000[-m_xy];
668  const T* H_101 = &H_100[1];
669  const T* H_110 = &H_100[sizex];
670  const T* H_111 = &H_110[1];
671 
672  const T* H0_10 = &H000[-sizex];
673  const T* H0_11 = &H0_10[1];
674 
675  const T* H00_1 = &H000[-1];
676  const T* H001 = &H000[ 1];
677  const T* H002 = &H000[ 2];
678 
679 
680  const T* H010 = &H000[sizex];
681  const T* H011 = &H010[ 1];
682  const T* H012 = &H010[ 2];
683  const T* H01_1 = &H010[-1];
684 
685  const T* H020 = &H010[sizex];
686  const T* H021 = &H020[ 1];
687 
688  const T* H100 = &H000[m_xy];
689 
690  const T* H1_10 = &H100[sizex];
691  const T* H1_11 = &H1_10[1];
692 
693  const T* H10_1 = &H100[-1];
694  const T* H101 = &H100[ 1];
695  const T* H102 = &H100[ 2];
696 
697  const T* H110 = &H100[sizex];
698  const T* H111 = &H110[ 1];
699  const T* H112 = &H110[ 2];
700  const T* H11_1 = &H110[-1];
701 
702  const T* H120 = &H110[sizex];
703  const T* H121 = &H120[ 1];
704 
705  const T* H200 = &H100[m_xy];
706  const T* H201 = &H200[1];
707  const T* H210 = &H200[sizex];
708  const T* H211 = &H210[1];
709 
710  // use trilinear interpolation to calc the gradient
711  return T3DVector<Out> (
712  Out((zm * (ym * (dx * (*H002 - *H000) + xm * (*H001 - *H00_1))+
713  dy * (dx * (*H012 - *H010) + xm * (*H011 - *H01_1)))+
714  dz * (ym * (dx * (*H102 - *H100) + xm * (*H101 - *H10_1))+
715  dy * (dx * (*H112 - *H110) + xm * (*H111 - *H11_1)))) * 0.5),
716 
717  Out((zm * (ym * (xm * (*H010 - *H0_10) + dx * (*H011 - *H0_11))+
718  dy * (xm * (*H020 - *H000) + dx * (*H021 - *H001)))+
719  dz * (ym * (xm * (*H110 - *H1_10) + dx * (*H111 - *H1_11))+
720  dy * (xm * (*H120 - *H100) + dx * (*H121 - *H101)))) * 0.5),
721  Out((zm * (ym * (xm * (*H100 - *H_100) + dx * (*H101 - *H_101))+
722  dy * (xm * (*H110 - *H_110) + dx * (*H111 - *H_111)))+
723  dz * (ym * (xm * (*H200 - *H000) + dx * (*H201 - *H001))+
724  dy * (xm * (*H210 - *H010) + dx * (*H211 - *H011)))) * 0.5));
725  }
726  return T3DVector<Out>();
727 
728 }
729 
730 #ifdef __GNUC__
731 #pragma GCC diagnostic push
732 #ifndef __clang__
733 #pragma GCC diagnostic ignored "-Wattributes"
734 #endif
735 #endif
736 
737 #define DECLARE_EXTERN(TYPE) \
738  extern template class EXPORT_3D T3DDatafield<TYPE>;
739 
740 
741 DECLARE_EXTERN(double);
742 DECLARE_EXTERN(float);
743 DECLARE_EXTERN(unsigned int);
744 DECLARE_EXTERN(int);
745 DECLARE_EXTERN(short);
746 DECLARE_EXTERN(unsigned short);
747 DECLARE_EXTERN(unsigned char );
748 DECLARE_EXTERN(signed char);
749 
750 #ifdef LONG_64BIT
751 DECLARE_EXTERN(signed long);
752 DECLARE_EXTERN(unsigned long);
753 #endif
754 
757 
758 extern template class EXPORT_3D CTParameter<C3DBounds>;
759 extern template class EXPORT_3D CTParameter<C3DFVector>;
760 extern template class EXPORT_3D TTranslator<C3DFVector>;
761 extern template class EXPORT_3D TAttribute<C3DFVector>;
762 
763 
764 #undef DECLARE_EXTERN
765 
766 #ifdef __GNUC__
767 #pragma GCC diagnostic pop
768 #endif
769 
771 
772 #endif
T2DDatafield< T > get_data_plane_yz(size_t x) const
T3DDatafield< T >::range_iterator iterator
void put_data_plane_xy(size_t z, const T2DDatafield< T > &p)
range_iterator begin_range(const C3DBounds &begin, const C3DBounds &end)
T3DVector< Out > get_gradient(const T3DVector< float > &p) const
T2DDatafield< T > get_data_plane_xz(size_t y) const
a 3D iterator that knows its position in the 3D grid ans supports iterating over sub-ranges ...
Definition: 3d/iterator.hh:189
void swap(T3DDatafield &other)
swap the data ofthis 3DDatafield with another one
Generic string vs. attribute translator singleton.
Definition: attributes.hh:509
void write_yslice_flat(size_t y, const std::vector< atomic_type > &buffer)
void make_single_ref()
void put_data_line_y(int x, int z, const std::vector< T > &buffer)
range_iterator_with_boundary_flag begin_range_with_boundary_flags(const C3DBounds &begin, const C3DBounds &end)
T z
vector element
Definition: 3d/vector.hh:55
a 3D iterator that knows its position in the 3D grid, has a flag indicating whether it is on a bounda...
Definition: 3d/iterator.hh:43
value_type get_trilin_interpol_val_at(const T3DVector< float > &p) const
A templated class of a 3D data field.
Definition: 3d/datafield.hh:84
iterator begin_at(size_t x, size_t y, size_t z)
const_iterator begin_at(size_t x, size_t y, size_t z) const
void read_zslice_flat(size_t z, std::vector< atomic_type > &buffer) const
T3DDatafield< unsigned char > C3DUBDatafield
a data field of float values
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
size_type size() const
T3DDatafield< bool > C3DBitDatafield
a data field of float values
#define EXPORT_3D
Definition: defines3d.hh:44
reference operator[](int i)
CTParameter< C3DBounds > C3DBoundsParameter
3D size parameter type
void read_yslice_flat(size_t y, std::vector< atomic_type > &buffer) const
const_reference operator()(const C3DBounds &l) const
value_type strip_avg()
TTranslator< C3DFVector > C3DFVectorTranslator
Range get_range(const C3DBounds &start, const C3DBounds &end)
T3DDatafield< unsigned int > C3DUIDatafield
a data field of 32 bit unsigned int values
void get_data_line_z(int x, int y, std::vector< T > &buffer) const
Generic type of a complex paramter.
Definition: parameter.hh:163
const_reference operator[](int i) const
T3DDatafield< long > C3DLDatafield
a data field of 32 bit signed int values
EXPORT_2D C2DFVectorfield get_gradient(const C2DImage &image)
T y
vector element
Definition: 3d/vector.hh:53
void write_xslice_flat(size_t x, const std::vector< atomic_type > &buffer)
T3DDatafield< T >::const_range_iterator iterator
CTParameter< C3DFVector > C3DFVectorParameter
3D vector parameter type
const C3DBounds & get_size() const
#define DECLARE_EXTERN(TYPE)
size_t get_plane_size_xy() const
const_reference operator()(size_t x, size_t y, size_t z) const
A class to hold data on a regular 2D grid.
Definition: 2d/datafield.hh:53
T2DDatafield< T > get_data_plane_xy(size_t z) const
void put_data_line_z(int x, int y, const std::vector< T > &buffer)
void write_zslice_flat(size_t z, const std::vector< atomic_type > &buffer)
const_iterator end() const
#define DECLARE_EXTERN_ITERATORS(TYPE)
Definition: 3d/datafield.hh:41
T3DDatafield & operator=(const T3DDatafield &org)
void read_xslice_flat(size_t x, std::vector< atomic_type > &buffer) const
void put_data_plane_xz(size_t y, const T2DDatafield< T > &p)
range_iterator_with_boundary_flag end_range_with_boundary_flags(const C3DBounds &begin, const C3DBounds &end)
void mask(const TMask &m)
a shortcut data type
iterator begin()
T3DDatafield< int > C3DIDatafield
a data field of 32 bit signed int values
virtual ~T3DDatafield()
make sure the destructor is virtual
value_type get_avg()
range_iterator end_range(const C3DBounds &begin, const C3DBounds &end)
value_type get_interpol_val_at(const T3DVector< float > &p) const
void get_data_line_x(int y, int z, std::vector< T > &buffer) const
reference operator()(size_t x, size_t y, size_t z)
value_type operator()(const T3DVector< float > &pos) const
iterator end()
void put_data_plane_yz(size_t x, const T2DDatafield< T > &p)
void get_data_line_y(int x, int z, std::vector< T > &buffer) const
value_type get_block_avrg(const C3DBounds &Start, const C3DBounds &BlockSize) const
reference operator()(const C3DBounds &l)
const_iterator begin() const
T3DDatafield< unsigned long > C3DULDatafield
a data field of 32 bit unsigned int values
void put_data_line_x(int y, int z, const std::vector< T > &buffer)
T3DDatafield< float > C3DFDatafield
a data field of float values
T x
vector element
Definition: 3d/vector.hh:51
bool holds_unique_data() const
Class of an attribute that holds data of type T.
Definition: attributes.hh:116
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36