diff --git a/include/IRFFT.h b/include/IRFFT.h index e2313c1..aae4e9a 100644 --- a/include/IRFFT.h +++ b/include/IRFFT.h @@ -31,10 +31,11 @@ #define FFT_HXX_ // FFTW includes: -#include +// #include // ITK includes: #include +#include // system includes: #include @@ -44,11 +45,12 @@ namespace itk_fft { + using TransformDirectionEnum = itk::ComplexToComplexFFTImageFilterEnums::TransformDirection; //---------------------------------------------------------------- // set_num_fftw_threads // // set's number of threads used by fftw, returns previous value: - extern std::size_t set_num_fftw_threads(std::size_t num_threads); + // extern std::size_t set_num_fftw_threads(std::size_t num_threads); //---------------------------------------------------------------- // itk_image_t @@ -65,13 +67,23 @@ namespace itk_fft // typedef std::complex fft_complex_t; + //---------------------------------------------------------------- + // itk_complex_image_t + // + typedef itk::Image itk_complex_image_t; + + //---------------------------------------------------------------- + // itk_complex_imageptr_t + // + typedef itk_complex_image_t::Pointer itk_complex_imageptr_t; + //---------------------------------------------------------------- // fft_data_t // class fft_data_t { public: - fft_data_t(): data_(NULL), nx_(0), ny_(0) {} + fft_data_t(): image_(), nx_(0), ny_(0) {} fft_data_t(const unsigned int w, const unsigned int h); explicit fft_data_t(const itk_imageptr_t & real); fft_data_t(const itk_imageptr_t & real, @@ -105,8 +117,8 @@ namespace itk_fft void apply_lp_filter(const double r, const double s = 0); // accessors: - inline const fft_complex_t * data() const { return data_; } - inline fft_complex_t * data() { return data_; } + inline const itk_complex_image_t * data() const { return image_.get(); } + inline itk_complex_image_t * data() { return image_.get(); } inline unsigned int nx() const { return nx_; } inline unsigned int ny() const { return ny_; } @@ -121,14 +133,14 @@ namespace itk_fft inline const fft_complex_t & at(const unsigned int & x, const unsigned int & y) const - { return data_[y + ny_ * x]; } + { return image_->GetPixel({ x, y }); } inline fft_complex_t & at(const unsigned int & x, const unsigned int & y) - { return data_[y + ny_ * x]; } + { return image_->GetPixel({ x, y }); } protected: - fft_complex_t * data_; + itk_complex_imageptr_t image_; unsigned int nx_; unsigned int ny_; }; @@ -144,7 +156,7 @@ namespace itk_fft // fft // extern bool - fft(const fft_data_t & in, fft_data_t & out, int sign = FFTW_FORWARD); + fft(const fft_data_t & in, fft_data_t & out, TransformDirectionEnum sign = TransformDirectionEnum::FORWARD); //---------------------------------------------------------------- // ifft @@ -168,239 +180,239 @@ namespace itk_fft } - //---------------------------------------------------------------- - // fn_fft_c_t - // - typedef fft_complex_t(*fn_fft_c_t)(const fft_complex_t & in); - - //---------------------------------------------------------------- - // fn_fft_cc_t - // - typedef fft_complex_t(*fn_fft_cc_t)(const fft_complex_t & a, - const fft_complex_t & b); - - //---------------------------------------------------------------- - // elem_by_elem - // - extern void - elem_by_elem(fn_fft_c_t f, - const fft_data_t & in, - fft_data_t & out); - - //---------------------------------------------------------------- - // elem_by_elem - // - extern void - elem_by_elem(fft_complex_t(*f)(const float & a, - const fft_complex_t & b), - const float & a, - const fft_data_t & b, - fft_data_t & c); - - //---------------------------------------------------------------- - // elem_by_elem - // - extern void - elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, - const float & b), - const fft_data_t & a, - const float & b, - fft_data_t & c); - - //---------------------------------------------------------------- - // elem_by_elem - // - extern void - elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, - const fft_complex_t & b), - const fft_complex_t & a, - const fft_data_t & b, - fft_data_t & c); - - //---------------------------------------------------------------- - // elem_by_elem - // - extern void - elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, - const fft_complex_t & b), - const fft_data_t & a, - const fft_complex_t & b, - fft_data_t & c); - - //---------------------------------------------------------------- - // elem_by_elem - // - extern void - elem_by_elem(fn_fft_cc_t f, - const fft_data_t & a, - const fft_data_t & b, - fft_data_t & c); - - - //---------------------------------------------------------------- - // conj - // - // element-by-element complex conjugate: - // - inline void - conj(const fft_data_t & in, fft_data_t & out) - { elem_by_elem(std::conj, in, out); } - - //---------------------------------------------------------------- - // conj - // - inline fft_data_t - conj(const fft_data_t & in) - { - fft_data_t out; - conj(in, out); - return out; - } - - - //---------------------------------------------------------------- - // sqrt - // - // element-by-element square root: - // - inline void - sqrt(const fft_data_t & in, fft_data_t & out) - { elem_by_elem(std::sqrt, in, out); } - - //---------------------------------------------------------------- - // sqrt - // - inline fft_data_t - sqrt(const fft_data_t & in) - { - fft_data_t out; - sqrt(in, out); - return out; - } - - - //---------------------------------------------------------------- - // _mul - // - template - inline fft_complex_t - _mul(const a_t & a, const b_t & b) - { return a * b; } - - //---------------------------------------------------------------- - // mul - // - // element-by-element multiplication, c = a * b: - // - template - inline void - mul(const a_t & a, const b_t & b, fft_data_t & c) - { elem_by_elem(_mul, a, b, c); } - - //---------------------------------------------------------------- - // mul - // - template - inline fft_data_t - mul(const a_t & a, const b_t & b) - { - fft_data_t c; - mul(a, b, c); - return c; - } - - - //---------------------------------------------------------------- - // _div - // - template - inline fft_complex_t - _div(const a_t & a, const b_t & b) - { return a / b; } - - //---------------------------------------------------------------- - // div - // - // element-by-element division, c = a / b: - // - template - inline void - div(const a_t & a, const b_t & b, fft_data_t & c) - { elem_by_elem(_div, a, b, c); } - - //---------------------------------------------------------------- - // div - // - template - inline fft_data_t - div(const a_t & a, const b_t & b) - { - fft_data_t c; - div(a, b, c); - return c; - } - - - //---------------------------------------------------------------- - // _add - // - template - inline fft_complex_t - _add(const a_t & a, const b_t & b) - { return a + b; } - - //---------------------------------------------------------------- - // add - // - // element-by-element addition, c = a + b: - // - template - inline void - add(const a_t & a, const b_t & b, fft_data_t & c) - { elem_by_elem(_add, a, b, c); } - - //---------------------------------------------------------------- - // add - // - template - inline fft_data_t - add(const a_t & a, const b_t & b) - { - fft_data_t c; - add(a, b, c); - return c; - } - - - //---------------------------------------------------------------- - // _sub - // - template - inline fft_complex_t - _sub(const a_t & a, const b_t & b) - { return a - b; } - - //---------------------------------------------------------------- - // sub - // - // element-by-element subtraction, c = a - b: - // - template - inline void - sub(const a_t & a, const b_t & b, fft_data_t & c) - { elem_by_elem(_sub, a, b, c); } - - //---------------------------------------------------------------- - // sub - // - template - inline fft_data_t - sub(const a_t & a, const b_t & b) - { - fft_data_t c; - sub(a, b, c); - return c; - } + // //---------------------------------------------------------------- + // // fn_fft_c_t + // // + // typedef fft_complex_t(*fn_fft_c_t)(const fft_complex_t & in); + + // //---------------------------------------------------------------- + // // fn_fft_cc_t + // // + // typedef fft_complex_t(*fn_fft_cc_t)(const fft_complex_t & a, + // const fft_complex_t & b); + + // //---------------------------------------------------------------- + // // elem_by_elem + // // + // extern void + // elem_by_elem(fn_fft_c_t f, + // const fft_data_t & in, + // fft_data_t & out); + + // //---------------------------------------------------------------- + // // elem_by_elem + // // + // extern void + // elem_by_elem(fft_complex_t(*f)(const float & a, + // const fft_complex_t & b), + // const float & a, + // const fft_data_t & b, + // fft_data_t & c); + + // //---------------------------------------------------------------- + // // elem_by_elem + // // + // extern void + // elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + // const float & b), + // const fft_data_t & a, + // const float & b, + // fft_data_t & c); + + // //---------------------------------------------------------------- + // // elem_by_elem + // // + // extern void + // elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + // const fft_complex_t & b), + // const fft_complex_t & a, + // const fft_data_t & b, + // fft_data_t & c); + + // //---------------------------------------------------------------- + // // elem_by_elem + // // + // extern void + // elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + // const fft_complex_t & b), + // const fft_data_t & a, + // const fft_complex_t & b, + // fft_data_t & c); + + // //---------------------------------------------------------------- + // // elem_by_elem + // // + // extern void + // elem_by_elem(fn_fft_cc_t f, + // const fft_data_t & a, + // const fft_data_t & b, + // fft_data_t & c); + + + // //---------------------------------------------------------------- + // // conj + // // + // // element-by-element complex conjugate: + // // + // inline void + // conj(const fft_data_t & in, fft_data_t & out) + // { elem_by_elem(std::conj, in, out); } + + // //---------------------------------------------------------------- + // // conj + // // + // inline fft_data_t + // conj(const fft_data_t & in) + // { + // fft_data_t out; + // conj(in, out); + // return out; + // } + + + // //---------------------------------------------------------------- + // // sqrt + // // + // // element-by-element square root: + // // + // inline void + // sqrt(const fft_data_t & in, fft_data_t & out) + // { elem_by_elem(std::sqrt, in, out); } + + // //---------------------------------------------------------------- + // // sqrt + // // + // inline fft_data_t + // sqrt(const fft_data_t & in) + // { + // fft_data_t out; + // sqrt(in, out); + // return out; + // } + + + // //---------------------------------------------------------------- + // // _mul + // // + // template + // inline fft_complex_t + // _mul(const a_t & a, const b_t & b) + // { return a * b; } + + // //---------------------------------------------------------------- + // // mul + // // + // // element-by-element multiplication, c = a * b: + // // + // template + // inline void + // mul(const a_t & a, const b_t & b, fft_data_t & c) + // { elem_by_elem(_mul, a, b, c); } + + // //---------------------------------------------------------------- + // // mul + // // + // template + // inline fft_data_t + // mul(const a_t & a, const b_t & b) + // { + // fft_data_t c; + // mul(a, b, c); + // return c; + // } + + + // //---------------------------------------------------------------- + // // _div + // // + // template + // inline fft_complex_t + // _div(const a_t & a, const b_t & b) + // { return a / b; } + + // //---------------------------------------------------------------- + // // div + // // + // // element-by-element division, c = a / b: + // // + // template + // inline void + // div(const a_t & a, const b_t & b, fft_data_t & c) + // { elem_by_elem(_div, a, b, c); } + + // //---------------------------------------------------------------- + // // div + // // + // template + // inline fft_data_t + // div(const a_t & a, const b_t & b) + // { + // fft_data_t c; + // div(a, b, c); + // return c; + // } + + + // //---------------------------------------------------------------- + // // _add + // // + // template + // inline fft_complex_t + // _add(const a_t & a, const b_t & b) + // { return a + b; } + + // //---------------------------------------------------------------- + // // add + // // + // // element-by-element addition, c = a + b: + // // + // template + // inline void + // add(const a_t & a, const b_t & b, fft_data_t & c) + // { elem_by_elem(_add, a, b, c); } + + // //---------------------------------------------------------------- + // // add + // // + // template + // inline fft_data_t + // add(const a_t & a, const b_t & b) + // { + // fft_data_t c; + // add(a, b, c); + // return c; + // } + + + // //---------------------------------------------------------------- + // // _sub + // // + // template + // inline fft_complex_t + // _sub(const a_t & a, const b_t & b) + // { return a - b; } + + // //---------------------------------------------------------------- + // // sub + // // + // // element-by-element subtraction, c = a - b: + // // + // template + // inline void + // sub(const a_t & a, const b_t & b, fft_data_t & c) + // { elem_by_elem(_sub, a, b, c); } + + // //---------------------------------------------------------------- + // // sub + // // + // template + // inline fft_data_t + // sub(const a_t & a, const b_t & b) + // { + // fft_data_t c; + // sub(a, b, c); + // return c; + // } } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e46209b..95b69a8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -20,6 +20,7 @@ set(Nornir_SRCS IRAABBox.cxx IRTerminator.cxx IRStdMutex.cxx + IRFFT.cxx ) itk_module_add_library(Nornir ${Nornir_SRCS}) diff --git a/src/IRFFT.cxx b/src/IRFFT.cxx new file mode 100644 index 0000000..79dae26 --- /dev/null +++ b/src/IRFFT.cxx @@ -0,0 +1,693 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* +This program 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 2 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 this program; if not, write to the Free Software +Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : fft.cxx +// Author : Pavel A. Koshevoy +// Created : 2005/06/03 10:16 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Wrapper class and helper functions for working with +// FFTW3 and ITK images. + +// system includes: +#ifndef _USE_MATH_DEFINES +#define _USE_MATH_DEFINES +#endif +#include +#include + +// ITK includes: +#include +#include +#include +#include +#include +#include + +// local includes: +#include "IRFFT.h" +#include "IRMutexInterface.h" +#include "itkIRUtils.h" + + +namespace itk_fft +{ + // static the_mutex_interface_t * mutex = NULL; + // //---------------------------------------------------------------- + // // fftw_mutex + // // + // static the_mutex_interface_t * + // fftw_mutex() + // { + // if(mutex == NULL) + // mutex = the_mutex_interface_t::create(); + // return mutex; + // } + + //---------------------------------------------------------------- + // NUM_FFTW_THREADS + // + static std::size_t NUM_FFTW_THREADS = 1; + + //---------------------------------------------------------------- + // set_num_fftw_threads + // + // set's number of threads used by fftw, returns previous value: + // std::size_t set_num_fftw_threads(std::size_t num_threads) + // { + // the_lock_t lock(fftw_mutex()); + // std::size_t prev = NUM_FFTW_THREADS; + // if (num_threads > 0) + // { + + // NUM_FFTW_THREADS = num_threads; + // } + + // return prev; + // } + + + // //---------------------------------------------------------------- + // // fft_cache_t + // // + // class fft_cache_t + // { + // public: + // fft_cache_t(): + // fwd_(NULL), + // inv_(NULL), + // w_(0), + // h_(0), + // h_complex_(0), + // h_padded_(0), + // buffer_(NULL) + // {} + + // ~fft_cache_t() + // { + // clear(); + // } + + // void clear() + // { + // if (buffer_) + // { + // // fftw is not thread safe: + // // the_lock_t lock(fftw_mutex()); + + // fftwf_destroy_plan(fwd_); + // fwd_ = NULL; + + // fftwf_destroy_plan(inv_); + // inv_ = NULL; + + // fftwf_free(buffer_); + // buffer_ = NULL; + // } + + // w_ = 0; + // h_ = 0; + // h_complex_ = 0; + // h_padded_ = 0; + // } + + // void update(const unsigned int & w, const unsigned int & h) + // { + // if (w_ == w && h_ == h) + // { + // return; + // } + + // // fftw is not thread safe: + // // the_lock_t lock(fftw_mutex()); + + // // if (buffer_) + // // { + // // fftwf_destroy_plan(fwd_); + // // fftwf_destroy_plan(inv_); + // // fftwf_free(buffer_); + // // } + + // w_ = w; + // h_ = h; + + // h_complex_ = h_ / 2 + 1; + // h_padded_ = h_complex_ * 2; + + // buffer_ = (float *)(fftwf_malloc(w_ * h_padded_ * sizeof(float))); + // fftwf_plan_with_nthreads(NUM_FFTW_THREADS); + // fwd_ = fftwf_plan_dft_r2c_2d(w_, + // h_, + // buffer_, + // (fftwf_complex *)buffer_, + // FFTW_DESTROY_INPUT | FFTW_ESTIMATE); + + // fft_data_t dummy_in(4, 4); + // fft_data_t dummy_out(4, 4); + // fftwf_plan_with_nthreads(NUM_FFTW_THREADS); + // inv_ = fftwf_plan_dft_2d(w_, + // h_, + // (fftwf_complex *)(dummy_in.data()), + // (fftwf_complex *)(dummy_out.data()), + // FFTW_BACKWARD, + // FFTW_ESTIMATE); + // } + + // fftwf_plan fwd_; // analysis, forward fft + // fftwf_plan inv_; // synthesis, inverse fft + // unsigned int w_; + // unsigned int h_; + // unsigned int h_complex_; + // unsigned int h_padded_; + // float * buffer_; + // }; + + //---------------------------------------------------------------- + // tss + // + // static boost::thread_specific_ptr tss; + using ForwardFFTFilterType = itk::ForwardFFTImageFilter; + static thread_local ForwardFFTFilterType::Pointer forwardFFTFilter = ForwardFFTFilterType::New(); + using InverseFFTFilterType = itk::InverseFFTImageFilter; + static thread_local InverseFFTFilterType::Pointer inverseFFTFilter = InverseFFTFilterType::New(); + + //---------------------------------------------------------------- + // fft_data_t::fft_data_t + // + fft_data_t::fft_data_t(const itk_image_t::Pointer & real): + image_(), + nx_(0), + ny_(0) + { + setup(real); + } + + //---------------------------------------------------------------- + // fft_data_t::fft_data_t + // + fft_data_t::fft_data_t(const unsigned int w, const unsigned int h): + image_(), + nx_(0), + ny_(0) + { + resize(w, h); + } + + //---------------------------------------------------------------- + // fft_data_t::fft_data_t + // + fft_data_t::fft_data_t(const itk_image_t::Pointer & real, + const itk_image_t::Pointer & imag): + image_(), + nx_(0), + ny_(0) + { + setup(real, imag); + } + + //---------------------------------------------------------------- + // fft_data_t::fft_data_t + // + fft_data_t::fft_data_t(const fft_data_t & data): + image_(), + nx_(0), + ny_(0) + { + *this = data; + } + + //---------------------------------------------------------------- + // fft_data_t::operator = + // + fft_data_t & + fft_data_t::operator = (const fft_data_t & data) + { + assert(this != &data); + + if (data.image_ != nullptr) + { + using DuplicatorType = itk::ImageDuplicator; + DuplicatorType::Pointer duplicator = DuplicatorType::New(); + duplicator->SetInputImage(data.image_); + duplicator->Update(); + image_ = duplicator->GetOutput(); + resize(data.nx_, data.ny_); + } + else + { + cleanup(); + } + + return *this; + } + + //---------------------------------------------------------------- + // fft_data_t::cleanup + // + void + fft_data_t::cleanup() + { + // if (data_ != NULL) fftwf_free((fft_complex_t *)(data_)); + // data_ = NULL; + nx_ = 0; + ny_ = 0; + } + + //---------------------------------------------------------------- + // fft_data_t::resize + // + void + fft_data_t::resize(const unsigned int w, const unsigned int h) + { + const unsigned int new_sz = w * h; + const unsigned int old_sz = nx_ * ny_; + if (old_sz == new_sz) return; + + // if (data_ != NULL) fftwf_free((fft_complex_t *)(data_)); + // data_ = (fft_complex_t *)(fftwf_malloc(new_sz * sizeof(fft_complex_t))); + // assert(data_ != NULL); + image_->SetRegions({ w, h }); + image_->Allocate(); + + nx_ = w; + ny_ = h; + } + + //---------------------------------------------------------------- + // fft_data_t::fill + // + void + fft_data_t::fill(const float real, const float imag) + { + const fft_complex_t value(real, imag); + this->image_->FillBuffer(value); + } + + //---------------------------------------------------------------- + // fft_data_t::setup + // + void + fft_data_t::setup(const itk_image_t::Pointer & real, + const itk_image_t::Pointer & imag) + { + typedef itk::ImageRegionConstIteratorWithIndex itex_t; + typedef itk_image_t::IndexType index_t; + + itk::Size<2> size = real->GetLargestPossibleRegion().GetSize(); + resize(size[0], size[1]); + + // shortcut: + fftwf_complex * hack = reinterpret_cast(data_); + + // iterate over the real component image: + itex_t itex(real, real->GetLargestPossibleRegion()); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + const index_t index = itex.GetIndex(); + const unsigned int x = index[0]; + const unsigned int y = index[1]; + const unsigned int i = y + ny_ * x; + + hack[i][0] = itex.Get(); + hack[i][1] = 0.0; + } + + // iterate over the imaginary component image: + if (imag.GetPointer() == NULL) return; + itex = itex_t(imag, real->GetLargestPossibleRegion()); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + const index_t index = itex.GetIndex(); + const unsigned int x = index[0]; + const unsigned int y = index[1]; + const unsigned int i = y + ny_ * x; + + hack[i][1] = itex.Get(); + } + } + + //---------------------------------------------------------------- + // fft_data_t::component + // + itk_image_t::Pointer + fft_data_t::component(const bool imag) const + { + typedef itk::ImageRegionIteratorWithIndex itex_t; + typedef itk_image_t::IndexType index_t; + + itk_image_t::Pointer out = itk_image_t::New(); + itk::Size<2> size; + size[0] = nx_; + size[1] = ny_; + out->SetRegions(size); + out->Allocate(); + + // shortcut: + const fftwf_complex * hack = reinterpret_cast(data_); + + // iterate over the image: + itex_t itex(out, out->GetLargestPossibleRegion()); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + const index_t index = itex.GetIndex(); + const unsigned int x = index[0]; + const unsigned int y = index[1]; + const unsigned int i = y + ny_ * x; + + itex.Set(hack[i][imag]); + } + + return out; + } + + static std::vector syLookup; + static std::vector sySquaredLookup; + + //---------------------------------------------------------------- + // fft_data_t::apply_lp_filter + // + void + fft_data_t::apply_lp_filter(const double r, const double s) + { + if (r > ::sqrt(2.0)) return; + + const unsigned int hx = nx_ / 2; + const unsigned int hy = ny_ / 2; + + const double r0 = (r - s); + const double r1 = (r + s); + const double dr = r1 - r0; + + const double r0sqr = r0 * r0; + const double r1sqr = r1 * r1; + + //The use of this static cache needs to be wrapped in a lock, copied, or some other protection if multithreaded calls with different image sizes are used. + if((unsigned int)syLookup.size() != ny_) + { + syLookup.clear(); + sySquaredLookup.clear(); + + syLookup.reserve(ny_); + sySquaredLookup.reserve(ny_); + for (unsigned int y = 0; y < ny_; y++) + { + double sy = 2.0 * (double((y + hy) % ny_) - double(hy)) / double(ny_); + double y2 = sy * sy; + syLookup.push_back(sy); + sySquaredLookup.push_back(sy); + } + } + + float * data_ = this->image_->GetBufferPointer(); + unsigned int i = 0; + for (unsigned int x = 0; x < nx_; x++) + { + double sx = 2.0 * (double((x + hx) % nx_) - double(hx)) / double(nx_); + double x2 = sx * sx; + + for (unsigned int y = 0; y < ny_; y++) + { + double sy = syLookup[y]; + double y2 = sySquaredLookup[y]; + + double d2 = x2 + y2; + double v; + + if(d2 < r0sqr) + v = 1.0; + else if(d2 > r1sqr) + v = 0.0; + else + { + double d = ::sqrt(x2 + y2); + v = (1.0 + cos(M_PI * (d - r0) / dr)) / 2.0; + } + + data_[i++] *= v; + } + } + } + + + //---------------------------------------------------------------- + // fft + // + bool + fft(itk_image_t::ConstPointer & in, fft_data_t & out) + { + typedef itk::ImageRegionConstIteratorWithIndex itex_t; + typedef itk_image_t::IndexType index_t; + + itk::Size<2> size = in->GetLargestPossibleRegion().GetSize(); + const unsigned int w = size[0]; + const unsigned int h = size[1]; + + fft_cache_t * cache = tss.get(); + if (!cache) + { + cache = new fft_cache_t(); + tss.reset(cache); + } + cache->update(w, h); + + // iterate over the image: + itex_t itex(in, in->GetLargestPossibleRegion()); + for (itex.GoToBegin(); !itex.IsAtEnd(); ++itex) + { + const index_t index = itex.GetIndex(); + const unsigned int x = index[0]; + const unsigned int y = index[1]; + const unsigned int i = y + cache->h_padded_ * x; + + cache->buffer_[i] = itex.Get(); + } + + if (!cache->fwd_) return false; + fftwf_execute_dft_r2c(cache->fwd_, + cache->buffer_, + (fftwf_complex *)(cache->buffer_)); + + // shortcuts: + const unsigned int h_complex = cache->h_complex_; + fft_complex_t * buffer = (fft_complex_t *)cache->buffer_; + + // fill in the rest of the output data: + out.resize(w, h); + for (unsigned int x = 0; x < w; x++) + { + for (unsigned int y = 0; y < h_complex; y++) + { + unsigned int i = y + h_complex * x; + out(x, y) = buffer[i]; + } + + for (unsigned int y = h_complex; y < h; y++) + { + unsigned int i = (h - y) + h_complex * ((w - x) % w); + out(x, y) = std::conj(buffer[i]); + } + } + + return true; + } + + + //---------------------------------------------------------------- + // ifft + // + bool + ifft(const fft_data_t & in, fft_data_t & out) + { + out.resize(in.nx(), in.ny()); + + fft_cache_t * cache = tss.get(); + if (!cache) + { + return false; + } + + fftwf_execute_dft(cache->inv_, + (fftwf_complex *)(in.data()), + (fftwf_complex *)(out.data())); + return true; + } + + //---------------------------------------------------------------- + // elem_by_elem + // + void + elem_by_elem(fn_fft_c_t f, + const fft_data_t & in, + fft_data_t & out) + { + const unsigned nx = in.nx(); + const unsigned ny = in.ny(); + + out.resize(nx, ny); + fft_complex_t * dout = out.data(); + const fft_complex_t * din = in.data(); + + const unsigned int size = nx * ny; + for (unsigned int i = 0; i < size; i++) + { + dout[i] = f(din[i]); + } + } + + + //---------------------------------------------------------------- + // elem_by_elem + // + void + elem_by_elem(fft_complex_t(*f)(const double & a, + const fft_complex_t & b), + const double & a, + const fft_data_t & b, + fft_data_t & c) + { + assert(&b != &c); + + const unsigned nx = b.nx(); + const unsigned ny = b.ny(); + c.resize(nx, ny); + + fft_complex_t * dc = c.data(); + const fft_complex_t * db = b.data(); + + const unsigned int size = nx * ny; + for (unsigned int i = 0; i < size; i++) + { + dc[i] = f(a, db[i]); + } + } + + + //---------------------------------------------------------------- + // elem_by_elem + // + void + elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + const double & b), + const fft_data_t & a, + const double & b, + fft_data_t & c) + { + assert(&a != &c); + + const unsigned nx = a.nx(); + const unsigned ny = a.ny(); + c.resize(nx, ny); + + fft_complex_t * dc = c.data(); + const fft_complex_t * da = a.data(); + + const unsigned int size = nx * ny; + for (unsigned int i = 0; i < size; i++) + { + dc[i] = f(da[i], b); + } + } + + + //---------------------------------------------------------------- + // elem_by_elem + // + void + elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + const fft_complex_t & b), + const fft_complex_t & a, + const fft_data_t & b, + fft_data_t & c) + { + assert(&b != &c); + + const unsigned nx = b.nx(); + const unsigned ny = b.ny(); + c.resize(nx, ny); + + fft_complex_t * dc = c.data(); + const fft_complex_t * db = b.data(); + + const unsigned int size = nx * ny; + for (unsigned int i = 0; i < size; i++) + { + dc[i] = f(a, db[i]); + } + } + + + //---------------------------------------------------------------- + // elem_by_elem + // + void + elem_by_elem(fft_complex_t(*f)(const fft_complex_t & a, + const fft_complex_t & b), + const fft_data_t & a, + const fft_complex_t & b, + fft_data_t & c) + { + assert(&a != &c); + + const unsigned nx = a.nx(); + const unsigned ny = a.ny(); + c.resize(nx, ny); + + fft_complex_t * dc = c.data(); + const fft_complex_t * da = a.data(); + + const unsigned int size = nx * ny; + for (unsigned int i = 0; i < size; i++) + { + dc[i] = f(da[i], b); + } + } + + + //---------------------------------------------------------------- + // elem_by_elem + // + void + elem_by_elem(fn_fft_cc_t f, + const fft_data_t & a, + const fft_data_t & b, + fft_data_t & c) + { + assert(&a != &c); + assert(&b != &c); + + const unsigned nx = a.nx(); + const unsigned ny = a.ny(); + assert(nx == b.nx() && ny == b.ny()); + + c.resize(nx, ny); + fft_complex_t * dc = c.data(); + + const fft_complex_t * da = a.data(); + const fft_complex_t * db = b.data(); + + const unsigned int size = nx * ny; + for (unsigned int i = 0; i < size; i++) + { + dc[i] = f(da[i], db[i]); + } + } + +} // namespace itk_fft