cstplan.hh
Go to the documentation of this file.
/* -*- mia-c++ -*-
*
* This file is part of MIA - a toolbox for medical image analysis
* Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
*
* MIA is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MIA; if not, see <http://www.gnu.org/licenses/>.
*
*/
#include <stdexcept>
#include <string>
#include <vector>
#include <fftw3.h>
#include <mia/core/msgstream.hh>
NS_MIA_BEGIN
template <typename T>
class TCSTPlan
{
public:
typedef typename T::value_type value_type;
TCSTPlan(fftwf_r2r_kind forward, std::vector<int> size);
virtual ~TCSTPlan();
template <typename D>
void execute(const D& in_data, D& out_data) const;
const std::vector<int>& get_size() const;
private:
mutable value_type *m_in;
mutable value_type *m_out;
std::vector<int> m_size;
float m_scale;
size_t m_n;
fftwf_plan m_forward_plan;
fftwf_plan m_backward_plan;
virtual void do_execute(value_type *buffer) const = 0;
struct Scale {
Scale(float scale): m_scale(scale)
{
}
value_type operator ()(const value_type& x) const
{
return m_scale * x;
}
float m_scale;
};
};
template <typename T>
const std::vector<int>& TCSTPlan<T>::get_size() const
{
return m_size;
}
template <typename T>
TCSTPlan<T>::TCSTPlan(fftwf_r2r_kind forward, std::vector<int> size):
m_size(size),
m_scale(1.0)
{
std::string msg;
size_t rank = m_size.size();
fftwf_r2r_kind backward;
switch (forward) {
case FFTW_RODFT00:
for (size_t i = 0; i < rank; ++i)
m_scale *= 0.5 / ( size[i] + 1.0);
break;
case FFTW_REDFT00:
for (size_t i = 0; i < rank; ++i)
m_scale *= 0.5 / ( size[i] - 1.0);
break;
default:
for (size_t i = 0; i < rank; ++i)
m_scale *= 0.5 / size[i];
break;
}
switch (forward) {
case FFTW_RODFT00:
case FFTW_REDFT00:
case FFTW_RODFT11:
case FFTW_REDFT11:
backward = forward;
break;
case FFTW_RODFT10:
backward = FFTW_RODFT01;
break;
case FFTW_REDFT10:
backward = FFTW_REDFT01;
break;
case FFTW_RODFT01:
backward = FFTW_RODFT10;
break;
case FFTW_REDFT01:
backward = FFTW_REDFT10;
break;
default:
throw std::invalid_argument("TCSTPlan: unknown transformtion requested");
}
std::vector<fftwf_r2r_kind> fw_kind(rank, forward);
std::vector<fftwf_r2r_kind> bw_kind(rank, backward);
const int howmany = sizeof(value_type) / sizeof(float);
cvdebug() << "howmany = " << howmany << "\n";
m_n = 1;
for (size_t i = 0; i < rank; ++i)
m_n *= m_size[i];
if (NULL == (m_in = (value_type *) fftwf_malloc(sizeof(T) * m_n))) {
msg = "unable to allocate FFTW data";
goto in_fail;
}
if ( NULL == (m_out = (value_type *) fftwf_malloc(sizeof(T) * m_n))) {
msg = "unable to allocate FFTW data";
goto out_fail;
}
if (0 == (m_forward_plan = fftwf_plan_many_r2r(rank, &size[0], howmany,
m_in, NULL, howmany, 1,
m_out, NULL, howmany, 1,
&fw_kind[0], FFTW_ESTIMATE))) {
msg = "unable to create FFTW forward plan";
goto plan_fw_fail;
}
if (0 == (m_backward_plan = fftwf_plan_many_r2r(rank, &size[0], howmany,
m_out, NULL, howmany, 1,
m_in, NULL, howmany, 1,
&bw_kind[0], FFTW_ESTIMATE))) {
msg = "unable to create FFTW backward plan";
goto plan_bw_fail;
}
return;
plan_bw_fail:
fftwf_destroy_plan(m_forward_plan);
plan_fw_fail:
fftwf_free(m_out);
out_fail:
fftwf_free(m_in);
in_fail:
throw std::runtime_error(msg);
}
template <typename T>
TCSTPlan<T>::~TCSTPlan()
{
fftwf_destroy_plan(m_backward_plan);
fftwf_destroy_plan(m_forward_plan);
fftwf_free(m_out);
fftwf_free(m_in);
}
template <typename T>
template <typename D>
void TCSTPlan<T>::execute(const D& in_data, D& out_data) const
{
assert(m_n == in_data.size());
assert(m_n == out_data.size());
copy(in_data.begin(), in_data.end(), m_in);
fftwf_execute( m_forward_plan);
do_execute(m_out);
fftwf_execute(m_backward_plan);
transform( m_in, m_in + m_n, out_data.begin(), Scale(m_scale) );
}
NS_MIA_END
TCSTPlan::value_type
T::value_type value_type
Definition cstplan.hh:37
TCSTPlan::execute
void execute(const D &in_data, D &out_data) const
Definition cstplan.hh:193
TCSTPlan::get_size
const std::vector< int > & get_size() const
Definition cstplan.hh:74
TCSTPlan::TCSTPlan
TCSTPlan(fftwf_r2r_kind forward, std::vector< int > size)
Definition cstplan.hh:80
TCSTPlan::~TCSTPlan
virtual ~TCSTPlan()
Definition cstplan.hh:183
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition defines.hh:33
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition defines.hh:36
msgstream.hh
cvdebug
CDebugSink & cvdebug()
Definition msgstream.hh:226