FFT2D_free.cpp

#include "FFT2D_free.hh"
//#include "fftpack5/fftpack5.h"
#include "fftpack/fftpack.h"
#include <iostream>

using namespace BIAS;
using namespace std;

template <class StorageType> 
FFT2D_free<StorageType>::
FFT2D_free()
{
  WorkSize0_ = WorkSize1_ = 0;
  WorkArray0_ = WorkArray1_ =  NULL;
  TmpArray0_ = TmpArray1_ = NULL;
  TmpArraySize0_ = TmpArraySize1_ = 0;  
  Normalize_ = 0;
}


template <class StorageType> 
FFT2D_free<StorageType>::
~FFT2D_free()
{
  
}


template <class StorageType> 
void FFT2D_free<StorageType>::
Init(int sizex, int sizey)
{
  // if size matches, do nothing, if it it does not match, clear and re-init

  if (sizex != WorkSize0_) {
    delete[] TmpArray0_;
    delete[] WorkArray0_;
    TmpArray0_ = NULL;
    WorkArray0_ = NULL;
  }
  if (sizey != WorkSize1_) {
    delete[] TmpArray1_;
    delete[] WorkArray1_;
    TmpArray1_ = NULL;
    WorkArray1_ = NULL;
  }


  if (WorkArray0_ == NULL) {
    WorkSize0_ = sizex;
    TmpArraySize0_ = 4 * WorkSize0_+15;
    TmpArray0_ = new float[TmpArraySize0_];
    WorkArray0_ = new float[2*WorkSize0_];
    cffti_(&WorkSize0_,TmpArray0_);
  }
  if (WorkArray1_ == NULL) {
    WorkSize1_ = sizey;
    TmpArraySize1_ = 4 * WorkSize1_+15;
    TmpArray1_ = new float[TmpArraySize1_];
    WorkArray1_ = new float[2*WorkSize1_];
    cffti_(&WorkSize1_,TmpArray1_);
  }
  Normalize_ = sizex * sizey;
}



template <class StorageType> 
void FFT2D_free<StorageType>::
ComplexColumnToWorkArray1_(const Image<float> &ComplexImg,
                   unsigned int col)
{
  const float **ida=ComplexImg.GetImageDataArray();
  unsigned int count = ComplexImg.GetHeight();
  for (unsigned int j=0;j<count;j++){
    WorkArray1_[j*2] = ida[j][col*2];      
    WorkArray1_[j*2+1] = ida[j][col*2+1];
  }
}

template <class StorageType> 
void FFT2D_free<StorageType>::
ComplexRowToWorkArray0_(const Image<float> &ComplexImg,
                unsigned int row)
{
  const float **ida=ComplexImg.GetImageDataArray();
  memcpy(WorkArray0_,ida[row],ComplexImg.GetWidth()*2*sizeof(float));
  
}



template <class StorageType> 
void FFT2D_free<StorageType>::
WorkArray0ToComplexRow_(Image<float> &ComplexImg,unsigned int row)
{
  float **out_ida=ComplexImg.GetImageDataArray();
  memcpy(out_ida[row],WorkArray0_,ComplexImg.GetWidth()*2*sizeof(float));

}


template <class StorageType> 
void FFT2D_free<StorageType>::
WorkArray1ToComplexColumn_(Image<float> &ComplexImg,unsigned int col)
{
  float **out_ida=ComplexImg.GetImageDataArray();
  unsigned int count = ComplexImg.GetHeight(); 
    for (unsigned int j=0;j<count;j++){
      out_ida[j][col*2] = WorkArray1_[j*2];
      out_ida[j][col*2+1] = WorkArray1_[j*2+1];
    }

}


template <class StorageType> 
void FFT2D_free<StorageType>::
Forward(const Image<StorageType> &in, Image<float> &ComplexOut)
{
  if (in.GetChannelCount() != 1) 
    BEXCEPTION("Can not FFT on more then one channel");
  // first, extract ech row ad transform it
  unsigned int NumRows=in.GetHeight();
  unsigned int NumCols=in.GetWidth();
  unsigned int NewNumCols = NumCols/2 +1;
  Image<float> TmpImg(NumCols,NumRows,2);
  ComplexOut.Init(NewNumCols,NumRows,2);
  const  StorageType **in_ida=in.GetImageDataArray();

  Init(NumCols, NumRows);

  // compute fft for each row and put into ComplexOut
  for (unsigned int j=0;j<NumRows; j++) {
    memset(WorkArray0_,0,2*WorkSize0_*sizeof(float)); // 
    for (unsigned int i=0;i<NumCols;i++) 
      WorkArray0_[i*2] = in_ida[j][i];
    cfftf_(&WorkSize0_, WorkArray0_,TmpArray0_);
    WorkArray0ToComplexRow_(TmpImg, j);
  }

  //   now for each col, same procedure again
  // restrict to non-redundant half (see www.fftw3.org) 
  for (unsigned int i=0;i<NewNumCols; i++) {
    ComplexColumnToWorkArray1_(TmpImg,i);
    cfftf_(&WorkSize1_,WorkArray1_,TmpArray1_);
    WorkArray1ToComplexColumn_(ComplexOut,i);
  }
}

template <class StorageType> 
void FFT2D_free<StorageType>::
Reverse(const Image<float> &ComplexIn, Image<StorageType> &res)
{
  unsigned int NumRows=ComplexIn.GetHeight();
  unsigned int NumCols=ComplexIn.GetWidth();
  unsigned int NewNumCols = (NumCols-1)*2;
  res.Init(NewNumCols,NumRows,2);
  Image<float> tmpImg(NewNumCols,NumRows,2);
  // first the columns
  Init( NewNumCols,NumRows);
  for (unsigned int i=0;i<NumCols; i++) {
    ComplexColumnToWorkArray1_(ComplexIn,i);
    cfftb_(&WorkSize1_,WorkArray1_,TmpArray1_);
    WorkArray1ToComplexColumn_(tmpImg,i);
    WorkArray1ToComplexColumn_(tmpImg,NewNumCols-i-1);
  }



  //// now the rows, then into res
  StorageType **out_ida=res.GetImageDataArray();
  for (unsigned int j=0;j<NumRows; j++) {
    ComplexRowToWorkArray0_(tmpImg,j);
    cfftb_(&WorkSize0_,WorkArray0_,TmpArray0_);
    for (unsigned int i=0;i<NumCols*2; i++)
      out_ida[j][i] = (StorageType)WorkArray0_[i];
  }
}

template <class StorageType> 
void FFT2D_free<StorageType>::
Normalize(Image<float> &Image)
{
  float *data=Image.GetImageData();
  unsigned int count = Image.GetWidth() *Image.GetHeight() *
    Image.GetChannelCount();
  for (unsigned int i=0;i<count;i++)
    data[i] = data[i]/float(Normalize_);

}


template <class StorageType> 
void FFT2D_free<StorageType>::
GetMagnitude(const Image<float> &ComplexIn, Image<float> &value)
{
  if (!value.IsEmpty()) value.Release();
  value.Init(ComplexIn.GetWidth(),ComplexIn.GetHeight());
  
  const float *in=ComplexIn.GetImageData();
  float *out = value.GetImageData();
  unsigned int count = ComplexIn.GetWidth() *ComplexIn.GetHeight()*2;
  unsigned int j=0;
  for (unsigned int i=0;i<count;i+=2){
    out[j++] = sqrt(in[i]*in[i] + in[i+1]*in[i+1]);
  }



}


namespace BIAS{
template class BIASImage_EXPORT FFT2D_free< unsigned char>;

#if defined(BUILD_IMAGE_FLOAT)
template class BIASImage_EXPORT FFT2D_free< float>;
#endif
#if defined(BUILD_IMAGE_SHORT)
template class BIASImage_EXPORT FFT2D_free<short>;
#endif

}