FFT2D.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
Multimediale Systeme der Informationsverarbeitung
Institut fuer Informatik
Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


#include "FFT2D.hh"
#include <fftw3.h>
#include <Base/Debug/Error.hh>

#include <complex>
#include <iomanip>

using namespace BIAS;
using namespace std;

template <class InputStorageType, class OutputStorageType>
FFT2D<InputStorageType, OutputStorageType>::
FFT2D()
{
  _SizeX=0;
  _SizeY=0;
  _OutNum=0;
  _in=NULL;
  _out=NULL;
  _p_forward=NULL;
}

template <class InputStorageType, class OutputStorageType>
FFT2D<InputStorageType, OutputStorageType>::~FFT2D()
{
  Release();
  fftw_cleanup();
}

template <class InputStorageType, class OutputStorageType>
void FFT2D<InputStorageType, OutputStorageType>::
Init(int width, int height)
{
#ifdef BIAS_DEBUG
  if (_SizeX!=0 || _SizeY!=0 ||  _OutNum!=0 || _out || _in || _p_forward){
    BIASERR("FFT2D: call Release() before calling Init() a second time");
    BIASABORT;
  }
#endif
  _SizeX=height;
  _SizeY=width;
  _Size = _SizeX*_SizeY;
  _OutSizeY = (_SizeY>>1)+1;
  _OutNum=_SizeX * _OutSizeY;
  _out =(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _OutNum);
  _in = (double*)fftw_malloc(sizeof(double) * _Size);
  // this is slow a few seconds
  _p_forward = fftw_plan_dft_r2c_2d(_SizeX, _SizeY, _in, _out, FFTW_MEASURE);
  _p_reverse = fftw_plan_dft_c2r_2d(_SizeX, _SizeY, _out, _in, FFTW_ESTIMATE);
}






template <class InputStorageType, class OutputStorageType>
void FFT2D<InputStorageType, OutputStorageType>::
Release()
{
  if (_p_forward) fftw_destroy_plan(_p_forward);
  if (_in) fftw_free(_in); 
  if (_out) fftw_free(_out);
  _in=NULL;
  _out=NULL;
  _p_forward=NULL;
  _SizeX=_SizeY=_Size=_OutNum=0;
}


template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
Forward_(const Image<InputStorageType>& src)
{
#ifdef BIAS_DEBUG
  if (_SizeX==0 || _SizeY==0 ){
    BIASERR("FFT2D: call Init() before calling Filter() ");
    BIASABORT;
  }
#endif

  BIASASSERT(src.GetChannelCount() ==1);
  BIASASSERT((int)src.GetWidth() == _SizeY);
  BIASASSERT((int)src.GetHeight() == _SizeX);
  // first, copy and cast src into _in
  double *inp = _in;
  const InputStorageType *imgp = src.GetImageData();
  for ( int i=0;i<_Size; i++) {
    inp[i] = (double) imgp[i];
  }

  // run fft via _plan
  fftw_execute(_p_forward);
  return 0;
}


template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
Reverse_(const Image<OutputStorageType>& src)
{
#ifdef BIAS_DEBUG
  if (_SizeX==0 || _SizeY==0 ){
    BIASERR("FFT2D: call Init() before calling Filter() ");
    BIASABORT;
  }
#endif

  BIASASSERT(src.GetChannelCount() ==2);
  BIASASSERT((int)src.GetWidth() == _OutSizeY);
  BIASASSERT((int)src.GetHeight() == _SizeX);
  // first, copy and cast src into _in
  fftw_complex *outp = _out;
  const OutputStorageType *imgp = src.GetImageData();
  for ( int i=0;i<_OutNum; i++) {
    outp[i][0] = (double) imgp[2*i];
    outp[i][1] = (double) imgp[2*i+1];
  }

  // run fft via _plan
  fftw_execute(_p_reverse);
  return 0;
}




/** dst.GetChannelCount()==2*src.GetCHannelCount() */
template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
Filter(const Image<InputStorageType>& src, Image<OutputStorageType>& dst)
{
  Forward_(src);

  if (!dst.IsEmpty()) dst.Release();
  dst.Init(_OutSizeY, _SizeX, 2);
  // copy and cast back int dst, two channel mode
  fftw_complex *outp = _out;
  OutputStorageType *oimgp = dst.GetImageData();
  for ( int i=0;i<_OutNum; i++){
    oimgp[2*i] = (OutputStorageType) outp[i][0];
    oimgp[2*i+1] = (OutputStorageType) outp[i][1];
  }
  return 0;
}



/** dstX.GetChannelCount()==src.GetCHannelCount() */
template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
Filter(const Image<InputStorageType>& src, Image<OutputStorageType>& dst1, 
       Image<OutputStorageType>& dst2)
{
  Forward_(src);

  if (!dst1.IsEmpty()) dst1.Release();
  dst1.Init(_OutSizeY, _SizeX,1);
  if (!dst2.IsEmpty()) dst2.Release();
  dst2.Init(_OutSizeY, _SizeX,1);
  // copy and cast back int dst, two channel mode
  fftw_complex *outp = _out;
  OutputStorageType *oimgp1 = dst1.GetImageData();
  OutputStorageType *oimgp2 = dst2.GetImageData();
  for ( int i=0;i<_OutNum; i++){
    oimgp1[i] = (OutputStorageType) (outp[i][0]);
    oimgp2[i] = (OutputStorageType) (outp[i][1]);
  }

  return 0;
}


/** dstX.GetChannelCount()==src.GetCHannelCount() */
template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
TransformAbs(const Image<InputStorageType>& src, Image<OutputStorageType>& dst)
{
  Forward_(src);
  if (!dst.IsEmpty()) dst.Release();
  dst.Init(_OutSizeY ,_SizeX,1);
  // copy and cast back int dst, abs mode
  fftw_complex *outp = _out;
  OutputStorageType *oimgp = dst.GetImageData();
  for ( int i=0;i<_OutNum; i++){

    oimgp[i] = 
      (OutputStorageType) (sqrt((outp[i][0] * outp[i][0]) + (outp[i][1] * outp[i][1]))/_Size);
  }

  return 0;
}

/** dstX.GetChannelCount()==src.GetCHannelCount() */
template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
TransformLogAbs(const Image<InputStorageType>& src, Image<OutputStorageType>& dst)
{
  Forward_(src);
  if (!dst.IsEmpty()) dst.Release();
  dst.Init( _OutSizeY ,_SizeX,1);
  // copy and cast back int dst, abs mode
  fftw_complex *outp = _out;
  OutputStorageType *oimgp = dst.GetImageData();
  for ( int i=0;i<_OutNum; i++){
    oimgp[i] = (OutputStorageType) log(sqrt((outp[i][0]* outp[i][0]) + (outp[i][1]* outp[i][1])));
  }

  return 0;
}

/** dstX.GetChannelCount()==src.GetCHannelCount() */
template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
TransformPhase(const Image<InputStorageType>& src, Image<OutputStorageType>& dst)
{
  Forward_(src);
  if (!dst.IsEmpty()) dst.Release();
  dst.Init( _OutSizeY ,_SizeX,1);
  // copy and cast back int dst, abs mode
  fftw_complex *outp = _out;
  OutputStorageType *oimgp = dst.GetImageData();
  for ( int i=0;i<_OutNum; i++){
    oimgp[i] = (OutputStorageType) atan2f((float)outp[i][0], (float)outp[i][1]);
  }

  return 0;
}

#include <Base/Common/BIASpragmaStart.hh>
template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
Normalize(Image<OutputStorageType>& dst)
{
  OutputStorageType *d = dst.GetImageData();
  unsigned int count = dst.GetWidth() * dst.GetHeight() *dst.GetChannelCount();
  for (unsigned int i=0; i<count; i++)
  {
    // TODO
    ///\todo fix the cast by explicit template JW
    d[i] /= (OutputStorageType) _Size;
  }
  return 0;
}

template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
CrossPowerSpectrum(const Image<InputStorageType>& src1, const Image<InputStorageType>& src2, Image<OutputStorageType>& dst1){

  fftw_complex *_out2;
  _out2 =(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * _OutNum);
  Forward_(src2);
  //copy second FFT to _out2;
  fftw_complex *outp1 = _out;
  fftw_complex *outp2 = _out2;
  for ( int i=0;i<_OutNum; i++){
    outp2[i][0] =  outp1[i][0];
    outp2[i][1] =  outp1[i][1];
  }
  Forward_(src1);

  if (!dst1.IsEmpty()) dst1.Release();
  //dst1.Init(_OutSizeY, _SizeX,2);
  dst1.Init(_OutSizeY, _SizeX,2);
  // copy and cast back int dst, two channel mode
  double abs;
  outp2 = _out2;
  outp1 = _out;
  OutputStorageType *oimgp1 = dst1.GetImageData();

  for ( int i=0;i<_OutNum; i++){
    abs = (double)(outp2[i][0] * outp2[i][0] + outp2[i][1] * outp2[i][1]);
    oimgp1[2*i] = (OutputStorageType) ((outp1[i][0] * outp2[i][0] + outp1[i][1] * outp2[i][1])/abs);
    oimgp1[2*i+1] = (OutputStorageType) ((outp2[i][0] * outp1[i][1] - outp1[i][0] * outp2[i][1])/abs);  
  }

  return 0;
}

#include <Base/Common/BIASpragmaEnd.hh>



template <class InputStorageType, class OutputStorageType>
int FFT2D<InputStorageType, OutputStorageType>::
TransformReverse(const Image<OutputStorageType>& src,
                 Image<InputStorageType>& dst)
{
  Reverse_(src);
  if (!dst.IsEmpty()) dst.Release();
  dst.Init(_SizeY, _SizeX,1);
  // copy and cast back int dst, abs mode
  double *inp = _in;
  InputStorageType *oimgp = dst.GetImageData();
  for ( int i=0;i<_Size; i++){
    oimgp[i] = (InputStorageType) inp[i];
  }
  return 0;
}

template <class InputStorageType, class OutputStorageType>
void FFT2D<InputStorageType, OutputStorageType>::
GetBordersValid_(int& /*border_x*/, int& /*border_y*/) const
{
  BIASERR("no parameter support, yet");
  BIASABORT;
}



//////////////////////////////////////////////////////////////////////////
//    explicit  instantiation    (must be at end of file)
//////////////////////////////////////////////////////////////////////////

#if defined(BUILD_IMAGE_FLOAT)
template class BIASImage_EXPORT FFT2D<float, float>;
#endif
#if defined(BUILD_IMAGE_DOUBLE)
template class BIASImage_EXPORT FFT2D<double, double>;
#endif
#if defined(BUILD_IMAGE_UCHAR) && defined(BUILD_IMAGE_FLOAT)
template class BIASImage_EXPORT FFT2D<unsigned char, float>;
#endif
#if defined(BUILD_IMAGE_UCHAR) && defined(BUILD_IMAGE_DOUBLE)
template class BIASImage_EXPORT FFT2D<unsigned char, double>;
#endif