ExampleLineMatcher.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 ExampleLineMatcher.cpp
    @brief This little example demonstrates the usage of line matching
    @ingroup g_examples
    @author woelk 
*/

#include <Matcher2D/RegionMatcher.hh>
#include <Base/Image/Image.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Image/ImageConvert.hh>

//#include <Filter/Gradient.hh>
//#include <Filter/GradientSimple.hh>
#include <Filter/GradientSobel3x3.hh>

#include <Base/Math/Random.hh>
#include <Base/Debug/TimeMeasure.hh>

#include <iostream>
#include <iomanip>
#include <math.h>

using namespace BIAS;
using namespace std;



// storage type(s) for filter and image templates
//#define InST char
#define InST float
//#define OutST char
#define OutST float



int main(int argc, char*argv[])
{
    Image<unsigned char> inputUC;
    Image<InST> im1, im2;

    Image<OutST> grad1, grad2;
    unsigned int p1[2], start[2], end[2];
    unsigned int halfwinsize=7;
    unsigned int epsilon=2;
    double gradientscale=0.5;
    unsigned int result[2];
    double correlation=-1.0;
    unsigned int offset=epsilon+halfwinsize;
    int grad_threshold=2;
    int res;
    bool UseSAD=false;

    RegionMatcher matcher;
    //matcher.SetDebugLevel(D_REGION_MATCHER_LINE);
    //Gradient grad;
    //GradientSimple<int, int> grad;
    GradientSobel3x3<InST, OutST> grad;

    Random rand;
    double angle=rand.GetUniformDistributed(0.0, M_PI);
    double length1, length2;
    int count = 10000;
    int maxlength, maxlength2;
    bool debug=false;
    int debuglevel=0;
    TimeMeasure timerNCC, timerSAD;

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

    if (ImageIO::Load(argv[1], inputUC)!=0){
        BIASERR("error loading image "<<argv[1]);
        return -1;
    }
    //BIAS::ImageConvert::ConvertST(inputUC, im1, BIAS::ImageBase::ST_char);
    BIAS::ImageConvert::ConvertST(inputUC, im1, BIAS::ImageBase::ST_float);
    im1=inputUC;
    im2=im1;

    //grad.Sobel3Norm2(im1, grad1);
    grad.Filter(im1,grad1);
    grad2=grad1;

    //if (ImageIO::Save("grad", grad1)!=0){
    if (ImageIO::Save("grad", grad1)!=0){
        BIASERR("error writing image grad.mip");
        return -1;
    }

    for (int i=0; i< count; i++){

        do{
            p1[0]=rand.GetUniformDistributedInt((int)offset, 
                (int)im1.GetWidth()-1-(int)offset);
            p1[1]=rand.GetUniformDistributedInt((int)offset, 
                (int)im1.GetHeight()-1-(int)offset);
        } while (grad1.GetImageDataArray()[p1[1]][p1[0]]<grad_threshold );

        if (debuglevel>=2) 
            cout << i<<":\tp1: "<<p1[0]<<" "<<p1[1]<<endl;
        maxlength=
            ((int)p1[0]-(int)offset>(int)im1.GetWidth()-(int)offset-(int)p1[0]-1) ? 
            (int)im1.GetWidth()-(int)offset-(int)p1[0]-1 : (int)p1[0]-(int)offset ;
        maxlength2=
            ((int)p1[1]-(int)offset>(int)im1.GetHeight()-(int)offset-(int)p1[1]-1) ? 
            (int)im1.GetHeight()-(int)offset-(int)p1[1]-1 : (int)p1[1]-(int)offset ;
        maxlength=(maxlength>maxlength2)?maxlength2:maxlength;
        if (debuglevel>=2)       
            cout << "p1: "<<p1[0]<<" "<<p1[1]<<"\tmaxlength: "<<maxlength<<endl;
        length1=rand.GetUniformDistributedInt(0, maxlength);
        start[0]=(unsigned int)rint((int)p1[0]+sin(angle)*length1);
        start[1]=(unsigned int)rint((int)p1[1]+cos(angle)*length1);
        length2=rand.GetUniformDistributedInt(0, maxlength);
        end[0]=(unsigned int)rint((int)p1[0]-sin(angle)*length2);
        end[1]=(unsigned int)rint((int)p1[1]-cos(angle)*length2);

        if (debuglevel>=2)     
            cout << "searching between "
            <<start[0]<<" "<<start[1]<<" and "<<end[0]<<" "<<end[1]<<endl;

        int mycount=1;
        //if (UseSAD){
        timerSAD.Start();
        for (int k=0; k<mycount; k++){             
            if (matcher.LineMatcherSAD<OutST>(
                p1, start, end, 
                (const InST **)im1.GetImageDataArray(), 
                (const InST **)im2.GetImageDataArray(), 
                (const InST **)grad1.GetImageDataArray(),
                (const InST **)grad2.GetImageDataArray(), 
                halfwinsize, 
                epsilon, gradientscale, result, 
                correlation)
                !=0)
            {
                BIASERR("error in LineMatcherNCC");
                return -2;
            };
        }
        timerSAD.Stop();
        cout << "SAD took "<<timerSAD.GetUserTime()/mycount<<" ms"<<endl;
        // } else {
        timerNCC.Start();
        for (int k=0; k<mycount; k++){
            if (matcher.LineMatcherNCC(p1, start, end, 
                                       (const InST **)im1.GetImageDataArray(), 
                (const InST **)im2.GetImageDataArray(), 
                (const InST **)grad1.GetImageDataArray(),
                (const InST **)grad2.GetImageDataArray(), halfwinsize, 
                epsilon, gradientscale, result, 
                correlation)!=0){
                    BIASERR("error in LineMatcherNCC");
                    return -2;
                };
        }
        timerNCC.Stop();
        cout << "NCC took "<<timerNCC.GetUserTime()/mycount<<" ms"<<endl;
        //}

        if (debuglevel>=1)      
            cout << "correspondence to "<<p1[0]<<" "<<p1[1]<<" searched between "
            <<start[0]<<" "<<start[1]<<" and "<<end[0]<<" "<<end[1]<<" is at "
            <<result[0]<<" "<<result[1]<<" with correlation "<<correlation
            <<endl;

        if ((result[0]!=p1[0])||(result[1]!=p1[1])){
            cerr << "\n\n***********\n\nerror, found wrong result\n\n***********\n\n"<<endl;
            debug=true;
        } else {
            if (UseSAD){
                if (correlation>0.1) debug=true;
                else debug=false;
            } else {
                if (correlation<1.0) debug=true;
                else debug = false;
            }
        }
        if (debug){
            cerr << "staring debug: "<<endl;
            matcher.SetDebugLevel(D_REGION_MATCHER_LINE);
            if (UseSAD){
                if ((res =matcher.LineMatcherSAD(p1, start, end, 
                    (const InST **)im1.GetImageDataArray(), 
                    (const InST **)im2.GetImageDataArray(), 
                    (const InST **)grad1.GetImageDataArray(),
                    (const InST **)grad2.GetImageDataArray(), 
                    halfwinsize, 
                    epsilon, gradientscale, result, 
                    correlation))!=0){
                        BIASERR("error in LineMatcherSAD "<<res);
                    };
            } else {
                if ((res =matcher.LineMatcherNCC(p1, start, end, 
                    (const InST **)im1.GetImageDataArray(), 
                    (const InST **)im2.GetImageDataArray(), 
                    (const InST **)grad1.GetImageDataArray(),
                    (const InST **)grad2.GetImageDataArray(), 
                    halfwinsize, 
                    epsilon, gradientscale, result, 
                    correlation))!=0){
                        BIASERR("error in LineMatcherNCC "<<res);
                    };
            }

            matcher.SetDebugLevel(0);
        }

    }; // for 
    cout << "sum NCC took "<<timerNCC.GetUserTime()<<" ms"<<endl
        << "sum SAD took "<<timerSAD.GetUserTime()<<" ms"<<endl;

    return 0;
}