Basic Image AlgorithmS Library 2.8.0

EvaluateAlignment.cpp

00001 /*
00002 This file is part of the BIAS library (Basic ImageAlgorithmS).
00003 
00004 Copyright (C) 2003-2009    (see file CONTACT for details)
00005   Multimediale Systeme der Informationsverarbeitung
00006   Institut fuer Informatik
00007   Christian-Albrechts-Universitaet Kiel
00008 
00009 
00010 BIAS is free software; you can redistribute it and/or modify
00011 it under the terms of the GNU Lesser General Public License as published by
00012 the Free Software Foundation; either version 2.1 of the License, or
00013 (at your option) any later version.
00014 
00015 BIAS is distributed in the hope that it will be useful,
00016 but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 GNU Lesser General Public License for more details.
00019 
00020 You should have received a copy of the GNU Lesser General Public License
00021 along with BIAS; if not, write to the Free Software
00022 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00023 */
00024 
00025 
00026 /**
00027    @example EvaluateAlignment.cpp
00028    @relates ImageAlignment, TextureTransformHomography,TextureTransformAffine
00029    @brief Alignment evaluation
00030    @ingroup g_examples
00031    @author woelk
00032 */
00033 
00034 #include <Base/Image/Image.hh>
00035 #include <Base/ImageUtils/ImageDraw.hh>
00036 #include <Base/Image/ImageIO.hh>
00037 #include <Base/Image/ImageConvert.hh>
00038 #include <Base/Geometry/HomgPoint2D.hh>
00039 #include <Base/Debug/TimeMeasure.hh>
00040 #include <Image/TextureTransformHomography.hh>
00041 #include <Image/TextureTransformAffine.hh>
00042 #include <Base/Math/Random.hh>
00043 #include "MathAlgo/Lapack.hh"
00044 #include <Matcher2D/ImageAlignment.hh>
00045 #include <Image/HomographyMapping.hh>
00046 #include <Geometry/RMatrix.hh>
00047 #include <iostream>
00048 #include <iomanip>
00049 #include "Base/Common/FileHandling.hh"
00050 #include "Utils/ThreeDOut.hh"
00051 using namespace BIAS;
00052 using namespace std;
00053 
00054 
00055 // loads an image and converts to single channel float image
00056 int LoadImage(const string& name, Image<float>& im)
00057 {
00058   ImageBase baseim;
00059 
00060   // load image
00061   if (ImageIO::Load(name, baseim)!=0){
00062     BIASERR("error loading image "<<name);
00063     return -1;
00064   }
00065   // convert to float
00066   if (ImageConvert::ConvertST(baseim, im, ImageBase::ST_float)!=0){
00067     BIASERR("error converting image "<<name);
00068     return -2;
00069   }
00070   // convert to grey if necessary
00071   if (im.GetChannelCount()!=1){
00072     if (ImageConvert::ToGrey(im, im)!=0){
00073       BIASERR("error converting image to grey "<<name);
00074       return -3;
00075     }
00076   }
00077   im.SetMetaData(*baseim.GetMetaData());
00078 
00079   return 0;
00080 }
00081 
00082 
00083 // the main function
00084 int main(int argc, char*argv[])
00085 {
00086   // read parameter from command line
00087   if (argc<2) {
00088     cerr << argv[0] << " <image> "<< endl;
00089     return -1;
00090   }
00091   // load first image
00092   Image<float> im[2];
00093   if (LoadImage(argv[1], im[0]) != 0) {
00094     BIASERR("error loading image " << argv[0]);
00095     return -1;
00096   }
00097   im[1] = im[0];
00098   //  Gauss<float, float> smoother;
00099   // smoother.SetSigma(2.0);
00100   //smoother.Filter(im[1], im[0]);
00101   im[1].SetZero();
00102 
00103   Random R;
00104   Vector<double> par, origpar;
00105   Matrix<double> cov;
00106   TextureTransformAffine *TT = new TextureTransformAffine;
00107   cout<<"Testing AFFINITY"<<endl;
00108   par.newsize(6); par[0] = 0.0; par[1] = 0; par[2] = 0; par[3] = 0;
00109   par[4] = 0; par[5] = 0;
00110 
00111   cov.newsize(par.Size(), par.Size());
00112   cov.SetIdentity();
00113 
00114   cout<<"using parameters "<<par<<endl;
00115   origpar = par;
00116   TT->SetParameters(par);
00117 
00118   Vector2<double> thepoint((im[0].GetWidth()-1.0)/2.0,
00119                            (im[0].GetHeight()-1.0)/2.0);
00120 
00121 
00122 
00123 
00124   cout<<"setting origin !"<<thepoint<<endl;
00125   TT->SetOrigin(thepoint);
00126 
00127   // get the point to track
00128   HomgPoint2D p[2];
00129   p[0][2] = p[1][2] = 1.0;
00130   p[0][0] = im[0].GetWidth()/2.0;
00131   p[0][1] = im[0].GetHeight()/2.0;
00132 
00133   LocalAffineFrame LAF;
00134   LAF.SetT(thepoint);
00135   double d=20;// halfwinsize
00136   LAF.SetA(Matrix<double>(2,2,MatrixIdentity)*d);
00137 
00138 
00139   ImageAlignment IA;
00140   IA.SetTextureTransform(*TT);
00141   //IA.AddDebugLevel(D_IMAGEALIGNMENT_PROGRESS);
00142   //IA.AddDebugLevel(D_IMAGEALIGNMENT_PARAMETER);
00143   //IA.AddDebugLevel(D_IMAGEALIGNMENT_PERPIXEL);
00144   IA.AddDebugLevel(D_IMAGEALIGNMENT_IMAGES);
00145 
00146   IA.SetMaxIterations(50);
00147   IA.SetUpdateThreshold(0.001);
00148   IA.SetDampening(0.0001);
00149 
00150 
00151 
00152   KMatrix K(MatrixIdentity);
00153   K[1][1] =K[0][0] = im[0].GetWidth()*2.0;//2
00154    //K[1][1] = im[0].GetHeight()/2.0;
00155   K[0][2] = (im[0].GetWidth()-1.0)/2.0;
00156   K[1][2] = (im[0].GetHeight()-1.0)/2.0;
00157   cout << "K is "<<K<<endl;
00158   KMatrix Kinv;
00159   K.GetInverse(Kinv);
00160   //ImageIO::Save("im0.mip", im[0]);
00161   ImageIO::Save("im0.mip", im[0]);
00162 
00163   Matrix<double> PredictedCov(6,6,MatrixIdentity);
00164   for (unsigned int i=0;i<4;i++) PredictedCov[i][i] = 0.01;
00165   Matrix<double> EmpiricCov(6,6,MatrixIdentity);
00166   for (unsigned int i=0; i<4; i++)
00167     EmpiricCov[i][i] = 0.5/(d*d*2.0);
00168   for (unsigned int i=4; i<6; i++)
00169     EmpiricCov[i][i] = 0.5;
00170 
00171 
00172   ofstream logfile("trackingerror.log");
00173   logfile<<"# log written by "<<argv[0]<<" hws="<<d<<endl<<flush;
00174   logfile<<"# angle empiricdistance preddistance distance cornercovdistance result"<<endl<<flush;
00175   PMatrix P1(MatrixIdentity);
00176   P1 = K * P1;
00177   P1.InvalidateDecomposition();
00178   ThreeDOut vrml;
00179   vrml.AddPMatrix(P1,im[0].GetWidth(), im[0].GetHeight(), RGBAuc(0,255,0,255),0.1);
00180   vrml.AddPoint(Vector3<double>(1,0,1), RGBAuc(255,0,0,255));
00181   vrml.AddPoint(Vector3<double>(0,0,1), RGBAuc(255,255,0,255));
00182   for (unsigned int angle=0; angle<80; angle+=2) {
00183     cout<<"angle is now "<<angle<<" degree."<<endl;
00184     // prepare alignment
00185     PyramidImage<float> i1;
00186 
00187     i1.Init(im[0]);
00188     im[1].SetZero();
00189 
00190 
00191 
00192 
00193 
00194 
00195 
00196     RMatrix Rot;
00197     Vector3<double> axis(0,1,0);
00198     Rot.Set(axis,double(angle)/180.0*M_PI);
00199     cout<<"Rot is "<<Rot<<endl;
00200     Vector3<double> C(0,0,-1);
00201     C = Rot * C - C;
00202     cout <<"C is "<<C<<endl;
00203     Vector3<double> normal(0,0,1);
00204     Vector3<double> s(0,0,1);
00205     PMatrix P2(K,Rot,C);
00206     vrml.AddPMatrix(P2,im[0].GetWidth(), im[0].GetHeight(), RGBAuc(0,0,255,255),0.1);
00207     HMatrix H;
00208     H.SetIdentity();
00209     H.MapAcrossPlane(P2, P1, HomgPlane3D(0,0,1,-1));
00210     /*
00211     H = s.ScalarProduct(normal) * Rot  +
00212       C.OuterProduct(normal);
00213 
00214     H = K * H * Kinv;
00215     */
00216     HomographyMapping<float, float> HM;
00217     HM.SetHomography(H);
00218     HM.Map(im[0], im[1], MapAnisotropic, false);
00219 
00220 
00221     // some noise on image 2
00222     double pixelnoise = 0.0;
00223     if (pixelnoise>0.0) {
00224       float *pData = im[1].GetImageData();
00225       const float *pDataEnd = im[1].GetImageData()+im[1].GetPixelCount();
00226       do {
00227         *pData = R.GetNormalDistributed(*pData, pixelnoise);
00228       } while (pData++ < pDataEnd);
00229     }
00230     stringstream ss;
00231     ss<<"im1-"<<FileHandling::toString(angle)<<".mip";
00232 
00233     //ImageIO::Save(ss.str(), im[1]);
00234     ImageIO::Save(ss.str(), im[1]);
00235 
00236     PyramidImage<float> i2;
00237     i2.Init(im[1]);
00238 
00239     cout<<"Setting alignment pixels "<<endl<<flush;
00240     IA.SetAlignmentPixelsRectangular(LAF,2*d+1);
00241     cout<<"Aligning ..."<<endl<<flush;
00242 
00243 
00244     HMatrix HLin;
00245     H.GetLinearized(HomgPoint2D(thepoint)).GetInverse(HLin);
00246 
00247     //cout<<"Linearized is "<<HLin<<" again: "<<endl
00248     //    <<HLin.GetLinearized(HomgPoint2D(thepoint))<<endl;
00249 
00250     Vector<double> gtpar(6);
00251     par = gtpar = TT->ExtractParameters(HLin);
00252 
00253 
00254     // align images
00255     cov = EmpiricCov;
00256     int  result = IA.AutoAlign(i1, i2, par, cov);
00257     //cout<<"result is "<<par<<endl;
00258 
00259     if (result!=0) {
00260       BIASERR("alignment failed !!!"<<result);
00261     }
00262     cout<<"angle is "<<angle<<endl;
00263     if (angle%10==0) getchar();
00264     //int  result = //IA.AutoAlign(i1, i2, par, cov);
00265     //  IA.StrictPyramidAlign(i1, i2, par, cov);
00266 
00267     vector<BIAS::Vector<double> > cornerParams(5, Vector<double>(6));
00268     double linearizationerror = 0;
00269     for (unsigned int corner=0; corner<5; corner++) {
00270       HomgPoint2D cornerpoint;
00271       switch (corner) {
00272       case 0: cornerpoint = HomgPoint2D(thepoint[0]-d,thepoint[1]-d); break;
00273       case 1: cornerpoint = HomgPoint2D(thepoint[0]-d,thepoint[1]+d); break;
00274       case 2: cornerpoint = HomgPoint2D(thepoint[0]+d,thepoint[1]-d); break;
00275       case 3: cornerpoint = HomgPoint2D(thepoint[0]+d,thepoint[1]+d); break;
00276       default:;
00277       }
00278       linearizationerror +=
00279         H.ComputeLinearizationError(HomgPoint2D(thepoint),cornerpoint);
00280       HMatrix HLinTmp;
00281       H.GetLinearized(cornerpoint).GetInverse(HLinTmp);
00282       cornerParams[corner] = TT->ExtractParameters(HLinTmp);
00283       //cout<<"affine transform at corner is "<<cornerParams[corner]<<endl;
00284     }
00285 
00286     BIAS::Vector<double> meanCornerParam(6,0.0);
00287     for (unsigned int i=0; i<cornerParams.size(); i++) {
00288       meanCornerParam += cornerParams[i];
00289     }
00290     meanCornerParam *= 1.0/double(cornerParams.size());
00291     Matrix<double> cornerCov(6,6,MatrixZero);
00292     for (unsigned int i=0; i<cornerParams.size(); i++) {
00293       cornerCov +=
00294         BIAS::Vector<double>(meanCornerParam - cornerParams[i]).OuterProduct( BIAS::Vector<double>(meanCornerParam - cornerParams[i]));
00295     }
00296     cornerCov *= 1.0/double(cornerParams.size()-1.0);
00297     //cout<<"cornercov is "<<cornerCov<<endl;
00298     for (unsigned int i=0; i<6; i++) cornerCov[i][i] += 1e-8;
00299 
00300     linearizationerror /= 4.0;
00301     cout<<"average linearization error is "<<linearizationerror<<endl;
00302     PredictedCov.SetIdentity();
00303     PredictedCov *= 1e-8;
00304     for (unsigned int i=0; i<4; i++)
00305       PredictedCov[i][i] += linearizationerror*linearizationerror/(d*d*2.0);
00306     for (unsigned int i=4; i<6; i++)
00307       PredictedCov[i][i] += linearizationerror*linearizationerror;
00308 
00309     gtpar -= par;
00310 
00311 
00312     double empiricdistance = MahalanobisDistance(EmpiricCov,gtpar);
00313     cout<<"empiric distance is "<<empiricdistance<<endl;
00314     double preddistance = MahalanobisDistance(PredictedCov,gtpar);
00315     cout<<"compensated distance is "<<preddistance<<endl;
00316     for (unsigned int i=0; i<6; i++) cov[i][i] += 1e-8;
00317     double distance = MahalanobisDistance(cov,gtpar);
00318     cout<<"obtained data distance is "<<distance<<endl;
00319     double cornercovdistance = MahalanobisDistance(cornerCov,gtpar);
00320     cout<<"corner cov distance is "<<cornercovdistance<<endl;
00321     if (result>=0) {
00322       logfile<<angle<<" "<<empiricdistance<<" "<<preddistance<<" "<<distance<<" "<<cornercovdistance<<" "<<result<<endl<<flush;
00323 
00324     }
00325 
00326   }
00327   vrml.VRMLOut("scene.wrl");
00328   return 0;
00329 }
00330 
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends