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-2014 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 
44 template <class T>
46 
47  typedef ::std::vector<T> data_array;
48 
49  typedef std::shared_ptr<data_array> ref_data_type;
50 
52  C3DBounds m_size;
53 
55  size_t m_xy;
56 
58  ref_data_type m_data;
59 
61  static const T Zero;
62 
63  static const unsigned int m_elements;
64 
65 public:
66 
69  void make_single_ref();
70 
75  bool holds_unique_data()const {
76  return m_data.unique();
77  }
78 
79 
81 
83  typedef typename std::vector<T>::iterator iterator;
84  typedef typename std::vector<T>::const_iterator const_iterator;
85  typedef typename std::vector<T>::const_reference const_reference;
86  typedef typename std::vector<T>::reference reference;
87  typedef typename std::vector<T>::const_pointer const_pointer;
88  typedef typename std::vector<T>::pointer pointer;
89  typedef typename std::vector<T>::value_type value_type;
90  typedef typename std::vector<T>::size_type size_type;
91  typedef typename std::vector<T>::difference_type difference_type;
92  typedef typename atomic_data<T>::type atomic_type;
93  typedef range3d_iterator<iterator> range_iterator;
94  typedef range3d_iterator<const_iterator> const_range_iterator;
95 
96  typedef range3d_iterator_with_boundary_flag<iterator> range_iterator_with_boundary_flag;
97  typedef range3d_iterator_with_boundary_flag<const_iterator> const_range_iterator_with_boundary_flag;
98 
99  typedef C3DBounds dimsize_type;
101 
102  T3DDatafield();
103 
105  T3DDatafield(const C3DBounds& _Size);
106 
111  T3DDatafield(const C3DBounds& size, const T *data);
112 
113 
118  T3DDatafield(const C3DBounds& size, const data_array& data);
119 
120 
122  T3DDatafield(const T3DDatafield<T>& org);
123 
125  virtual ~T3DDatafield();
126 
131  template <typename Out>
133 
135  template <typename Out>
136  T3DVector<Out> get_gradient(size_t x, size_t y, size_t z) const;
137 
139  template <typename Out>
140  T3DVector<Out> get_gradient(int index) const;
141 
143  value_type get_interpol_val_at(const T3DVector<float >& p) const;
144 
145  /* some rough interpolation using barycentric coordinates, needs less addition and
146  multiplications then tri-linear interp. but is usally of low quality
147  \remark this function may vanish
148  value_type get_barycent_interpol_val_at(const T3DVector<float >& p) const;
149  */
150 
152  value_type get_trilin_interpol_val_at(const T3DVector<float >& p) const;
153 
156  value_type get_block_avrg(const C3DBounds& Start, const C3DBounds& BlockSize) const;
157 
162  T3DDatafield& operator = (const T3DDatafield& org);
163 
165  const C3DBounds& get_size() const
166  {
167  return m_size;
168  }
169 
171  void clear();
172 
174  size_type size()const
175  {
176  return m_data->size();
177  }
178 
180  void swap(T3DDatafield& other);
181 
183  value_type get_avg();
184 
187  value_type strip_avg();
188 
189 
191  value_type operator()(const T3DVector<float >& pos)const;
192 
194  const_reference operator()(size_t x, size_t y, size_t z) const
195  {
196  // Look if we are inside, and give reference, else give the zero
197  if (x < m_size.x && y < m_size.y && z < m_size.z) {
198  const std::vector<T>& data = *m_data;
199  return data[x+ m_size.x * (y + m_size.y * z)];
200  }
201  return Zero;
202  }
203 
204 
206  const_reference operator()(const C3DBounds& l)const
207  {
208  return (*this)(l.x,l.y,l.z);
209  }
210 
212  reference operator()(size_t x, size_t y, size_t z)
213  {
214  // Look if we are inside, and give reference, else throw exception
215  // since write access is wanted
216  assert(x < m_size.x && y < m_size.y && z < m_size.z);
217  return (*m_data)[x + m_size.x *(y + m_size.y * z)];
218  }
219 
220 
221 
223  reference operator()(const C3DBounds& l)
224  {
225  return (*this)(l.x,l.y,l.z);
226  }
227 
229  void get_data_line_x(int y, int z, std::vector<T>& buffer)const;
230 
232  void get_data_line_y(int x, int z, std::vector<T>& buffer)const;
233 
235  void get_data_line_z(int x, int y, std::vector<T>& buffer)const;
236 
238  void put_data_line_x(int y, int z, const std::vector<T> &buffer);
239 
241  void put_data_line_y(int x, int z, const std::vector<T> &buffer);
242 
244  void put_data_line_z(int x, int y, const std::vector<T> &buffer);
245 
247  template <class TMask>
248  void mask(const TMask& m);
249 
263  void read_xslice_flat(size_t x, std::vector<atomic_type>& buffer) const;
264 
277  void read_yslice_flat(size_t y, std::vector<atomic_type>& buffer) const;
278 
291  void read_zslice_flat(size_t z, std::vector<atomic_type>& buffer) const;
292 
297  void write_zslice_flat(size_t z, const std::vector<atomic_type>& buffer);
298 
299 
304  void write_yslice_flat(size_t y, const std::vector<atomic_type>& buffer);
305 
310  void write_xslice_flat(size_t x, const std::vector<atomic_type>& buffer);
311 
317  T2DDatafield<T> get_data_plane_xy(size_t z)const;
318 
324  T2DDatafield<T> get_data_plane_yz(size_t x)const;
325 
331  T2DDatafield<T> get_data_plane_xz(size_t y)const;
332 
338  void put_data_plane_xy(size_t z, const T2DDatafield<T>& p);
339 
345  void put_data_plane_yz(size_t x, const T2DDatafield<T>& p);
346 
352  void put_data_plane_xz(size_t y, const T2DDatafield<T>& p);
353 
355  const_iterator begin()const
356  {
357  return m_data->begin();
358  }
359 
363  const_iterator begin_at(size_t x, size_t y, size_t z)const
364  {
365  return m_data->begin() + (z * m_size.y + y) * m_size.x + x;
366  }
367 
368 
372  const_iterator end()const
373  {
374  return m_data->end();
375  }
376 
380  iterator begin()
381  {
382  make_single_ref();
383  return m_data->begin();
384  }
385 
388  range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end);
389 
391  range_iterator end_range(const C3DBounds& begin, const C3DBounds& end);
392 
393 
396  const_range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end)const;
397 
399  const_range_iterator end_range(const C3DBounds& begin, const C3DBounds& end)const;
400 
401 
410  iterator begin_at(size_t x, size_t y, size_t z)
411  {
412  make_single_ref();
413  return m_data->begin() + (z * m_size.y + y) * m_size.x + x;
414  }
415 
419  iterator end()
420  {
421  make_single_ref();
422  return m_data->end();
423  }
424 
426  const_reference operator[](int i)const
427  {
428  return (*m_data)[i];
429  }
430 
434  reference operator[](int i)
435  {
436  assert(m_data.unique());
437  return (*m_data)[i];
438  }
439 
440 
442  size_t get_plane_size_xy()const
443  {
444  return m_xy;
445  };
446 
447 private:
448 };
449 
450 
453 
456 
459 
460 
463 
466 
469 
472 
475 
478 
480 DECLARE_TYPE_DESCR(C3DBounds);
481 DECLARE_TYPE_DESCR(C3DFVector);
483 
484 // some implementations
485 
486 template <class T>
487 template <typename Out>
488 T3DVector<Out> T3DDatafield<T>::get_gradient(size_t x, size_t y, size_t z) const
489 {
490  const std::vector<T>& data = *m_data;
491  const int sizex = m_size.x;
492  // Look if we are inside the used space
493  if (x - 1 < m_size.x - 2 && y - 1 < m_size.y - 2 && z - 1 < m_size.z - 2) {
494 
495  // Lookup all neccessary Values
496  const T *H = &data[x + m_size.x * (y + m_size.y * z)];
497 
498  return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
499  Out((H[sizex] - H[-sizex]) * 0.5),
500  Out((H[m_xy] - H[-m_xy]) * 0.5));
501  }
502 
503  return T3DVector<Out>();
504 }
505 
506 
507 template <class T>
508 template <typename Out>
510 {
511  const int sizex = m_size.x;
512  // Lookup all neccessary Values
513  const T *H = &(*m_data)[hardcode];
514 
515  return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
516  Out((H[sizex] - H[-sizex]) * 0.5),
517  Out((H[m_xy] - H[-m_xy]) * 0.5));
518 }
519 
520 
524 template <>
525 template <typename Out>
527 {
528 
529  // Lookup all neccessary Values
530  return T3DVector<Out> (Out(((*m_data)[hardcode + 1] - (*m_data)[hardcode -1]) * 0.5),
531  Out(((*m_data)[hardcode + m_size.x] - (*m_data)[hardcode -m_size.x]) * 0.5),
532  Out(((*m_data)[hardcode + m_xy] - (*m_data)[hardcode -m_xy]) * 0.5));
533 }
534 
535 template <class T>
536 template <typename Out>
538 {
539  // This will become really funny
540  const int sizex = m_size.x;
541  const std::vector<T>& data = *m_data;
542  // Calculate the int coordinates near the POI
543  // and the distances
544  size_t x = size_t (p.x);
545  float dx = p.x - x;
546  float xm = 1 - dx;
547  size_t y = size_t (p.y);
548  float dy = p.y - y;
549  float ym = 1 - dy;
550  size_t z = size_t (p.z);
551  float dz = p.z - z;
552  float zm = 1 - dz;
553 
554  // Look if we are inside the used space
555  if (x-1 < m_size.x-3 && y -1 < m_size.y-3 && z - 1 < m_size.z-3 ) {
556  // Lookup all neccessary Values
557  const T *H000 = &data[x + sizex * y + m_xy * z];
558 
559  const T* H_100 = &H000[-m_xy];
560  const T* H_101 = &H_100[1];
561  const T* H_110 = &H_100[sizex];
562  const T* H_111 = &H_110[1];
563 
564  const T* H0_10 = &H000[-sizex];
565  const T* H0_11 = &H0_10[1];
566 
567  const T* H00_1 = &H000[-1];
568  const T* H001 = &H000[ 1];
569  const T* H002 = &H000[ 2];
570 
571 
572  const T* H010 = &H000[sizex];
573  const T* H011 = &H010[ 1];
574  const T* H012 = &H010[ 2];
575  const T* H01_1 = &H010[-1];
576 
577  const T* H020 = &H010[sizex];
578  const T* H021 = &H020[ 1];
579 
580  const T* H100 = &H000[m_xy];
581 
582  const T* H1_10 = &H100[sizex];
583  const T* H1_11 = &H1_10[1];
584 
585  const T* H10_1 = &H100[-1];
586  const T* H101 = &H100[ 1];
587  const T* H102 = &H100[ 2];
588 
589  const T* H110 = &H100[sizex];
590  const T* H111 = &H110[ 1];
591  const T* H112 = &H110[ 2];
592  const T* H11_1 = &H110[-1];
593 
594  const T* H120 = &H110[sizex];
595  const T* H121 = &H120[ 1];
596 
597  const T* H200 = &H100[m_xy];
598  const T* H201 = &H200[1];
599  const T* H210 = &H200[sizex];
600  const T* H211 = &H210[1];
601 
602  // use trilinear interpolation to calc the gradient
603  return T3DVector<Out> (
604  Out((zm * (ym * (dx * (*H002 - *H000) + xm * (*H001 - *H00_1))+
605  dy * (dx * (*H012 - *H010) + xm * (*H011 - *H01_1)))+
606  dz * (ym * (dx * (*H102 - *H100) + xm * (*H101 - *H10_1))+
607  dy * (dx * (*H112 - *H110) + xm * (*H111 - *H11_1)))) * 0.5),
608 
609  Out((zm * (ym * (xm * (*H010 - *H0_10) + dx * (*H011 - *H0_11))+
610  dy * (xm * (*H020 - *H000) + dx * (*H021 - *H001)))+
611  dz * (ym * (xm * (*H110 - *H1_10) + dx * (*H111 - *H1_11))+
612  dy * (xm * (*H120 - *H100) + dx * (*H121 - *H101)))) * 0.5),
613  Out((zm * (ym * (xm * (*H100 - *H_100) + dx * (*H101 - *H_101))+
614  dy * (xm * (*H110 - *H_110) + dx * (*H111 - *H_111)))+
615  dz * (ym * (xm * (*H200 - *H000) + dx * (*H201 - *H001))+
616  dy * (xm * (*H210 - *H010) + dx * (*H211 - *H011)))) * 0.5));
617  }
618  return T3DVector<Out>();
619 
620 }
621 
623 
624 #endif
T3DVector< Out > get_gradient(const T3DVector< float > &p) const
a 3D iterator that knows its position in the 3D grid ans supports iterating over sub-ranges ...
Definition: 3d/iterator.hh:189
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
A templated class of a 3D data field.
Definition: 3d/datafield.hh:45
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
T3DDatafield< unsigned char > C3DUBDatafield
a data field of float values
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:43
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
const_reference operator()(const C3DBounds &l) const
T3DDatafield< unsigned int > C3DUIDatafield
a data field of 32 bit unsigned int values
Generic type of a complex paramter.
Definition: parameter.hh:164
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
CTParameter< C3DFVector > C3DFVectorParameter
3D vector parameter type
const C3DBounds & get_size() const
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:52
const_iterator end() const
iterator begin()
T3DDatafield< int > C3DIDatafield
a data field of 32 bit signed int values
reference operator()(size_t x, size_t y, size_t z)
iterator end()
reference operator()(const C3DBounds &l)
const_iterator begin() const
T3DDatafield< unsigned long > C3DULDatafield
a data field of 32 bit unsigned int values
T3DDatafield< float > C3DFDatafield
a data field of float values
T x
vector element
Definition: 3d/vector.hh:51
bool holds_unique_data() const
Definition: 3d/datafield.hh:75
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:46