FFT2D_Tiles.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_Tiles.hh"
#include <fftw3.h>
#include <Base/Debug/Error.hh>


using namespace BIAS;
using namespace std;

template <class InputStorageType, class OutputStorageType>
FFT2D_Tiles<InputStorageType, OutputStorageType>::
FFT2D_Tiles()
{
  this->_TileSizeX = 0;
  this->_TileSizeY = 0;
  this->_NumTilesX = 0;
  this->_NumTilesY = 0;
}

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

template <class InputStorageType, class OutputStorageType>
void FFT2D_Tiles<InputStorageType, OutputStorageType>::
Init(int width, int height, int tilewidth, int tileheight)
{
#ifdef BIAS_DEBUG
  if (this->_SizeX!=0 || this->_SizeY!=0 ||  this->_OutNum!=0 || this->_out || this->_in || this->_p_forward){
    BIASERR("FFT2D_Tiles: call Release() before calling Init() a second time");
    BIASEXIT(EXIT_FAILURE, "FFT2D_Tiles: call Release() before calling Init() a second time");
  }
#endif
  this->_SizeX=height;
  this->_SizeY=width;

  this->_TileSizeX = tileheight;
  this->_TileSizeY = tilewidth;
  this->_NumTilesX = this->_SizeX / this->_TileSizeX;
  if (this->_SizeX % this->_TileSizeX >0) this->_NumTilesX++;
  this->_NumTilesY = this->_SizeY / this->_TileSizeY;
  if (this->_SizeY % this->_TileSizeY >0) this->_NumTilesY++;



  this->_Size = this->_TileSizeX*this->_TileSizeY;
  this->_OutSizeY = (this->_TileSizeY>>1)+1;
  this->_OutNum=this->_TileSizeX * this->_OutSizeY;
  this->_out =(fftw_complex*)fftw_malloc(sizeof(fftw_complex) * this->_OutNum);
  this->_in = (double*)fftw_malloc(sizeof(double) * this->_Size);
  // this is slow a few seconds
  this->_p_forward = fftw_plan_dft_r2c_2d(this->_TileSizeX, this->_TileSizeY, this->_in, this->_out, FFTW_MEASURE);
  this->_p_reverse = fftw_plan_dft_c2r_2d(this->_TileSizeX, this->_TileSizeY, this->_out, this->_in, FFTW_MEASURE);
}






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


template <class InputStorageType, class OutputStorageType>
int FFT2D_Tiles<InputStorageType, OutputStorageType>::
Forward_(const Image<InputStorageType>& /*src*/)
{
  BIASERR("not implemented");
  return 0;
}


template <class InputStorageType, class OutputStorageType>
int FFT2D_Tiles<InputStorageType, OutputStorageType>::
Reverse_(const Image<OutputStorageType>& /*src*/)
{
  BIASERR("not implemented");
  return 0;
}




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

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



/** dstX.GetChannelCount()==src.GetCHannelCount() */
template <class InputStorageType, class OutputStorageType>
int FFT2D_Tiles<InputStorageType, OutputStorageType>::
Filter(const Image<InputStorageType>& /*src*/, 
       Image<OutputStorageType>& /*dst1*/, 
       Image<OutputStorageType>& /*dst2*/)
{
  BIASERR("not implemented");
  return 0;
}


/** dst.GetChannelCount()==src.GetCHannelCount() */
template <class InputStorageType, class OutputStorageType>
int FFT2D_Tiles<InputStorageType, OutputStorageType>::
TransformAbs(const Image<InputStorageType>& src, Image<OutputStorageType>& dst)
{
  Image<InputStorageType> tmpImg;
  InputStorageType **Iida ;

  if ( ((int) src.GetHeight()  != this->_TileSizeX * this->_NumTilesX)
    || ((int) src.GetWidth() != this->_TileSizeY * this->_NumTilesY)) {
      tmpImg = src;
      tmpImg.ZeroPad( this->_TileSizeX * this->_NumTilesX, this->_TileSizeY * this->_NumTilesY );
      Iida = tmpImg.GetImageDataArray();
    }
  else {
    tmpImg = src;
    Iida = tmpImg.GetImageDataArray();
  }
  // prepare Image
  if (!dst.IsEmpty()) dst.Release();
  dst.Init( this->_OutSizeY * this->_NumTilesY,this->_TileSizeX * this->_NumTilesX,1);

  OutputStorageType **Oida = dst.GetImageDataArray();
  // for all tiles
  for ( int ty=0; ty<this->_NumTilesY; ty++) {
    for ( int tx=0; tx<this->_NumTilesX; tx++) {
      double *inp = this->_in;
      fftw_complex *outp = this->_out;

      // copy this tile int _in;
      for ( int y = ty*this->_TileSizeY; y<(ty+1)*this->_TileSizeY; y++)
        for ( int x = tx*this->_TileSizeX; x<(tx+1)*this->_TileSizeX; x++) {
          *inp++ = Iida[y][x];
        }
        // fft it
        fftw_execute(this->_p_forward);
        // copy result back
        for ( int y = ty*this->_OutSizeY; y<(ty+1)*this->_OutSizeY; y++)
          for ( int x = tx*this->_TileSizeX; x<(tx+1)*this->_TileSizeX; x++) {
            // remeber: x is height and y is width
            Oida[x][y] = (OutputStorageType) sqrt((*outp[0]* *outp[0]) + 
              (*outp[1]* *outp[1]));
            outp++;
          }
    }
  }

  return 0;
}





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

  return 0;
}

//////////////////////////////////////////////////////////////////////////
//                 instantiation
//////////////////////////////////////////////////////////////////////////
//#define FILTER_INSTANTIATION_CLASS FFT2D_Tiles
//#define FILTER_INSTANTIATION_NO_UNSIGNED_OUTPUT
//#include "Filter/Filterinst.hh"


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