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
33template <typename T>
35{
36public:
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
48private:
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
73template <typename T>
74const std::vector<int>& TCSTPlan<T>::get_size() const
75{
76 return m_size;
77}
78
79template <typename T>
80TCSTPlan<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;
172plan_bw_fail:
173 fftwf_destroy_plan(m_forward_plan);
174plan_fw_fail:
175 fftwf_free(m_out);
176out_fail:
177 fftwf_free(m_in);
178in_fail:
179 throw std::runtime_error(msg);
180}
181
182template <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
191template <typename T>
192template <typename D>
193void 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