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

/**
   @example ExampleAlignment.cpp
   @relates ImageAlignment, TextureTransformHomography,TextureTransformAffine
   @brief Alignment example
   @ingroup g_examples
   @author woelk
*/

#include <Base/Image/Image.hh>
#include <Base/ImageUtils/ImageDraw.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Image/ImageConvert.hh>
#include <Base/Geometry/HomgPoint2D.hh>
#include <Base/Debug/TimeMeasure.hh>
#include <Image/TextureTransformHomography.hh>
#include <Image/TextureTransformSimilar.hh>
#include <Image/TextureTransformAffine.hh>
#include <Image/TextureTransformDisplacement.hh>
#include <Image/TextureTransformDisparity.hh>
#include <Base/Math/Random.hh>

#include <Matcher2D/ImageAlignment.hh>
#include <Image/TextureMapping.hh>
#include <iostream>
#include <iomanip>

using namespace BIAS;
using namespace std;

// loads an image and converts to single channel float image
int LoadImage(const string& name, Image<float>& im)
{
  ImageBase baseim;

  // load image
  if (ImageIO::Load(name, baseim)!=0){
    BIASERR("error loading image "<<name);
    return -1;
  }
  // convert to float
  if (ImageConvert::ConvertST(baseim, im, ImageBase::ST_float)!=0){
    BIASERR("error converting image "<<name);
    return -2;
  }
  // convert to grey if necessary
  if (im.GetChannelCount()!=1){
    if (ImageConvert::ToGrey(im, im)!=0){
      BIASERR("error converting image to grey "<<name);
      return -3;
    }
  }
  im.SetMetaData(*baseim.GetMetaData());

  return 0;
}


// the main function
int main(int argc, char*argv[])
{
  // read parameter from command line
  if (argc<2) {
    cerr << argv[0] << " <image> "<< endl;
    return -1;
  }
  // load first image
  Image<float> im[2];
  if (LoadImage(argv[1], im[0]) != 0) {
    BIASERR("error loading image " << argv[1]);
    return -1;
  }
  //  if (LoadImage(argv[2], im[1]) != 0) {
  //    BIASERR("error loading image " << argv[2]);
  //    return -1;
  //  }

  im[1] = im[0];
  //Gauss<float, float> smoother;
  //smoother.SetSigma(2.0);
  //smoother.Filter(im[1], im[0]);
  im[1].SetZero();

  // get the point to track
  HomgPoint2D p[2];
  p[0][2] = p[1][2] = 1.0;
  p[0][0] = im[0].GetWidth()/2.0;
  p[0][1] = im[0].GetHeight()/2.0;

  LocalAffineFrame LAF;
  Vector2<double> thepoint(p[0][0],p[0][1] );
  LAF.SetT(thepoint);
  LAF.SetA(Matrix<double>(2,2,MatrixIdentity)*(im[0].GetWidth()/16.0));

  Random R;
  Vector<double> par, origpar;
  Matrix<double> cov;
  TextureTransform *TT = NULL;
  int mode = 3;
  switch (mode) {
  case 4:TT = new TextureTransformHomography;
    cout<<"Testing HOMOGRAPHY"<<endl;
    par.newsize(6); par[0] = 0.1; par[1] = 0.1; par[2] = 0.1; par[3] = 0.1;
    par[4] = 1; par[5] = -1; par[6] = 0.01; par[7] = -0.01;
    break;
  case 3:TT = new TextureTransformAffine;
    cout<<"Testing AFFINITY"<<endl;
    par.newsize(6); par[0] = 0.2; par[1] = -0.1; par[2] = -0.2; par[3] = 0.3;
    par[4] = 1; par[5] = -1;
    par[0] = cos(0.4)-1.0; par[1] = sin(0.4); par[2] = -sin(0.4); par[3] = cos(0.4)-1.0;
    par[4] = -2; par[5] = 2;
    break;
  case 2:TT = new TextureTransformSimilar;
    cout<<"Testing SIMILARITY"<<endl;
    par.newsize(4); par[0] = -0.3; par[1] = 0.3; par[2] = 1; par[3] = -1; break;
  case 1:TT = new TextureTransformDisplacement;
    cout<<"Testing DISPLACEMENT"<<endl;
    par.newsize(2); par[0] = -1; par[1] = 1; break;
  default:TT = new TextureTransformDisparity;
    cout<<"Testing DISPARITY"<<endl;
    par.newsize(1); par[0]= 1; break;
  }
  cov.newsize(par.Size(), par.Size());
  cov.SetIdentity();

  cout<<"using parameters "<<par<<endl;
  origpar = par;
  TT->SetParameters(par);
  TextureTransformSimilar* pTTS = dynamic_cast<TextureTransformSimilar*>(TT);
  if (pTTS!=NULL) {
    Vector2<double> thepoint2(thepoint);
    //thepoint2[1] = 0;thepoint2[0] = 0;
    cout<<"setting origin !"<<thepoint2<<endl;
    pTTS->SetOrigin(thepoint2);
  }
  TextureTransformAffine* pTTA = dynamic_cast<TextureTransformAffine*>(TT);
  if (pTTA!=NULL) {
    Vector2<double> thepoint2(thepoint);
    //thepoint2[1] = 0;thepoint2[0] = 0;
    cout<<"setting origin !"<<thepoint2<<endl;
    pTTA->SetOrigin(thepoint2);
  }

  ImageAlignment IA;
  IA.SetTextureTransform(*TT);
  IA.AddDebugLevel(D_IMAGEALIGNMENT_PROGRESS);
  //IA.AddDebugLevel(D_IMAGEALIGNMENT_PARAMETER);
  //IA.AddDebugLevel(D_IMAGEALIGNMENT_PERPIXEL);
  IA.AddDebugLevel(D_IMAGEALIGNMENT_IMAGES);

  IA.SetMaxIterations(50);
  IA.SetUpdateThreshold(0.01);


  // map with random parameters
  TextureMapping<float, float> TM;

  Vector<double> invParams =  TT->GetInverseParameters();
  TT->SetParameters(invParams);

  im[1] = im[0];
  im[0].SetZero();
  TM.SetTextureTransform(*TT);
  TM.Map(im[1], im[0], MapBilinear);

  // change brightness of image 2
  bool brightnesschange = true;
  if (false) {//(brightnesschange) {
    float *pData = im[1].GetImageData();
    const float *pDataEnd = im[1].GetImageData()+im[1].GetPixelCount();
    do {
      *pData = 0.9f * *pData + 20.0f;
    } while (pData++ < pDataEnd);
  }
  IA.SetAffineBrightnessInvariance(brightnesschange);
  IA.SetDampening(1.0);
  // some noise on image 2
  double pixelnoise = 0.0;
  if (pixelnoise>0.0) {
    float *pData = im[1].GetImageData();
    const float *pDataEnd = im[1].GetImageData()+im[1].GetPixelCount();
    do {
      *pData = R.GetNormalDistributed(*pData, pixelnoise);
    } while (pData++ < pDataEnd);
  }
  //ImageIO::Save("im0.mip", im[0]);
  //ImageIO::Save("im1.mip", im[1]);
  ImageIO::Save("im0.mip", im[0]);
  ImageIO::Save("im1.mip", im[1]);

  // some noise on parameters to start estimate
  cov.SetIdentity();
  cov *= 1e-10;
  double parameternoise = 0.2;
  if (parameternoise>0.0) {
    for (int i=par.Size()-1; i>=0; i--)  {
      par[i] *= R.GetUniformDistributed(1.0-parameternoise, 1.0+parameternoise);
      cov[i][i] = par[i]*par[i]*parameternoise*parameternoise + 1e-10;
    }
  }

  // prepare alignment
  PyramidImage<float> i1,i2;
  i1.Init(im[0], 4);
  i2.Init(im[1], 4);

  cout<<"Setting alignment pixels"<<endl<<flush;
  IA.SetAlignmentPixelsRectangular(LAF,80);
  cout<<"Aligning ..."<<endl<<flush;

  Vector<double> startpar = par;
  // align images
  //int  result = IA.Align(i1, i2, par, cov, -1);
  int  result = //IA.AutoAlign(i1, i2, par, cov);
    IA.StrictPyramidAlign(i1, i2, par, cov);

  // show results
  double scale=0, offset=0;
  IA.GetLastBrightnessChange(scale, offset);
  cout << "result: " << result << ", "
       << p[0] << " -> " << par << ", "<<endl;
  cout<<"original parameters were "<<origpar<<endl;
  cout<<"optimization started at "<<startpar<<endl;
  cout<<"Brightness change is "<<std::setprecision(3)<<scale<<" "<<offset<<endl;
  return 0;
}