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-2013 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::shared_ptr<std::vector<T> > ref_data_type;
48 
50  C3DBounds m_size;
51 
53  size_t m_xy;
54 
56  ref_data_type m_data;
57 
59  static const T Zero;
60 
61  static const unsigned int m_elements;
62 
63 public:
64 
67  void make_single_ref();
68 
73  bool holds_unique_data()const {
74  return m_data.unique();
75  }
76 
77 
79 
81  typedef typename std::vector<T>::iterator iterator;
82  typedef typename std::vector<T>::const_iterator const_iterator;
83  typedef typename std::vector<T>::const_reference const_reference;
84  typedef typename std::vector<T>::reference reference;
85  typedef typename std::vector<T>::const_pointer const_pointer;
86  typedef typename std::vector<T>::pointer pointer;
87  typedef typename std::vector<T>::value_type value_type;
88  typedef typename std::vector<T>::size_type size_type;
89  typedef typename std::vector<T>::difference_type difference_type;
90  typedef typename atomic_data<T>::type atomic_type;
91  typedef range3d_iterator<iterator> range_iterator;
92  typedef range3d_iterator<const_iterator> const_range_iterator;
93  typedef C3DBounds dimsize_type;
95 
96  T3DDatafield();
97 
99  T3DDatafield(const C3DBounds& _Size);
100 
105  T3DDatafield(const C3DBounds& size, const T *data);
106 
107 
109  T3DDatafield(const T3DDatafield<T>& org);
110 
112  virtual ~T3DDatafield();
113 
118  template <typename Out>
120 
122  template <typename Out>
123  T3DVector<Out> get_gradient(size_t x, size_t y, size_t z) const;
124 
126  template <typename Out>
127  T3DVector<Out> get_gradient(int index) const;
128 
130  value_type get_interpol_val_at(const T3DVector<float >& p) const;
131 
132  /* some rough interpolation using barycentric coordinates, needs less addition and
133  multiplications then tri-linear interp. but is usally of low quality
134  \remark this function may vanish
135  value_type get_barycent_interpol_val_at(const T3DVector<float >& p) const;
136  */
137 
139  value_type get_trilin_interpol_val_at(const T3DVector<float >& p) const;
140 
143  value_type get_block_avrg(const C3DBounds& Start, const C3DBounds& BlockSize) const;
144 
149  T3DDatafield& operator = (const T3DDatafield& org);
150 
152  const C3DBounds& get_size() const
153  {
154  return m_size;
155  }
156 
158  void clear();
159 
161  size_type size()const
162  {
163  return m_data->size();
164  }
165 
167  void swap(T3DDatafield& other);
168 
170  value_type get_avg();
171 
174  value_type strip_avg();
175 
176 
178  value_type operator()(const T3DVector<float >& pos)const;
179 
181  const_reference operator()(size_t x, size_t y, size_t z) const
182  {
183  // Look if we are inside, and give reference, else give the zero
184  if (x < m_size.x && y < m_size.y && z < m_size.z) {
185  const std::vector<T>& data = *m_data;
186  return data[x+ m_size.x * (y + m_size.y * z)];
187  }
188  return Zero;
189  }
190 
191 
193  const_reference operator()(const C3DBounds& l)const
194  {
195  return (*this)(l.x,l.y,l.z);
196  }
197 
199  reference operator()(size_t x, size_t y, size_t z)
200  {
201  // Look if we are inside, and give reference, else throw exception
202  // since write access is wanted
203  assert(x < m_size.x && y < m_size.y && z < m_size.z);
204  return (*m_data)[x + m_size.x *(y + m_size.y * z)];
205  }
206 
207 
208 
210  reference operator()(const C3DBounds& l)
211  {
212  return (*this)(l.x,l.y,l.z);
213  }
214 
216  void get_data_line_x(int y, int z, std::vector<T>& buffer)const;
217 
219  void get_data_line_y(int x, int z, std::vector<T>& buffer)const;
220 
222  void get_data_line_z(int x, int y, std::vector<T>& buffer)const;
223 
225  void put_data_line_x(int y, int z, const std::vector<T> &buffer);
226 
228  void put_data_line_y(int x, int z, const std::vector<T> &buffer);
229 
231  void put_data_line_z(int x, int y, const std::vector<T> &buffer);
232 
234  template <class TMask>
235  void mask(const TMask& m);
236 
250  void read_xslice_flat(size_t x, std::vector<atomic_type>& buffer) const;
251 
264  void read_yslice_flat(size_t y, std::vector<atomic_type>& buffer) const;
265 
278  void read_zslice_flat(size_t z, std::vector<atomic_type>& buffer) const;
279 
284  void write_zslice_flat(size_t z, const std::vector<atomic_type>& buffer);
285 
286 
291  void write_yslice_flat(size_t y, const std::vector<atomic_type>& buffer);
292 
297  void write_xslice_flat(size_t x, const std::vector<atomic_type>& buffer);
298 
304  T2DDatafield<T> get_data_plane_xy(size_t z)const;
305 
311  T2DDatafield<T> get_data_plane_yz(size_t x)const;
312 
318  T2DDatafield<T> get_data_plane_xz(size_t y)const;
319 
325  void put_data_plane_xy(size_t z, const T2DDatafield<T>& p);
326 
332  void put_data_plane_yz(size_t x, const T2DDatafield<T>& p);
333 
339  void put_data_plane_xz(size_t y, const T2DDatafield<T>& p);
340 
342  const_iterator begin()const
343  {
344  return m_data->begin();
345  }
346 
350  const_iterator begin_at(size_t x, size_t y, size_t z)const
351  {
352  return m_data->begin() + (z * m_size.y + y) * m_size.x + x;
353  }
354 
355 
359  const_iterator end()const
360  {
361  return m_data->end();
362  }
363 
367  iterator begin()
368  {
369  make_single_ref();
370  return m_data->begin();
371  }
372 
375  range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end);
376 
378  range_iterator end_range(const C3DBounds& begin, const C3DBounds& end);
379 
380 
383  const_range_iterator begin_range(const C3DBounds& begin, const C3DBounds& end)const;
384 
386  const_range_iterator end_range(const C3DBounds& begin, const C3DBounds& end)const;
387 
388 
397  iterator begin_at(size_t x, size_t y, size_t z)
398  {
399  make_single_ref();
400  return m_data->begin() + (z * m_size.y + y) * m_size.x + x;
401  }
402 
406  iterator end()
407  {
408  make_single_ref();
409  return m_data->end();
410  }
411 
413  const_reference operator[](int i)const
414  {
415  return (*m_data)[i];
416  }
417 
421  reference operator[](int i)
422  {
423  assert(m_data.unique());
424  return (*m_data)[i];
425  }
426 
427 
429  size_t get_plane_size_xy()const
430  {
431  return m_xy;
432  };
433 
434 private:
435 };
436 
437 
440 
443 
446 
447 
450 
453 
456 
459 
462 
465 
467 DECLARE_TYPE_DESCR(C3DBounds);
468 DECLARE_TYPE_DESCR(C3DFVector);
470 
471 // some implementations
472 
473 template <class T>
474 template <typename Out>
475 T3DVector<Out> T3DDatafield<T>::get_gradient(size_t x, size_t y, size_t z) const
476 {
477  const std::vector<T>& data = *m_data;
478  const int sizex = m_size.x;
479  // Look if we are inside the used space
480  if (x - 1 < m_size.x - 2 && y - 1 < m_size.y - 2 && z - 1 < m_size.z - 2) {
481 
482  // Lookup all neccessary Values
483  const T *H = &data[x + m_size.x * (y + m_size.y * z)];
484 
485  return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
486  Out((H[sizex] - H[-sizex]) * 0.5),
487  Out((H[m_xy] - H[-m_xy]) * 0.5));
488  }
489 
490  return T3DVector<Out>();
491 }
492 
493 
494 template <class T>
495 template <typename Out>
497 {
498  const int sizex = m_size.x;
499  // Lookup all neccessary Values
500  const T *H = &(*m_data)[hardcode];
501 
502  return T3DVector<Out> (Out((H[1] - H[-1]) * 0.5),
503  Out((H[sizex] - H[-sizex]) * 0.5),
504  Out((H[m_xy] - H[-m_xy]) * 0.5));
505 }
506 
507 
511 template <>
512 template <typename Out>
514 {
515 
516  // Lookup all neccessary Values
517  return T3DVector<Out> (Out(((*m_data)[hardcode + 1] - (*m_data)[hardcode -1]) * 0.5),
518  Out(((*m_data)[hardcode + m_size.x] - (*m_data)[hardcode -m_size.x]) * 0.5),
519  Out(((*m_data)[hardcode + m_xy] - (*m_data)[hardcode -m_xy]) * 0.5));
520 }
521 
522 template <class T>
523 template <typename Out>
525 {
526  // This will become really funny
527  const int sizex = m_size.x;
528  const std::vector<T>& data = *m_data;
529  // Calculate the int coordinates near the POI
530  // and the distances
531  size_t x = size_t (p.x);
532  float dx = p.x - x;
533  float xm = 1 - dx;
534  size_t y = size_t (p.y);
535  float dy = p.y - y;
536  float ym = 1 - dy;
537  size_t z = size_t (p.z);
538  float dz = p.z - z;
539  float zm = 1 - dz;
540 
541  // Look if we are inside the used space
542  if (x-1 < m_size.x-3 && y -1 < m_size.y-3 && z - 1 < m_size.z-3 ) {
543  // Lookup all neccessary Values
544  const T *H000 = &data[x + sizex * y + m_xy * z];
545 
546  const T* H_100 = &H000[-m_xy];
547  const T* H_101 = &H_100[1];
548  const T* H_110 = &H_100[sizex];
549  const T* H_111 = &H_110[1];
550 
551  const T* H0_10 = &H000[-sizex];
552  const T* H0_11 = &H0_10[1];
553 
554  const T* H00_1 = &H000[-1];
555  const T* H001 = &H000[ 1];
556  const T* H002 = &H000[ 2];
557 
558 
559  const T* H010 = &H000[sizex];
560  const T* H011 = &H010[ 1];
561  const T* H012 = &H010[ 2];
562  const T* H01_1 = &H010[-1];
563 
564  const T* H020 = &H010[sizex];
565  const T* H021 = &H020[ 1];
566 
567  const T* H100 = &H000[m_xy];
568 
569  const T* H1_10 = &H100[sizex];
570  const T* H1_11 = &H1_10[1];
571 
572  const T* H10_1 = &H100[-1];
573  const T* H101 = &H100[ 1];
574  const T* H102 = &H100[ 2];
575 
576  const T* H110 = &H100[sizex];
577  const T* H111 = &H110[ 1];
578  const T* H112 = &H110[ 2];
579  const T* H11_1 = &H110[-1];
580 
581  const T* H120 = &H110[sizex];
582  const T* H121 = &H120[ 1];
583 
584  const T* H200 = &H100[m_xy];
585  const T* H201 = &H200[1];
586  const T* H210 = &H200[sizex];
587  const T* H211 = &H210[1];
588 
589  // use trilinear interpolation to calc the gradient
590  return T3DVector<Out> (
591  Out((zm * (ym * (dx * (*H002 - *H000) + xm * (*H001 - *H00_1))+
592  dy * (dx * (*H012 - *H010) + xm * (*H011 - *H01_1)))+
593  dz * (ym * (dx * (*H102 - *H100) + xm * (*H101 - *H10_1))+
594  dy * (dx * (*H112 - *H110) + xm * (*H111 - *H11_1)))) * 0.5),
595 
596  Out((zm * (ym * (xm * (*H010 - *H0_10) + dx * (*H011 - *H0_11))+
597  dy * (xm * (*H020 - *H000) + dx * (*H021 - *H001)))+
598  dz * (ym * (xm * (*H110 - *H1_10) + dx * (*H111 - *H1_11))+
599  dy * (xm * (*H120 - *H100) + dx * (*H121 - *H101)))) * 0.5),
600  Out((zm * (ym * (xm * (*H100 - *H_100) + dx * (*H101 - *H_101))+
601  dy * (xm * (*H110 - *H_110) + dx * (*H111 - *H_111)))+
602  dz * (ym * (xm * (*H200 - *H000) + dx * (*H201 - *H001))+
603  dy * (xm * (*H210 - *H010) + dx * (*H211 - *H011)))) * 0.5));
604  }
605  return T3DVector<Out>();
606 
607 }
608 
610 
611 #endif