My Project
cstplan.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-2017 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 
22 #include <stdexcept>
23 #include <string>
24 #include <vector>
25 
26 #include <fftw3.h>
27 
28 #include <mia/core/msgstream.hh>
29 
31 
32 
33 template <typename T>
34 class TCSTPlan
35 {
36 public:
37  typedef typename T::value_type value_type;
38 
39  TCSTPlan(fftwf_r2r_kind forward, std::vector<int> size);
40 
41  virtual ~TCSTPlan();
42 
43  template <typename D>
44  void execute(const D& in_data, D& out_data) const;
45 
46  const std::vector<int>& get_size() const;
47 
48 private:
49  mutable value_type *m_in;
50  mutable value_type *m_out;
51 
52  std::vector<int> m_size;
53  float m_scale;
54  size_t m_n;
55 
56  fftwf_plan m_forward_plan;
57  fftwf_plan m_backward_plan;
58 
59  virtual void do_execute(value_type *buffer) const = 0;
60 
61  struct Scale {
62  Scale(float scale): m_scale(scale)
63  {
64  }
65  value_type operator ()(const value_type& x) const
66  {
67  return m_scale * x;
68  }
69  float m_scale;
70  };
71 };
72 
73 template <typename T>
74 const std::vector<int>& TCSTPlan<T>::get_size() const
75 {
76  return m_size;
77 }
78 
79 template <typename T>
80 TCSTPlan<T>::TCSTPlan(fftwf_r2r_kind forward, std::vector<int> size):
81  m_size(size),
82  m_scale(1.0)
83 {
84  std::string msg;
85  size_t rank = m_size.size();
86  fftwf_r2r_kind backward;
87 
88  switch (forward) {
89  case FFTW_RODFT00:
90  for (size_t i = 0; i < rank; ++i)
91  m_scale *= 0.5 / ( size[i] + 1.0);
92 
93  break;
94 
95  case FFTW_REDFT00:
96  for (size_t i = 0; i < rank; ++i)
97  m_scale *= 0.5 / ( size[i] - 1.0);
98 
99  break;
100 
101  default:
102  for (size_t i = 0; i < rank; ++i)
103  m_scale *= 0.5 / size[i];
104 
105  break;
106  }
107 
108  switch (forward) {
109  case FFTW_RODFT00:
110  case FFTW_REDFT00:
111  case FFTW_RODFT11:
112  case FFTW_REDFT11:
113  backward = forward;
114  break;
115 
116  case FFTW_RODFT10:
117  backward = FFTW_RODFT01;
118  break;
119 
120  case FFTW_REDFT10:
121  backward = FFTW_REDFT01;
122  break;
123 
124  case FFTW_RODFT01:
125  backward = FFTW_RODFT10;
126  break;
127 
128  case FFTW_REDFT01:
129  backward = FFTW_REDFT10;
130  break;
131 
132  default:
133  throw std::invalid_argument("TCSTPlan: unknown transformtion requested");
134  }
135 
136  std::vector<fftwf_r2r_kind> fw_kind(rank, forward);
137  std::vector<fftwf_r2r_kind> bw_kind(rank, backward);
138  const int howmany = sizeof(value_type) / sizeof(float);
139  cvdebug() << "howmany = " << howmany << "\n";
140  m_n = 1;
141 
142  for (size_t i = 0; i < rank; ++i)
143  m_n *= m_size[i];
144 
145  if (NULL == (m_in = (value_type *) fftwf_malloc(sizeof(T) * m_n))) {
146  msg = "unable to allocate FFTW data";
147  goto in_fail;
148  }
149 
150  if ( NULL == (m_out = (value_type *) fftwf_malloc(sizeof(T) * m_n))) {
151  msg = "unable to allocate FFTW data";
152  goto out_fail;
153  }
154 
155  if (0 == (m_forward_plan = fftwf_plan_many_r2r(rank, &size[0], howmany,
156  m_in, NULL, howmany, 1,
157  m_out, NULL, howmany, 1,
158  &fw_kind[0], FFTW_ESTIMATE))) {
159  msg = "unable to create FFTW forward plan";
160  goto plan_fw_fail;
161  }
162 
163  if (0 == (m_backward_plan = fftwf_plan_many_r2r(rank, &size[0], howmany,
164  m_out, NULL, howmany, 1,
165  m_in, NULL, howmany, 1,
166  &bw_kind[0], FFTW_ESTIMATE))) {
167  msg = "unable to create FFTW backward plan";
168  goto plan_bw_fail;
169  }
170 
171  return;
172 plan_bw_fail:
173  fftwf_destroy_plan(m_forward_plan);
174 plan_fw_fail:
175  fftwf_free(m_out);
176 out_fail:
177  fftwf_free(m_in);
178 in_fail:
179  throw std::runtime_error(msg);
180 }
181 
182 template <typename T>
184 {
185  fftwf_destroy_plan(m_backward_plan);
186  fftwf_destroy_plan(m_forward_plan);
187  fftwf_free(m_out);
188  fftwf_free(m_in);
189 }
190 
191 template <typename T>
192 template <typename D>
193 void TCSTPlan<T>::execute(const D& in_data, D& out_data) const
194 {
195  assert(m_n == in_data.size());
196  assert(m_n == out_data.size());
197  copy(in_data.begin(), in_data.end(), m_in);
198  fftwf_execute( m_forward_plan);
199  do_execute(m_out);
200  fftwf_execute(m_backward_plan);
201  transform( m_in, m_in + m_n, out_data.begin(), Scale(m_scale) );
202 }
203 
T::value_type value_type
Definition: cstplan.hh:37
void execute(const D &in_data, D &out_data) const
Definition: cstplan.hh:193
const std::vector< int > & get_size() const
Definition: cstplan.hh:74
TCSTPlan(fftwf_r2r_kind forward, std::vector< int > size)
Definition: cstplan.hh:80
virtual ~TCSTPlan()
Definition: cstplan.hh:183
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition: defines.hh:33
#define NS_MIA_END
conveniance define to end the mia namespace
Definition: defines.hh:36
CDebugSink & cvdebug()
Definition: msgstream.hh:226