[Insight-developers] FFTW : undefined reference
Paul Koshevoy
koshevoy at sci.utah.edu
Thu Jun 29 12:11:34 EDT 2006
Hi Alan,
I've attached my fft wrapper classes. I don't want to claim that they do
everything -- I use them for fast translation estimation between 2
images. Below is a function to give you some context and help you see
how the wrappers are used.
NOTE: itk_image_t and fftw_data_t are declared in fft.hxx. The function
pads the two images on the right and the bottom to make sure both images
are of the same size (so that they can be multiplied together in the
frequency domain).
//----------------------------------------------------------------
// find_correlation
//
template <class image_t>
unsigned int
find_correlation(const image_t * fi,
const image_t * mi,
std::list<local_max_t> & max_list,
bool resampled_data,
const the_text_t & prefix = the_text_t(""),
const the_text_t & suffix = the_text_t(".tif"))
{
itk_image_t::Pointer z0 = cast<image_t, itk_image_t>(fi);
itk_image_t::Pointer z1 = cast<image_t, itk_image_t>(mi);
itk_image_t::SizeType max_sz = calc_padding<itk_image_t>(z0, z1);
z0 = pad<itk_image_t>(z0, max_sz);
z1 = pad<itk_image_t>(z1, max_sz);
fftw_data_t f0;
fft(z0, f0);
// resampled data requires less smoothing:
if (resampled_data) f0.apply_lp_filter(0.9, 0.1);
else f0.apply_lp_filter(0.5, 0.1);
fftw_data_t f1;
fft(z1, f1);
// resampled data requires less smoothing:
if (resampled_data) f1.apply_lp_filter(0.9, 0.1);
else f1.apply_lp_filter(0.5, 0.1);
const unsigned int & nx = f0.nx();
const unsigned int & ny = f0.ny();
fftw_data_t P(nx, ny);
for (unsigned int x = 0; x < nx; x++)
{
for (unsigned int y = 0; y < ny; y++)
{
#if 1
// Girod-Kuo:
P(x, y) =
(f1(x, y) * std::conj(f0(x, y))) /
(std::sqrt(f0(x, y) * std::conj(f0(x, y)) *
f1(x, y) * std::conj(f1(x, y))) + 1e-8);
#else
// plain correlation:
P(x, y) = f1(x, y) * std::conj(f0(x, y));
#endif
}
}
// resampled data produces less noisy PDF and requires less smoothing:
if (resampled_data) P.apply_lp_filter(0.6, 0.1);
else P.apply_lp_filter(0.4, 0.1);
// calculate the displacement probability density function:
fftw_data_t ifft_P;
bool ok = ifft(P, ifft_P);
assert(ok);
itk_image_t::Pointer PDF = ifft_P.real();
// look for the maxima in the PDF:
double area = double(max_sz[0] * max_sz[1]);
// a minimum of 5 pixels and a maximum of 64 pixels may be attributed
// to local maxima in the image:
// double fraction = std::min(64.0 / area, std::max(5.0 / area, 1e-3));
double fraction = std::min(64.0 / area, std::max(5.0 / area, 1e-2));
// the entire image should never be treated as a maxima cluster:
assert(fraction < 1.0);
// find the maxima clusters:
find_maxima_cm(max_list, PDF, 1.0 - fraction, prefix, suffix);
return max_list.size();
}
Alan wrote:
> Hi folks,
>
> Finally, Julien was right :
> http://public.kitware.com/pipermail/insight-users/2005-November/015757.html
> Nothing more, nothing less.
> Paul, I would be interested in testing your fftw wrapper class.
> Thx in advance.
>
> Alan
>
-------------- next part --------------
#ifndef FFT_HXX_
#define FFT_HXX_
// FFTW includes:
#include <fftw3.h>
// ITK includes:
#include <itkImage.h>
// system includes:
#include <complex>
#include <assert.h>
namespace itk_fftw
{
//----------------------------------------------------------------
// itk_image_t
//
typedef itk::Image<double, 2> itk_image_t;
//----------------------------------------------------------------
// fftw_complex_t
//
typedef std::complex<double> fftw_complex_t;
//----------------------------------------------------------------
// fftw_data_t
//
class fftw_data_t
{
public:
fftw_data_t(): data_(NULL), nx_(0), ny_(0) {}
fftw_data_t(const unsigned int w, const unsigned int h);
explicit fftw_data_t(const itk_image_t::Pointer & real);
fftw_data_t(const itk_image_t::Pointer & real,
const itk_image_t::Pointer & imag);
fftw_data_t(const fftw_data_t & data);
~fftw_data_t() { cleanup(); }
fftw_data_t & operator = (const fftw_data_t & data);
void cleanup();
void resize(const unsigned int w, const unsigned int h);
void fill(const double real, const double imag = 0.0);
void setup(const itk_image_t::Pointer & real,
const itk_image_t::Pointer & imag = itk_image_t::Pointer(NULL));
// ITK helpers:
itk_image_t::Pointer component(const bool imag = 0) const;
inline itk_image_t::Pointer real() const { return component(0); }
inline itk_image_t::Pointer imag() const { return component(1); }
// Apply a low-pass filter to this image. This function will
// zero-out high-frequency components, where the cutoff frequency
// is specified by radius r in [0, 1]. The sharpness of the
// cutoff may be controlled by parameter s, where s == 0 results in
// an ideal low-pass filter, and s == 1 is a low pass filter defined
// by a scaled and shifted cosine function, 1 at the origin,
// 0.5 at the cutoff frequency and 0 at twice the cutoff frequency.
void apply_lp_filter(const double r, const double s = 0.0);
// accessors:
inline const fftw_complex_t * data() const { return data_; }
inline fftw_complex_t * data() { return data_; }
inline const unsigned int nx() const { return nx_; }
inline const unsigned int ny() const { return ny_; }
inline const fftw_complex_t & operator() (const unsigned int & x,
const unsigned int & y) const
{ return at(x, y); }
inline fftw_complex_t & operator() (const unsigned int & x,
const unsigned int & y)
{ return at(x, y); }
inline const fftw_complex_t & at(const unsigned int & x,
const unsigned int & y) const
{ return data_[y + ny_ * x]; }
inline fftw_complex_t & at(const unsigned int & x,
const unsigned int & y)
{ return data_[y + ny_ * x]; }
protected:
fftw_complex_t * data_;
unsigned int nx_;
unsigned int ny_;
};
//----------------------------------------------------------------
// fft
//
extern bool
fft(const itk_image_t::Pointer & in, fftw_data_t & out);
//----------------------------------------------------------------
// fft
//
inline fftw_data_t
fft(const itk_image_t::Pointer & in)
{
fftw_data_t out;
bool ok = fft(in, out);
assert(ok);
return out;
}
//----------------------------------------------------------------
// fft
//
extern bool
fft(const fftw_data_t & in, fftw_data_t & out, int sign = FFTW_FORWARD);
//----------------------------------------------------------------
// fft
//
inline fftw_data_t
fft(const fftw_data_t & in)
{
fftw_data_t out;
bool ok = fft(in, out);
assert(ok);
return out;
}
//----------------------------------------------------------------
// ifft
//
extern bool
ifft(const fftw_data_t & in, fftw_data_t & out);
//----------------------------------------------------------------
// ifft
//
inline fftw_data_t
ifft(const fftw_data_t & in)
{
fftw_data_t out;
bool ok = ifft(in, out);
assert(ok);
return out;
}
//----------------------------------------------------------------
// fn_fftw_c_t
//
typedef fftw_complex_t(*fn_fftw_c_t)(const fftw_complex_t & in);
//----------------------------------------------------------------
// fn_fftw_cc_t
//
typedef fftw_complex_t(*fn_fftw_cc_t)(const fftw_complex_t & a,
const fftw_complex_t & b);
//----------------------------------------------------------------
// elem_by_elem
//
extern void
elem_by_elem(fn_fftw_c_t f,
const fftw_data_t & in,
fftw_data_t & out);
//----------------------------------------------------------------
// elem_by_elem
//
extern void
elem_by_elem(fftw_complex_t(*f)(const double & a,
const fftw_complex_t & b),
const double & a,
const fftw_data_t & b,
fftw_data_t & c);
//----------------------------------------------------------------
// elem_by_elem
//
extern void
elem_by_elem(fftw_complex_t(*f)(const fftw_complex_t & a,
const double & b),
const fftw_data_t & a,
const double & b,
fftw_data_t & c);
//----------------------------------------------------------------
// elem_by_elem
//
extern void
elem_by_elem(fftw_complex_t(*f)(const fftw_complex_t & a,
const fftw_complex_t & b),
const fftw_complex_t & a,
const fftw_data_t & b,
fftw_data_t & c);
//----------------------------------------------------------------
// elem_by_elem
//
extern void
elem_by_elem(fftw_complex_t(*f)(const fftw_complex_t & a,
const fftw_complex_t & b),
const fftw_data_t & a,
const fftw_complex_t & b,
fftw_data_t & c);
//----------------------------------------------------------------
// elem_by_elem
//
extern void
elem_by_elem(fn_fftw_cc_t f,
const fftw_data_t & a,
const fftw_data_t & b,
fftw_data_t & c);
//----------------------------------------------------------------
// conj
//
// element-by-element complex conjugate:
//
inline void
conj(const fftw_data_t & in, fftw_data_t & out)
{ elem_by_elem(std::conj, in, out); }
//----------------------------------------------------------------
// conj
//
inline fftw_data_t
conj(const fftw_data_t & in)
{
fftw_data_t out;
conj(in, out);
return out;
}
//----------------------------------------------------------------
// sqrt
//
// element-by-element square root:
//
inline void
sqrt(const fftw_data_t & in, fftw_data_t & out)
{ elem_by_elem(std::sqrt, in, out); }
//----------------------------------------------------------------
// sqrt
//
inline fftw_data_t
sqrt(const fftw_data_t & in)
{
fftw_data_t out;
sqrt(in, out);
return out;
}
//----------------------------------------------------------------
// _mul
//
template <class a_t, class b_t>
inline fftw_complex_t
_mul(const a_t & a, const b_t & b)
{ return a * b; }
//----------------------------------------------------------------
// mul
//
// element-by-element multiplication, c = a * b:
//
template <class a_t, class b_t>
inline void
mul(const a_t & a, const b_t & b, fftw_data_t & c)
{ elem_by_elem(_mul, a, b, c); }
//----------------------------------------------------------------
// mul
//
template <class a_t, class b_t>
inline fftw_data_t
mul(const a_t & a, const b_t & b)
{
fftw_data_t c;
mul<a_t, b_t>(a, b, c);
return c;
}
//----------------------------------------------------------------
// _div
//
template <class a_t, class b_t>
inline fftw_complex_t
_div(const a_t & a, const b_t & b)
{ return a / b; }
//----------------------------------------------------------------
// div
//
// element-by-element division, c = a / b:
//
template <class a_t, class b_t>
inline void
div(const a_t & a, const b_t & b, fftw_data_t & c)
{ elem_by_elem(_div, a, b, c); }
//----------------------------------------------------------------
// div
//
template <class a_t, class b_t>
inline fftw_data_t
div(const a_t & a, const b_t & b)
{
fftw_data_t c;
div<a_t, b_t>(a, b, c);
return c;
}
//----------------------------------------------------------------
// _add
//
template <class a_t, class b_t>
inline fftw_complex_t
_add(const a_t & a, const b_t & b)
{ return a + b; }
//----------------------------------------------------------------
// add
//
// element-by-element addition, c = a + b:
//
template <class a_t, class b_t>
inline void
add(const a_t & a, const b_t & b, fftw_data_t & c)
{ elem_by_elem(_add, a, b, c); }
//----------------------------------------------------------------
// add
//
template <class a_t, class b_t>
inline fftw_data_t
add(const a_t & a, const b_t & b)
{
fftw_data_t c;
add<a_t, b_t>(a, b, c);
return c;
}
//----------------------------------------------------------------
// _sub
//
template <class a_t, class b_t>
inline fftw_complex_t
_sub(const a_t & a, const b_t & b)
{ return a - b; }
//----------------------------------------------------------------
// sub
//
// element-by-element subtraction, c = a - b:
//
template <class a_t, class b_t>
inline void
sub(const a_t & a, const b_t & b, fftw_data_t & c)
{ elem_by_elem(_sub, a, b, c); }
//----------------------------------------------------------------
// sub
//
template <class a_t, class b_t>
inline fftw_data_t
sub(const a_t & a, const b_t & b)
{
fftw_data_t c;
sub<a_t, b_t>(a, b, c);
return c;
}
}
#endif // FFT_HXX_
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fft.cxx
Type: text/x-c++src
Size: 12019 bytes
Desc: not available
Url : http://www.itk.org/mailman/private/insight-developers/attachments/20060629/ce7783bd/fft.cxx
More information about the Insight-developers
mailing list