Basic Image AlgorithmS Library  2.8.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
EvaluateAlignment.cpp
1 /*
2 This file is part of the BIAS library (Basic ImageAlgorithmS).
3 
4 Copyright (C) 2003-2009 (see file CONTACT for details)
5  Multimediale Systeme der Informationsverarbeitung
6  Institut fuer Informatik
7  Christian-Albrechts-Universitaet Kiel
8 
9 
10 BIAS is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 2.1 of the License, or
13 (at your option) any later version.
14 
15 BIAS is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU Lesser General Public License for more details.
19 
20 You should have received a copy of the GNU Lesser General Public License
21 along with BIAS; if not, write to the Free Software
22 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 */
24 
25 
26 /**
27  @example EvaluateAlignment.cpp
28  @relates ImageAlignment, TextureTransformHomography,TextureTransformAffine
29  @brief Alignment evaluation
30  @ingroup g_examples
31  @author woelk
32 */
33 
34 #include <Base/Image/Image.hh>
35 #include <Base/ImageUtils/ImageDraw.hh>
36 #include <Base/Image/ImageIO.hh>
37 #include <Base/Image/ImageConvert.hh>
38 #include <Base/Geometry/HomgPoint2D.hh>
39 #include <Base/Debug/TimeMeasure.hh>
40 #include <Image/TextureTransformHomography.hh>
41 #include <Image/TextureTransformAffine.hh>
42 #include <Base/Math/Random.hh>
43 #include "MathAlgo/Lapack.hh"
44 #include <Matcher2D/ImageAlignment.hh>
45 #include <Image/HomographyMapping.hh>
46 #include <Geometry/RMatrix.hh>
47 #include <iostream>
48 #include <iomanip>
49 #include "Base/Common/FileHandling.hh"
50 #include "Utils/ThreeDOut.hh"
51 using namespace BIAS;
52 using namespace std;
53 
54 
55 // loads an image and converts to single channel float image
56 int LoadImage(const string& name, Image<float>& im)
57 {
58  ImageBase baseim;
59 
60  // load image
61  if (ImageIO::Load(name, baseim)!=0){
62  BIASERR("error loading image "<<name);
63  return -1;
64  }
65  // convert to float
66  if (ImageConvert::ConvertST(baseim, im, ImageBase::ST_float)!=0){
67  BIASERR("error converting image "<<name);
68  return -2;
69  }
70  // convert to grey if necessary
71  if (im.GetChannelCount()!=1){
72  if (ImageConvert::ToGrey(im, im)!=0){
73  BIASERR("error converting image to grey "<<name);
74  return -3;
75  }
76  }
77  im.SetMetaData(*baseim.GetMetaData());
78 
79  return 0;
80 }
81 
82 
83 // the main function
84 int main(int argc, char*argv[])
85 {
86  // read parameter from command line
87  if (argc<2) {
88  cerr << argv[0] << " <image> "<< endl;
89  return -1;
90  }
91  // load first image
92  Image<float> im[2];
93  if (LoadImage(argv[1], im[0]) != 0) {
94  BIASERR("error loading image " << argv[0]);
95  return -1;
96  }
97  im[1] = im[0];
98  // Gauss<float, float> smoother;
99  // smoother.SetSigma(2.0);
100  //smoother.Filter(im[1], im[0]);
101  im[1].SetZero();
102 
103  Random R;
104  Vector<double> par, origpar;
105  Matrix<double> cov;
107  cout<<"Testing AFFINITY"<<endl;
108  par.newsize(6); par[0] = 0.0; par[1] = 0; par[2] = 0; par[3] = 0;
109  par[4] = 0; par[5] = 0;
110 
111  cov.newsize(par.Size(), par.Size());
112  cov.SetIdentity();
113 
114  cout<<"using parameters "<<par<<endl;
115  origpar = par;
116  TT->SetParameters(par);
117 
118  Vector2<double> thepoint((im[0].GetWidth()-1.0)/2.0,
119  (im[0].GetHeight()-1.0)/2.0);
120 
121 
122 
123 
124  cout<<"setting origin !"<<thepoint<<endl;
125  TT->SetOrigin(thepoint);
126 
127  // get the point to track
128  HomgPoint2D p[2];
129  p[0][2] = p[1][2] = 1.0;
130  p[0][0] = im[0].GetWidth()/2.0;
131  p[0][1] = im[0].GetHeight()/2.0;
132 
133  LocalAffineFrame LAF;
134  LAF.SetT(thepoint);
135  double d=20;// halfwinsize
136  LAF.SetA(Matrix<double>(2,2,MatrixIdentity)*d);
137 
138 
139  ImageAlignment IA;
140  IA.SetTextureTransform(*TT);
141  //IA.AddDebugLevel(D_IMAGEALIGNMENT_PROGRESS);
142  //IA.AddDebugLevel(D_IMAGEALIGNMENT_PARAMETER);
143  //IA.AddDebugLevel(D_IMAGEALIGNMENT_PERPIXEL);
144  IA.AddDebugLevel(D_IMAGEALIGNMENT_IMAGES);
145 
146  IA.SetMaxIterations(50);
147  IA.SetUpdateThreshold(0.001);
148  IA.SetDampening(0.0001);
149 
150 
151 
153  K[1][1] =K[0][0] = im[0].GetWidth()*2.0;//2
154  //K[1][1] = im[0].GetHeight()/2.0;
155  K[0][2] = (im[0].GetWidth()-1.0)/2.0;
156  K[1][2] = (im[0].GetHeight()-1.0)/2.0;
157  cout << "K is "<<K<<endl;
158  KMatrix Kinv;
159  K.GetInverse(Kinv);
160  //ImageIO::Save("im0.mip", im[0]);
161  ImageIO::Save("im0.mip", im[0]);
162 
163  Matrix<double> PredictedCov(6,6,MatrixIdentity);
164  for (unsigned int i=0;i<4;i++) PredictedCov[i][i] = 0.01;
165  Matrix<double> EmpiricCov(6,6,MatrixIdentity);
166  for (unsigned int i=0; i<4; i++)
167  EmpiricCov[i][i] = 0.5/(d*d*2.0);
168  for (unsigned int i=4; i<6; i++)
169  EmpiricCov[i][i] = 0.5;
170 
171 
172  ofstream logfile("trackingerror.log");
173  logfile<<"# log written by "<<argv[0]<<" hws="<<d<<endl<<flush;
174  logfile<<"# angle empiricdistance preddistance distance cornercovdistance result"<<endl<<flush;
176  P1 = K * P1;
177  P1.InvalidateDecomposition();
178  ThreeDOut vrml;
179  vrml.AddPMatrix(P1,im[0].GetWidth(), im[0].GetHeight(), RGBAuc(0,255,0,255),0.1);
180  vrml.AddPoint(Vector3<double>(1,0,1), RGBAuc(255,0,0,255));
181  vrml.AddPoint(Vector3<double>(0,0,1), RGBAuc(255,255,0,255));
182  for (unsigned int angle=0; angle<80; angle+=2) {
183  cout<<"angle is now "<<angle<<" degree."<<endl;
184  // prepare alignment
186 
187  i1.Init(im[0]);
188  im[1].SetZero();
189 
190 
191 
192 
193 
194 
195 
196  RMatrix Rot;
197  Vector3<double> axis(0,1,0);
198  Rot.Set(axis,double(angle)/180.0*M_PI);
199  cout<<"Rot is "<<Rot<<endl;
200  Vector3<double> C(0,0,-1);
201  C = Rot * C - C;
202  cout <<"C is "<<C<<endl;
203  Vector3<double> normal(0,0,1);
204  Vector3<double> s(0,0,1);
205  PMatrix P2(K,Rot,C);
206  vrml.AddPMatrix(P2,im[0].GetWidth(), im[0].GetHeight(), RGBAuc(0,0,255,255),0.1);
207  HMatrix H;
208  H.SetIdentity();
209  H.MapAcrossPlane(P2, P1, HomgPlane3D(0,0,1,-1));
210  /*
211  H = s.ScalarProduct(normal) * Rot +
212  C.OuterProduct(normal);
213 
214  H = K * H * Kinv;
215  */
217  HM.SetHomography(H);
218  HM.Map(im[0], im[1], MapAnisotropic, false);
219 
220 
221  // some noise on image 2
222  double pixelnoise = 0.0;
223  if (pixelnoise>0.0) {
224  float *pData = im[1].GetImageData();
225  const float *pDataEnd = im[1].GetImageData()+im[1].GetPixelCount();
226  do {
227  *pData = R.GetNormalDistributed(*pData, pixelnoise);
228  } while (pData++ < pDataEnd);
229  }
230  stringstream ss;
231  ss<<"im1-"<<FileHandling::toString(angle)<<".mip";
232 
233  //ImageIO::Save(ss.str(), im[1]);
234  ImageIO::Save(ss.str(), im[1]);
235 
237  i2.Init(im[1]);
238 
239  cout<<"Setting alignment pixels "<<endl<<flush;
240  IA.SetAlignmentPixelsRectangular(LAF,2*d+1);
241  cout<<"Aligning ..."<<endl<<flush;
242 
243 
244  HMatrix HLin;
245  H.GetLinearized(HomgPoint2D(thepoint)).GetInverse(HLin);
246 
247  //cout<<"Linearized is "<<HLin<<" again: "<<endl
248  // <<HLin.GetLinearized(HomgPoint2D(thepoint))<<endl;
249 
250  Vector<double> gtpar(6);
251  par = gtpar = TT->ExtractParameters(HLin);
252 
253 
254  // align images
255  cov = EmpiricCov;
256  int result = IA.AutoAlign(i1, i2, par, cov);
257  //cout<<"result is "<<par<<endl;
258 
259  if (result!=0) {
260  BIASERR("alignment failed !!!"<<result);
261  }
262  cout<<"angle is "<<angle<<endl;
263  if (angle%10==0) getchar();
264  //int result = //IA.AutoAlign(i1, i2, par, cov);
265  // IA.StrictPyramidAlign(i1, i2, par, cov);
266 
267  vector<BIAS::Vector<double> > cornerParams(5, Vector<double>(6));
268  double linearizationerror = 0;
269  for (unsigned int corner=0; corner<5; corner++) {
270  HomgPoint2D cornerpoint;
271  switch (corner) {
272  case 0: cornerpoint = HomgPoint2D(thepoint[0]-d,thepoint[1]-d); break;
273  case 1: cornerpoint = HomgPoint2D(thepoint[0]-d,thepoint[1]+d); break;
274  case 2: cornerpoint = HomgPoint2D(thepoint[0]+d,thepoint[1]-d); break;
275  case 3: cornerpoint = HomgPoint2D(thepoint[0]+d,thepoint[1]+d); break;
276  default:;
277  }
278  linearizationerror +=
279  H.ComputeLinearizationError(HomgPoint2D(thepoint),cornerpoint);
280  HMatrix HLinTmp;
281  H.GetLinearized(cornerpoint).GetInverse(HLinTmp);
282  cornerParams[corner] = TT->ExtractParameters(HLinTmp);
283  //cout<<"affine transform at corner is "<<cornerParams[corner]<<endl;
284  }
285 
286  BIAS::Vector<double> meanCornerParam(6,0.0);
287  for (unsigned int i=0; i<cornerParams.size(); i++) {
288  meanCornerParam += cornerParams[i];
289  }
290  meanCornerParam *= 1.0/double(cornerParams.size());
291  Matrix<double> cornerCov(6,6,MatrixZero);
292  for (unsigned int i=0; i<cornerParams.size(); i++) {
293  cornerCov +=
294  BIAS::Vector<double>(meanCornerParam - cornerParams[i]).OuterProduct( BIAS::Vector<double>(meanCornerParam - cornerParams[i]));
295  }
296  cornerCov *= 1.0/double(cornerParams.size()-1.0);
297  //cout<<"cornercov is "<<cornerCov<<endl;
298  for (unsigned int i=0; i<6; i++) cornerCov[i][i] += 1e-8;
299 
300  linearizationerror /= 4.0;
301  cout<<"average linearization error is "<<linearizationerror<<endl;
302  PredictedCov.SetIdentity();
303  PredictedCov *= 1e-8;
304  for (unsigned int i=0; i<4; i++)
305  PredictedCov[i][i] += linearizationerror*linearizationerror/(d*d*2.0);
306  for (unsigned int i=4; i<6; i++)
307  PredictedCov[i][i] += linearizationerror*linearizationerror;
308 
309  gtpar -= par;
310 
311 
312  double empiricdistance = MahalanobisDistance(EmpiricCov,gtpar);
313  cout<<"empiric distance is "<<empiricdistance<<endl;
314  double preddistance = MahalanobisDistance(PredictedCov,gtpar);
315  cout<<"compensated distance is "<<preddistance<<endl;
316  for (unsigned int i=0; i<6; i++) cov[i][i] += 1e-8;
317  double distance = MahalanobisDistance(cov,gtpar);
318  cout<<"obtained data distance is "<<distance<<endl;
319  double cornercovdistance = MahalanobisDistance(cornerCov,gtpar);
320  cout<<"corner cov distance is "<<cornercovdistance<<endl;
321  if (result>=0) {
322  logfile<<angle<<" "<<empiricdistance<<" "<<preddistance<<" "<<distance<<" "<<cornercovdistance<<" "<<result<<endl<<flush;
323 
324  }
325 
326  }
327  vrml.VRMLOut("scene.wrl");
328  return 0;
329 }
330 
analytic properties of affine image warp
class HomgPoint2D describes a point with 2 degrees of freedom in projective coordinates.
Definition: HomgPoint2D.hh:67
void AddDebugLevel(const long int lv)
Definition: Debug.hh:355
int VRMLOut(const std::string &sFilename)
flush all 3d objects to a vrml file with name sFilename, this is the function most users would call ...
Definition: ThreeDOut.cpp:3670
void SetOrigin(const Vector2< double > &origin)
origin relative to which rotation and scale is performed
a 3x3 Matrix describing projective transformations between planes
Definition: HMatrix.hh:39
Maps image src to image sink with homography H (software implementation)
affine transformation of 2D image plane which relates image coordinate system and local affine featur...
Unified output of 3D entities via OpenGL or VRML.
Definition: ThreeDOut.hh:349
Matrix< T > & newsize(Subscript M, Subscript N)
Definition: cmat.h:269
MetaData * GetMetaData()
Definition: ImageBase.hh:436
float image storage type
Definition: ImageBase.hh:118
int AutoAlign(const PyramidImage< float > &I1, const PyramidImage< float > &I2, Vector< double > &params, Matrix< double > &Cov)
aligns image1 with image2 with respect to parameters params, scale space tracking is used recursively...
unsigned int Size() const
length of the vector
Definition: Vector.hh:134
void SetA(const BIAS::Matrix2x2< double > &A, const BIAS::Matrix< double > &cov=BIAS::Matrix< double >(4, 4, BIAS::MatrixZero))
return affine matrix
unsigned int GetWidth() const
Definition: ImageBase.hh:292
3D rotation matrix
Definition: RMatrix.hh:49
void SetAlignmentPixelsRectangular(const LocalAffineFrame &LAF, int numberOfPixels=-1)
choose any rectangular tracking region and resolution
Vector< T > & newsize(Subscript N)
Definition: vec.h:220
void SetUpdateThreshold(const double maxerr)
if a step smaller (norm L2) than this is taken, convergence is declared and the system finishes ...
void Set(const Vector3< ROTATION_MATRIX_TYPE > &w, const ROTATION_MATRIX_TYPE phi)
Set from rotation axis w and angle phi (in rad)
Vector< double > ExtractParameters(const Matrix3x3< double > &A) const
static int ConvertST(const BIAS::ImageBase &source, BIAS::ImageBase &dest, ImageBase::EStorageType targetST)
Function to convert the storage type of images e.g.
int GetInverse(Matrix3x3< T > &inv) const
Matrix inversion: inverts this and stores resulty in argument inv.
Definition: Matrix3x3.cpp:373
unsigned int GetChannelCount() const
returns the number of Color channels, e.g.
Definition: ImageBase.hh:362
void SetT(const BIAS::Vector2< double > &T, const BIAS::Matrix< double > &cov=BIAS::Matrix< double >(2, 2, BIAS::MatrixZero))
return offset
void Init(const Image< StorageType > &image, const unsigned py_size=0)
copy image into level 0 and create other levels according to parameters set so far (pyramidsize...
void MapAcrossPlane(const PMatrix &P1, const PMatrix &P2, const HomgPlane3D &scenePlane=HomgPlane3D(0, 0, 0, 1))
Set (*this) to map coordinates from camera p1 to p2 using the given scene plane.
Definition: HMatrix.cpp:81
HMatrix GetLinearized(const HomgPoint2D &x) const
returns 1st order Taylor expansion of homography at x
Definition: HMatrix.cpp:55
unsigned int GetHeight() const
Definition: ImageBase.hh:299
int Map(const Image< InputStorageType > &src, Image< OutputStorageType > &sink, InterpolationMethod=MapTrilinear, bool newSink=false, double SuperSampling=1.0)
backward mapping with various interpolations
Inverse Compositional Image Alignment (&quot;Registration&quot;)
BIAS::Vector4< unsigned char > RGBAuc
Definition: RGBA.hh:47
static int Save(const std::string &filename, const ImageBase &img, const enum TFileFormat FileFormat=FF_auto, const bool sync=BIAS_DEFAULT_SYNC, const int c_jpeg_quality=BIAS_DEFAULT_IMAGE_QUALITY, const bool forceNewID=BIAS_DEFAULT_FORCENEWID, const bool &writeMetaData=true)
Export image as file using extrnal libs.
Definition: ImageIO.cpp:725
void SetIdentity()
Converts matrix to identity matrix.
Definition: Matrix.cpp:383
static std::string toString(const T thenumber, int numzeros=DEFAULT_LEADING_ZEROS)
Converts number to string with leading zeros.
Definition: FileHandling.hh:98
const StorageType * GetImageData() const
overloaded GetImageData() from ImageBase
Definition: Image.hh:136
static int Load(const std::string &FileName, ImageBase &img)
first tries a call to Read MIP image and if that fails, tries to Import Image with all other availabl...
Definition: ImageIO.cpp:141
K describes the mapping from world coordinates (wcs) to pixel coordinates (pcs).
Definition: KMatrix.hh:48
void SetTextureTransform(const TextureTransform &T)
pass object type which holds texture transform
double ComputeLinearizationError(const HomgPoint2D &x0, const HomgPoint2D &y) const
compute distance between y mapped with homography or first order taylor approximation at development ...
Definition: HMatrix.cpp:332
unsigned int AddPMatrix(BIAS::PMatrix &P, const unsigned int &width, const unsigned int &height, const BIAS::RGBAuc &Color=RGBAuc_WHITE_OPAQUE, const double &dScale=DEF_P_SCALE, const std::string &name="")
decompose P and add to internal data structures.
Definition: ThreeDOut.cpp:1029
double MahalanobisDistance(const BIAS::Matrix< double > &Sigma, const BIAS::Vector< double > &V)
computes squared Mahalanobis distance V^T Sigma^-1 V efficiently using cholesky decomposition and exp...
Definition: Lapack.cpp:780
double GetNormalDistributed(const double mean, const double sigma)
on succesive calls return normal distributed random variable with mean and standard deviation sigma ...
Definition: Random.hh:71
describes a projective 3D -&gt; 2D mapping in homogenous coordinates
Definition: PMatrix.hh:88
unsigned int AddPoint(const BIAS::Vector3< double > &v, const BIAS::RGBAuc &Color=RGBAuc_WHITE_OPAQUE)
Definition: ThreeDOut.cpp:719
void SetDampening(const double &dampening)
adds a dampening d to hessian&#39;s diag
void SetMetaData(const MetaData &m)
Definition: ImageBase.hh:450
unsigned long int GetPixelCount() const
returns number of pixels in image
Definition: ImageBase.hh:402
This is the base class for images in BIAS.
Definition: ImageBase.hh:102
void SetHomography(const HMatrix &H)
set your homography H (source = H * sink) before calling Map()
class for producing random numbers from different distributions
Definition: Random.hh:51
void SetZero()
zeroes the image
Definition: ImageBase.hh:83
class BIASGeometryBase_EXPORT HomgPoint2D
void SetMaxIterations(const int maxiter)
number of iterations, before the system stops unsuccessfully
void SetIdentity()
set the elements of this matrix to the identity matrix (possibly overriding the inherited method) ...
Definition: Matrix3x3.hh:429
void SetParameters(const Vector< double > &p)
a11-1, a12, a21, a22-1, dx, dy
A homogeneous plane (in P^3) All points X on the plane p fulfill p &#39; * X = 0.
Definition: HomgPlane3D.hh:46
static int ToGrey(const ImageBase &source, ImageBase &dest)
wrapper for the templated function ToGrey