ExampleParabola.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003, 2004    (see file CONTACTS 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 ExampleParabola.cpp
    @brief This little example demonstrates the usage of parabola
    @ingroup g_examples
    @author woelk 
*/

#include <Matcher2D/CornerMatcher.hh>
#include <Base/Image/Image.hh>
#include <Base/ImageUtils/ImageDraw.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Geometry/HomgPoint2D.hh>
#include <Base/Debug/TimeMeasure.hh>
#include <Filter/Gradient.hh>

#include <iostream>
#include <iomanip>

using namespace BIAS;
using namespace std;

int main(int argc, char*argv[])
{
  Image<unsigned char> im[2], drawim;
  HomgPoint2D p[2], realpoint;

  double correlation, mincorrelation=1.0, error, errorsum=0.0;
  unsigned int halfsearchwinsize=2;
  unsigned int uip[2], rectsize=7;
  unsigned int halfnccwinsize=3;
  unsigned int newim, oldim, res;
  int counter;
  bool debug=false;
  string file;
  bool draw=false;
  bool UseLeastSquaresParabola=false;
  MetaData *md;
  AppData appdata;

  CornerMatcher matcher;
  TimeMeasure timer, tc;

  if (debug){
    matcher.SetDebugLevel(matcher.GetDebugLevel() |D_CORNER_MATCHER_SEARCH);
    matcher.SetDebugLevel(matcher.GetDebugLevel() |D_REGION_MATCHER_SEARCHNCC);
    matcher.SetDebugLevel(matcher.GetDebugLevel() |D_REGION_MATCHER_PARABOLA);
    //     matcher.SetDebugLevel(matcher.GetDebugLevel()|D_REGION_MATCHER_SEARCH_ODD);
    //     matcher.SetDebugLevel(matcher.GetDebugLevel()|D_REGION_MATCHER_SEARCH_EVEN);
  }

  if (argc<4){
    cerr << argv[0] << " : <x-coo> <y-coo> <image> <image> ..."<< endl;
    return -5;
  }

  p[0][0]=atof(argv[1]);
  p[0][1]=atof(argv[2]);
  p[0][2]=p[1][2]=1.0;

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

  counter=3;
  newim=0;
  while (counter < argc-1){
    counter++;
    newim=(newim==0)?(1):(0);
    oldim=(newim==0)?(1):(0);

    if (ImageIO::Load(argv[counter], im[newim])!=0){
      BIASERR("error loading image "<<argv[counter]);
      return -1;
    } else {
      if (debug)
        cout << "--------------- newim: " << argv[counter]
             << "  --- oldim: " << argv[counter-1] << endl;
      md=im[newim].GetMetaData();
      if (md->Find(AppData::MD_HomgPoint2D, "#[HomgPoint2D]", appdata)>=0){
        if (appdata.tag==AppData::MD_USE_ASCII){
          if (debug) cerr << appdata.stag <<" : "<<appdata.sdata<<endl;
          sscanf(appdata.sdata.c_str(), "[ %lf %lf %lf ]",
                 &realpoint[0], &realpoint[1], &realpoint[2]);
        } else {
          memcpy(realpoint.GetData(), appdata.data, appdata.length);
        }
        if (debug) cerr <<"realpoint is at "<<realpoint<<endl;
      } else {
        BIASERR("cannot find HomgPoint2D in meta data");
      }
    }

    if (p[oldim][0]<=halfnccwinsize+3 ||
        p[oldim][1]<=halfnccwinsize+3 ||
        p[oldim][0]+halfnccwinsize+3>=im[oldim].GetWidth()-1 ||
        p[oldim][1]+halfnccwinsize+3>=im[oldim].GetHeight()-1){
      cout << "point ("<<p[oldim][0]<<", "<<p[oldim][1]<<") too close to "
           <<"image border, stopping"<<endl;
      break;
    }

    timer.Start();
    tc.Start();
    if ((res=matcher.NCCSearch(p[oldim], p[oldim], im[oldim], im[newim],
                               halfnccwinsize, halfsearchwinsize, p[newim],
                               mincorrelation, correlation))!=0){
      BIASERR("error in NCCSearch");
    }

    if (res!=0){
      BIASERR("error matching");
      break;
    }

    if (debug)
      cerr << p[oldim] << " NCC -> " << p[newim]
           << " ("<<correlation<<")   Parabola -> ";

    if (UseLeastSquaresParabola){
      if ((res=matcher.ParabolaNCC5(p[oldim], p[newim], im[oldim], im[newim],
                                    halfnccwinsize, p[newim]))!=0){
        BIASERR("error in ParabolaNCC");
        break;
      }
    } else {
      if ((res=matcher.ParabolaNCC(p[oldim], p[newim], im[oldim], im[newim],
                                   halfnccwinsize, p[newim]))!=0){
        BIASERR("error in ParabolaNCC");
        break;
      }
    }
    tc.Stop();

    timer.Stop();

    tc.Print();
    tc.Reset();
    if (debug)
      cerr << p[newim] << endl;

    if (res!=0){
      BIASERR("error matching");
      break;
    }

    error=p[newim].Distance(realpoint);
    errorsum+=error;

    cout << setw(2) << counter-3 << " : (" << p[oldim][0]<<", "<<p[oldim][1]
         << ") -> ("<< p[newim][0]<<", "<<p[newim][1]<<")  corr: "<<correlation
         <<"   real: ("<<realpoint[0]<<", "<<realpoint[1]
         <<")   error: "<<error<<endl;

    if (draw){
      uip[0]=(unsigned int)rint(p[newim][0]);
      uip[1]=(unsigned int)rint(p[newim][1]);
      drawim=im[newim];
      ImageDraw<unsigned char>::RectangleCenter(drawim, uip, rectsize);

      file=argv[counter];
      file+=".tracked";
      //ImageIO::Save(file, drawim);
      ImageIO::Save(file, drawim);
    }
  }

  cerr <<"errorsum over "<<counter-3<<" images: "<<errorsum<<endl
       <<"  mean error "<<errorsum/(double)(counter-3)<<endl;

  cerr << "Parabola approximation over "<<argc-3<<" images took "
       <<timer.GetUserTime()<<" ms"
       << "\n  this is equal to "<<timer.GetUserTime()/double(argc-3)
       <<" ms per image and corner" << endl;

  return 0;
}