ExampleFMT2D.cpp

Example for using the fast fourier Transform (FFT) 2D

Author

MIP

/* 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*/


/**
 * @example ExampleFMT2D.cpp
   @relates FFT2D
   @brief Example for using the fast fourier Transform (FFT) 2D
   @ingroup g_examples
   @author  MIP
*/

#include <Image/FFT2D.hh>
#include <Base/Image/Image.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Image/ImageConvert.hh>
#include <Image/LogPolarMapping.hh>
#include <Image/HomographyMapping.hh>
#include <Base/Debug/TimeMeasure.hh>
#include <Geometry/RMatrix.hh>


using namespace std;
using namespace BIAS;

int cut(const Image<float>& src, Image<float>& dst){
  int width = src.GetWidth();
  int height = src.GetHeight();

  dst.Init(width-1,src.GetHeight(),1);
  const float *simgp = src.GetImageData();
  float *dimgp = dst.GetImageData();

  for(int i=0;i<height;i++){
    for(int j=1;j<width;j++){
      dimgp[(j-1)+i*(width-1)] = simgp[j+i*width];
    }
  }
  return 0;
}

int mirror(const Image<float>& src, Image<float>& dst){

  int width = src.GetWidth();
  int height = src.GetHeight();
  int qsize = width*height/2;
  if (!dst.IsEmpty()) dst.Release();
  //dst.Init(2*width,height,1);
  dst.Init(width,height,1);
  const float *simgp = src.GetImageData();
  float *oimgp = dst.GetImageData();
  for(int i=0;i<qsize;i++){

    oimgp[qsize+i] =  log(simgp[i]);
    oimgp[i] = log(simgp[qsize+i]);
  }
//Highpass Filter ??
/*
  for(int y =0; y<height;y++){
    for(int x=0;x<width;x++){
      oimgp[y*width+x] -= (10.0 * cos(0.5*M_PI*(double)x/(double)width) * cos(0.5*M_PI*(double)(y-height/2)/(double)(height/2)));
    }
  }
*/
/*
  oimgp[qsize] -= 8;
  oimgp[qsize+1] -= 6;
  oimgp[qsize+2] -= 4;
  oimgp[qsize+3] -= 2;
  oimgp[qsize-width] -= 6;
  oimgp[qsize-width+1] -= 4;
  oimgp[qsize-width+2] -= 2;
  oimgp[qsize+width] -= 6;
  oimgp[qsize+width+1] -= 4;
  oimgp[qsize+width+2] -= 2;
  oimgp[qsize-2*width] -= 4;
  oimgp[qsize-2*width+1] -= 2;
  oimgp[qsize-3*width] -= 2;
  oimgp[qsize+2*width] -= 4;
  oimgp[qsize+2*width+1] -= 2;
  oimgp[qsize+3*width] -= 2;
*/

  return 0;
}

int castfloat(const Image<unsigned char>& src, Image<float>& dst){
  int width = src.GetWidth();
  int height = src.GetHeight();
  int size = width*height;
  dst.Init( width,height,src.GetChannelCount());
  const unsigned char *imgp1 = src.GetImageData();
  float *oimgp1 = dst.GetImageData();
  for(int i=0;i<size;i++){
    oimgp1[i] = (float) imgp1[i];
  }
  return 0;
}

int castchar(const Image<float>& src, Image<unsigned char>& dst){
  int width = src.GetWidth();
  int height = src.GetHeight();
  int size = width*height;
  dst.Init( width,height,src.GetChannelCount());
  const float *imgp1 = src.GetImageData();
  unsigned char *oimgp1 = dst.GetImageData();
  for(int i=0;i<size;i++){
    oimgp1[i] = (unsigned char) imgp1[i];
  }
  return 0;
}

int main(int argc, char *argv[])
{

  TimeMeasure t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11;

  //load input images and initialize components
  t1.Start();
  if (argc != 3){
    cerr << argv[0] << "  <image> <image2>\n";
    return -1;
  }

  Image<unsigned char> im1, im2,im1grey,im2grey, resuc;
  Image<float> im1greyfloat, im2correctedfloat, im2greyfloat, mag1, mag1cut, mag1mirrored, mag2, mag2cut, mag2mirrored, cps, res;
  float max, max2;
  double angle,scale;
  unsigned int* coords = (unsigned int*) calloc(2,sizeof(unsigned int));
  unsigned int* coords2 = (unsigned int*) calloc(2,sizeof(unsigned int));
  LogPolarMapping<float, float> mapperlp;
  HomographyMapping<unsigned char, unsigned char> MapperHM;
  HMatrix H(MatrixIdentity);
  HMatrix HTrans(MatrixIdentity);
  HMatrix HTransInv(MatrixIdentity);
  HMatrix HRot(MatrixIdentity);
  HMatrix HScale(MatrixIdentity);

  if (ImageIO::Load(argv[1], im1)!=0){
    BIASERR("error loading image "<<argv[1]);
    return -2;
  } 
  if (ImageIO::Load(argv[2], im2)!=0){
    BIASERR("error loading image "<<argv[2]);
    return -2;
  }

  int width  = im1.GetWidth();
  int height = im1.GetHeight();

  if((width != (int)im2.GetWidth()) || 
     (height != (int)im2.GetHeight())){
    BIASERR("images must have same size!");
    return -2;
  }

  FFT2D<float, float> fft;
  mapperlp.SetImageCenter(width/2,height/2);
  Image<float> im1lp(width,height, 1);
  Image<float> im2lp(width,height, 1);
  Image<unsigned char>imgcorrected(width, height, 1);
  Image<unsigned char>imgcorrected2(width, height, 1);

  t10.Start();
  fft.Init(width,height);
  t10.Stop();
  
  ImageConvert::ToGrey(im1,im1grey);
  //ImageIO::Save("0_Grey1",im1grey);
  ImageIO::Save("0_Grey1",im1grey);
  ImageConvert::ToGrey(im2,im2grey);
  //ImageIO::Save("0_Grey2",im2grey);
  ImageIO::Save("0_Grey2",im2grey);
  im1lp.SetColorModel(im1grey.GetColorModel());
  im2lp.SetColorModel(im2grey.GetColorModel());
  t1.Stop();

  //cast image to float
  t2.Start();
  castfloat(im1grey,im1greyfloat);
  castfloat(im2grey,im2greyfloat);
  t2.Stop();

  //calculate magnitude for grey float input images
  t3.Start();
  fft.TransformAbs(im1greyfloat,mag1);
  //ImageIO::Save("1_mag1",mag1);
  ImageIO::Save("1_mag1",mag1);
  fft.TransformAbs(im2greyfloat,mag2);
  //ImageIO::Save("1_mag2",mag2);
  ImageIO::Save("1_mag2",mag2);
  t3.Stop();

  //cut image if too large and mirror magnitude to have CT at left border center position and width exactly half the original images width
  t4.Start();
  if(mag1.GetWidth()%2 != 0){
    cut(mag1,mag1cut);
    cut(mag2,mag2cut);
  }
  mirror(mag1cut,mag1mirrored);
  mirror(mag2cut,mag2mirrored);
  //ImageIO::Save("2_mag1mirrored",mag1mirrored);
  //ImageIO::Save("2_mag2mirrored",mag2mirrored);
  ImageIO::Save("2_mag1mirrored",mag1mirrored);
  ImageIO::Save("2_mag2mirrored",mag2mirrored);
  t4.Stop();

  //calculate Log-Polar spectrum for prepared magnitude
  t5.Start();
  mapperlp.Map(mag1mirrored, im1lp, MapBilinear);
  //if (ImageIO::Save("ImgLogPolar.mip", im1lp)!=0){
  if (ImageIO::Save("ImgLogPolar.mip", im1lp)!=0){
    BIASERR("error image");
  }

  mapperlp.Map(mag2mirrored, im2lp, MapBilinear);
  //if (ImageIO::Save("ImgLogPolar2.mip", im2lp)!=0){
  if (ImageIO::Save("ImgLogPolar2.mip", im2lp)!=0){
    BIASERR("error image");
  }
  t5.Stop();

  //calculate the Cross Power Spectrum for Log-Polar images
  t6.Start();
  fft.CrossPowerSpectrum(im1lp,im2lp,cps);
  //ImageIO::Save("3_CPS",cps);
  ImageIO::Save("3_CPS",cps);
  t6.Stop();

  //Transform reverse to find rotation and scale
  t7.Start();
  fft.TransformReverse(cps,res);
  //ImageIO::Save("4_Resolution",res);
  ImageIO::Save("4_Resolution",res);
  t7.Stop();

  //find the peak
  t8.Start();
  max = res.GetMaxPixelValue(0,coords);
  t8.Stop();

  //Correct the 2. source image with calculated rotation and scale
  t9.Start();
  angle = (double)coords[1] / (double) height * M_PI;
  if (angle > (0.5*M_PI)) angle -= M_PI;
  if(coords[0] > (unsigned int)(width/2))
    scale = pow((double)(width/2),((double)coords[0]-(double)width)/(double)width); 
  else
    scale = pow((double)(width/2),(double)coords[0]/(double) width);

  HTrans[0][2] = (double)(width/2);
  HTrans[1][2] = (double)(height/2);
  HTransInv[0][2] = (double)-(width/2);
  HTransInv[1][2] = (double)-(height/2);
  HRot[0][0]= cos(angle);   
  HRot[0][1]= -sin(angle);  
  HRot[1][0]= sin(angle);    
  HRot[1][1]= cos(angle); 
  HScale[0][0] = scale; 
  HScale[1][1] = scale;
  HScale[1][2] = 1.0;

  H =  HTrans * HScale * HRot * HTransInv;
  MapperHM.SetHomography(H);
  MapperHM.Map(im2grey, imgcorrected, MapBilinear);
  //ImageIO::Save("5_Corrected",imgcorrected);
  ImageIO::Save("5_Corrected",imgcorrected);
  t9.Stop();

  //Find the translation with the rotation and scale corrected image
  t11.Start();
  castfloat(imgcorrected,im2correctedfloat);
  fft.CrossPowerSpectrum(im1greyfloat,im2correctedfloat,cps);
  fft.TransformReverse(cps,res);
  //ImageIO::Save("6_Resolution2",res);
  ImageIO::Save("6_Resolution2",res);

  //Correct the image with the calculated translation
  max2 = res.GetMaxPixelValue(0,coords2);
  if(coords2[0] > (unsigned int)(width/2))
    HTrans[0][2] = width - (int)coords2[0];
  else
    HTrans[0][2] = -(int)coords2[0];
  if(coords2[1] > (unsigned int)(height/2))
    HTrans[1][2] = height - (int)coords2[1];
  else    
    HTrans[1][2] = -(int)coords2[1];

  MapperHM.SetHomography(HTrans);
  MapperHM.Map(imgcorrected, imgcorrected2, MapBilinear);
  //ImageIO::Save("7_Corrected2",imgcorrected2);
  ImageIO::Save("7_Corrected2",imgcorrected2);
  t11.Stop();  

  


  cout<<endl<<endl;
  cout <<"Loading Image and converting to grey took "<<t1.GetUserTime()<<" ms, "<<endl;
  cout <<"Floating took "<<t2.GetUserTime()<<" ms, "<<endl;
  cout <<"Init FFT took "<<t10.GetUserTime()<<" ms, "<<endl;
  cout <<"FFT took "<<t3.GetUserTime()<<" ms, "<<endl;
  cout <<"Mirror took "<<t4.GetUserTime()<<" ms, "<<endl;
  cout <<"LogPolar took: "<<t5.GetUserTime()<<" ms, "<<endl;
  cout <<"CPS took "<<t6.GetUserTime()<<" ms, "<<endl;
  cout <<"IFT took "<<t7.GetUserTime()<<" ms, "<<endl;
  cout <<"Finding Max took "<<t8.GetUserTime()<<" ms, "<<endl;
  cout <<"Correcting rotation and scale took "<<t9.GetUserTime()<<" ms, "<<endl;
  cout <<"Finding and correcting translation took "<<t11.GetUserTime()<<" ms, "<<endl;

  cout<<endl;

  cout<<"Max at: "<<coords[0]<<"/"<<coords[1]<<" value: "<<max<<endl<<endl;
  cout<<"Max2 at: "<<coords2[0]<<"/"<<coords2[1]<<" value: "<<max2<<endl<<endl;
  cout<<"Das entspricht (fuer ein 512x512 Bild):"<<endl;
  if(coords[0] > (unsigned int)(width/2))
    cout<<"Skalierungsfaktor: "<<pow(256.0,((double)coords[0]-width)/512.0)<<endl;
  else
  cout<<"Skalierungsfaktor: "<<pow(256.0,(double)coords[0]/512.0)<<endl;
  cout<<"Rotation: "<<180.0*angle/M_PI<<" Grad"<<endl;

  cout<<"Verschiebung: (";
  if(coords2[0] > (unsigned int)(width/2))
    cout<<width - (int)coords2[0];
  else
    cout<<-(int)coords2[0];
  cout<<",";
  if(coords2[1] > (unsigned int)(height/2))
    cout<<height - (int)coords2[1];
  else
    cout<<-(int)coords2[1];
  cout<<")"<<endl;
  return 0;
}