ExampleHomographyMapping.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 ExampleHomographyMapping.cpp
   @relates HomographyMapping
   @brief Example for HomographyMapping usage.
   @ingroup g_examples
   @author  MIP
*/

#include <Image/HomographyMapping.hh>
#include <Image/AffineMapping.hh>
#include <Image/DisplacementMapping.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Debug/TimeMeasure.hh>
#include <Base/Math/Random.hh>

using namespace BIAS;
using namespace std;

// #define AFFINE

int main(int argc, char *argv[])
{
  cout<<endl<<"Usage: "<<argv[0]
      <<" [testimage [ width height [homography.mat]]]"<<endl<<endl;
  cout<<"If you dont provide a test image, I will generate one :-)"
      <<endl<<endl;
  Image<unsigned char> im(640, 480, 1);
  if (argc>=2) {
    if (ImageIO::Load(argv[1], im)!=0){
      BIASERR("error loading "<<argv[1]);
      return -2;
    }

  }
  HMatrix H(MatrixIdentity), tmp;
  if (argc>4) {
    cout<<"Reading HMatrix "<<argv[4]<<endl;
    H.Load(argv[4]);
    cout<<"H is "<<H<<endl;
    H.GetInverse(tmp);
    H = tmp;
  } else {
    /*
    H[0][0]= 4.5 ;   H[0][1]= -0.08;      H[0][2] = 0;
    H[1][0]= 0.08;    H[1][1]= 4.5;     H[1][2] = 0;
    */

    // direction of enlargement
    double angle = -45.0 / 180.0 * M_PI;
    // outpuit direction
    double angle2 = 3.0 / 180.0 * M_PI;
    HMatrix HRot(MatrixIdentity), HRot2(MatrixIdentity);
    HRot[0][0]= cos(angle);   HRot[0][1] = sin(angle);
    HRot[1][0]= -sin(angle);    HRot[1][1]= cos(angle);


    HRot2[0][0]= cos(angle2);   HRot2[0][1] = sin(angle2);
    HRot2[1][0]= -sin(angle2);    HRot2[1][1]= cos(angle2);



    HMatrix HScale(MatrixIdentity);
    HScale[0][0] = 4.0; HScale[1][1] = 1.0;
    //HScale[0][0] = 2.0; HScale[1][1] = 2.0;
    // HScale[1][2] = -400;
    H = HRot2 * HRot.Transpose() * HScale * HRot;
    //H = HMatrix(MatrixIdentity);

    cout<<"Using default homography "<<H<<endl;
  }

  if (argc<2){
    cout<<"Generating test image orig.mip"<<endl;
    bool highFreqCross = true;
    HomgPoint2D point(0.0, 0.0, 1.0);
    HMatrix Hinv;
    H.GetInverse(Hinv);
    Hinv *= 1.0 / H[2][2];
    Random R;
    if (!im.IsEmpty()) im.Release();


    BackwardMapping<unsigned char, unsigned char>::
      GenerateTestImage(im,highFreqCross,5,245,Hinv);

    ImageIO::Save("orig.mip", im);





  } else {
    if (ImageIO::Load(argv[1], im)!=0){
      BIASERR("error loading "<<argv[1]);
      return -2;
    }
  }
  int width  = im.GetWidth()/2;
  int height = im.GetHeight()/2;

  if (argc>3) {
    width = atoi(argv[2]);
    height = atoi(argv[3]);
    cout<<"Using imagesize "<<width<<" x "<<height<<endl;
  } else {
    cout<<"Using input image size "<<width<<" x "<<height<<endl;
  }

  HMatrix HMove(MatrixIdentity), HMove2(MatrixIdentity), HMoveInv;
  HMove[0][2]= double(im.GetWidth())/-2.0;
  HMove[1][2] =double(im.GetHeight())/-2.0;
  HMove2[0][2]= double(width)/-2.0;   HMove2[1][2] =  double(height)/-2.0;
  HMove.GetInverse(HMoveInv);
  H = HMoveInv * H * HMove2;
  cout<<"ROI of source image is "<<*im.GetROI()
      <<" while image size is "<<width<<" x "<<height<<endl;


#ifdef AFFINE
  AffineMapping<unsigned char, unsigned char> Mapper;
#else
  HomographyMapping<unsigned char, unsigned char> Mapper, Mapper2;
#endif
  DisplacementMapping<unsigned char, unsigned char> Mapper_dis;

  // Mapper.SetPyramidSize(1);


  Image<unsigned char>
    imresnearest          (width, height, im.GetChannelCount()),
    imresbilinear         (width, height, im.GetChannelCount()),
    imrestrilinear        (width, height, im.GetChannelCount()),

    imresnearest_dismap   (width, height, im.GetChannelCount()),

    imresnearest_lookup   (width, height, im.GetChannelCount()),
    imresbilinear_lookup  (width, height, im.GetChannelCount()),
    imrestrilinear_lookup (width, height, im.GetChannelCount()),

    imresbicubic          (width, height, im.GetChannelCount()),
    imresanisotropic      (width, height, im.GetChannelCount()),
    imresdirect           (width, height, im.GetChannelCount());

    unsigned char *color = new unsigned char[im.GetChannelCount()];
    for(unsigned int i=0;i<im.GetChannelCount();i++) color[i]=(unsigned char)255;

    imresnearest.FillImageWithConstValue(color);
    imresbilinear.FillImageWithConstValue(color);
    imrestrilinear.FillImageWithConstValue(color);
    imresnearest_dismap.FillImageWithConstValue(color);
    imresnearest_lookup.FillImageWithConstValue(color);
    imresbilinear_lookup.FillImageWithConstValue(color);
    imrestrilinear_lookup.FillImageWithConstValue(color);
    imresbicubic.FillImageWithConstValue(color);
    imresanisotropic.FillImageWithConstValue(color);
    imresdirect.FillImageWithConstValue(color);

    delete[] color;
  //Nearest-Neigbour Displacement Map
  Image<float>
    NN_dismap             (width, height, 3);

#ifdef AFFINE
  Mapper.SetAffineTransformation(H[0][0],H[0][1],H[1][0], H[1][1],
                                 H[0][2],H[1][2]);
#else
  // trilinear mapping with automatic scale and pyramid size selection
  Mapper.SetHomography(H);

#endif


#ifdef AFFINE
  Mapper.MapDirect(im, imresdirect);
#endif


  //Mapper.GetDisplacementMap(NN_dismap, im);
  Mapper.GetDisplacementMap(NN_dismap, im.GetWidth(), im.GetHeight());
  Mapper_dis.SetDisplacementMap(NN_dismap);

  unsigned char* Data = imresnearest_dismap.GetImageData();

  int counter = 0;
  int counter_index = 0;
  double bias_x, bias_y;

  for(unsigned int j = 0; j < imresnearest_dismap.GetHeight(); j++){
    for(unsigned int i = 0; i < imresnearest_dismap.GetWidth(); i++){
      HomgPoint2D p_sink, p_src;
      p_sink[0] = i;
      p_sink[1] = j;
      p_sink[2] = 1.0;

      if(Mapper_dis.GetSourceCoordinates_(p_sink, p_src) != -1){


        im.TextureToBIASCoordinates(p_src[0], p_src[1], bias_x, bias_y);

        //BIG DEBUG CHECKING IF BIAS2TEXTURE and TEXTURE2BIAS work correctly
        /**
        HomgPoint2D p_src_p, p_bias;
        im.BIASToTextureCoordinates(bias_x, bias_y, p_src_p[0], p_src_p[1]);
        im.TextureToBIASCoordinates(p_src_p[0], p_src_p[1],
                                    p_bias[0], p_bias[1]);

        if(p_src[0] != p_src_p[0] || p_src[1] != p_src_p[1] ||
           bias_x != p_bias[0] || bias_y != p_bias[1]
           ){
          cout << endl;
          cout << "TEXT coords: " << p_src[0] << ", " << p_src[1] << endl;
          cout << "BIAS coords: " << bias_x << ", " << bias_y << endl;
          cout << "TEXT coords: " << p_src_p[0] << ", " << p_src_p[1] << endl;
          cout << "BIAS coords: " << p_bias[0] << ", " << p_bias[1] << endl;
        }
        **/

        int p_x = (int)rint(bias_x);
        int p_y = (int)rint(bias_y);
        int p_z = (int)(p_src[2]);

        if(p_z != -1 &&
           p_x >= 0 &&
           p_x < (int)im.GetWidth() &&
           p_y >= 0 &&
           p_y < (int)im.GetHeight()){
          counter_index = im.GetWidth()*p_y + p_x;
          *Data++ =(unsigned char)im.GetImageData()[counter_index];
        }
        else
          *Data++ = (unsigned char)(0);
      }
      counter++;
    }
  }


  cout<<"timing ..."<<endl<<flush;
  unsigned int numbermeasurements = 1;
  TimeMeasure t1,  t2,  t3,  t4,  t5,  t6;
  TimeMeasure t1l_init, t2l_init, t4l_init;//, t4l, t5l;
  TimeMeasure t1l, t2l, t4l;//, t4l, t5l;
  for (unsigned int i=0; i<numbermeasurements; i++) {

    //---NEAREST NEIGHBOUR---
    t1.Start();
    // NN mapping
    Mapper.Map(im, imresnearest, MapNearestNeighbor);
    t1.Stop();

    t1l_init.Start();
    Mapper.PrepareLookupTableMapping(im, imresnearest_lookup,
                                     MapNearestNeighbor);
    t1l_init.Stop();
    t1l.Start();
    // NN mapping
    Mapper.MapWithLookupTable(im, imresnearest_lookup, MapNearestNeighbor);
    t1l.Stop();

    //-----------------------

    //--------BILINEAR-------
    t2.Start();
    // fast backward mapping (bilinaer) disregarding sampling theroem
    Mapper.Map(im, imresbilinear, MapBilinear);
    t2.Stop();

    t2l_init.Start();
    Mapper.PrepareLookupTableMapping(im, imresbilinear_lookup,
                                     MapBilinear);

    t2l_init.Stop();
    t2l.Start();
    // trilinear mapping with automatic scale and pyramid size selection
    Mapper.MapWithLookupTable(im, imresbilinear_lookup, MapBilinear);
    t2l.Stop();
    //-----------------------
#ifdef TIMING
    BIASWARN("Anisoptropic filtering causes double free or corruption error.")
    t3.Start();
    // anisotropic
    Mapper.Map(im, imresanisotropic, MapAnisotropic);
    t3.Stop();
#endif //TIMING
    //-------TRILINEAR-------
    t4.Start();
    // trilinear mapping with automatic scale and pyramid size selection
    Mapper.Map(im, imrestrilinear, MapTrilinear, false,2.0);
    t4.Stop();

    t4l_init.Start();
    Mapper.PrepareLookupTableMapping(im, imrestrilinear_lookup,
                                     MapTrilinear);
    t4l_init.Stop();
    t4l.Start();
    // trilinear mapping with automatic scale and pyramid size selection
    Mapper.MapWithLookupTable(im, imrestrilinear_lookup, MapTrilinear);
    t4l.Stop();
    //-----------------------

    t5.Start();
    // bicubic mapping with automatic scale and pyramid size selection
    Mapper.Map(im, imresbicubic, MapBicubic);
    t5.Stop();
#ifdef AFFINE
    t6.Start();
    // direct affine mapping with automatic scale and pyramid size selection
    Mapper.MapDirect(im, imresdirect);
    t6.Stop();

#endif

  }
  cout <<"Anisotropic mapping took "
       <<t3.GetUserTime()/double(numbermeasurements)
       <<" ms, "<<endl

       <<"Bicubic mapping took "
       <<t5.GetUserTime()/double(numbermeasurements)
       <<" ms,"<<endl

       <<"NN-mapping (no lookups) took "
       <<t1.GetUserTime()/double(numbermeasurements)
       <<" ms"<< endl
       <<"NN-mapping (wt lookups) init took "
       <<t1l_init.GetUserTime()/double(numbermeasurements)
       <<" ms," << endl
       <<"NN-mapping (wt lookups) took "
       <<t1l.GetUserTime()/double(numbermeasurements)
       <<" ms," << endl

       <<"Bilinear mapping (no lookups) took "
       <<t2.GetUserTime()/double(numbermeasurements)
       <<" ms,"<<endl
       <<"Bilinear mapping (wt lookups) init took "
       <<t2l_init.GetUserTime()/double(numbermeasurements)
       <<" ms,"<<endl
       <<"Bilinear mapping (wt lookups) took "
       <<t2l.GetUserTime()/double(numbermeasurements)
       <<" ms,"<<endl

       <<"Trilinear mapping (no lookups) took "
       <<t4.GetUserTime()/double(numbermeasurements)
       <<" ms,"<<endl
       <<"Trilinear mapping (wt lookups) init took "
       <<t4l_init.GetUserTime()/double(numbermeasurements)
       <<" ms,"<<endl
       <<"Trilinear mapping (wt lookups) took "
       <<t4l.GetUserTime()/double(numbermeasurements)
       <<" ms,"<<endl
#ifdef AFFINE
       <<"Direct affine took "
       <<t6.GetUserTime()/double(numbermeasurements)<<" ms"
#endif
       <<endl;

  //if (ImageIO::Save("image_homography_anisotropic.mip", imresanisotropic)!=0){
  if (ImageIO::Save("image_homography_anisotropic.mip", imresanisotropic)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_trilinear.mip", imrestrilinear)!=0){
  if (ImageIO::Save("image_homography_trilinear.mip", imrestrilinear)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_trilinear_lookup.mip",
  //                   imrestrilinear_lookup)!=0){
  if (ImageIO::Save("image_homography_trilinear_lookup.mip",
                    imrestrilinear_lookup)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_bicubic.mip", imresbicubic)!=0){
  if (ImageIO::Save("image_homography_bicubic.mip", imresbicubic)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_bilinear.mip", imresbilinear)!=0){
  if (ImageIO::Save("image_homography_bilinear.mip", imresbilinear)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_bilinear_lookup.mip",
  //                   imresbilinear_lookup)!=0){
  if (ImageIO::Save("image_homography_bilinear_lookup.mip",
                     imresbilinear_lookup)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_nearest.mip", imresnearest)!=0){
  if (ImageIO::Save("image_homography_nearest.mip", imresnearest)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_nearest_lookup.mip",
  //                     imresnearest_lookup)!=0){
  if (ImageIO::Save("image_homography_nearest_lookup.mip",
                     imresnearest_lookup)!=0){
    BIASERR("error image");
    return -3;
  }
  //if (ImageIO::Save("image_homography_nearest_dismap.mip",
  //                   imresnearest_dismap)!=0){
  if (ImageIO::Save("image_homography_nearest_dismap.mip",
                     imresnearest_dismap)!=0){
    BIASERR("error image");
    return -3;
  }
#ifdef AFFINE
  if (ImageIO::Save("image_direct.mip", imresdirect)!=0){
    BIASERR("error image");
    return -3;
  }
#endif
  return 0;
}