3d/deformer.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 reg3d_deformer_hh
22 #define reg3d_deformer_hh
23 
24 #include <memory>
25 #include <tbb/parallel_for.h>
26 #include <tbb/blocked_range.h>
27 
28 
29 #include <mia/3d/image.hh>
30 #include <mia/3d/filter.hh>
31 #include <mia/3d/interpolator.hh>
32 #include <mia/core/threadedmsg.hh>
33 
34 
36 
43 struct FDeformer3D: public TFilter<P3DImage> {
45  m_vf(vf),
46  m_ipfac(ipfac)
47  {
48  }
49  template <typename T>
50  P3DImage operator () (const T3DImage<T>& image) const {
51  T3DImage<T> *timage = new T3DImage<T>(m_vf.get_size(), image);
52  P3DImage result(timage);
53  this->operator()(image, *timage);
54  return result;
55  }
56 
57  template <typename T>
58  P3DImage operator () (const T3DImage<T>& image, T3DImage<T>& result) const {
59  assert(result.get_size() == m_vf.get_size());
60  std::unique_ptr<T3DConvoluteInterpolator<T> > interp(m_ipfac.create(image.data()));
61  const auto& rinterp = *interp;
62 
63  auto callback = [this, &rinterp, &result](const tbb::blocked_range<size_t>& range){
64  CThreadMsgStream thread_stream;
65  CWeightCache cache = rinterp.create_cache();
66  for (auto z = range.begin(); z != range.end();++z) {
67  auto r = result.begin_at(0,0,z);
68  auto v = m_vf.begin_at(0,0,z);
69  for (size_t y = 0; y < result.get_size().y; ++y)
70  for (size_t x = 0; x < result.get_size().x; ++x, ++r, ++v)
71  *r = rinterp(C3DFVector(x - v->x, y - v->y, z - v->z), cache);
72  }
73  };
74  tbb::parallel_for(tbb::blocked_range<size_t>(0, result.get_size().z, 1), callback);
75  return P3DImage();
76  }
77 
78 private:
79  C3DFVectorfield m_vf;
80  C3DInterpolatorFactory m_ipfac;
81 };
82 
83 
85 
86 #endif