contrib/gel/vmal/vmal_homog2d.cxx
Go to the documentation of this file.
00001 // This is gel/vmal/vmal_homog2d.cxx
00002 #include "vmal_homog2d.h"
00003 #include <vnl/algo/vnl_svd.h>
00004 
00005 vmal_homog2d::vmal_homog2d()
00006 {
00007 }
00008 
00009 vmal_homog2d::~vmal_homog2d()
00010 {
00011 }
00012 
00013 void vmal_homog2d::compute_homo(const vcl_vector<vnl_double_3 > &pima1,
00014                                 const vcl_vector<vnl_double_3 > &pima2,
00015                                 vnl_double_3x3 &hmatrix)
00016 {
00017   int numpoints=pima1.size();
00018   vnl_matrix<double> A(2*numpoints,9);
00019   A.fill(0.0);
00020   for (int row=0,i=0; i<numpoints; ++i)
00021   {
00022     A(row, 0) = pima1[i][0] * pima2[i][2];
00023     A(row, 1) = pima1[i][1] * pima2[i][2];
00024     A(row, 2) = pima1[i][2] * pima2[i][2];
00025     A(row, 3) = 0;
00026     A(row, 4) = 0;
00027     A(row, 5) = 0;
00028     A(row, 6) = -pima1[i][0] * pima2[i][0];
00029     A(row, 7) = -pima1[i][1] * pima2[i][0];
00030     A(row, 8) = -pima1[i][2] * pima2[i][0];
00031     ++row;
00032 
00033     A(row, 0) = 0;
00034     A(row, 1) = 0;
00035     A(row, 2) = 0;
00036     A(row, 3) = pima1[i][0] * pima2[i][2];
00037     A(row, 4) = pima1[i][1] * pima2[i][2];
00038     A(row, 5) = pima1[i][2] * pima2[i][2];
00039     A(row, 6) = -pima1[i][0] * pima2[i][1];
00040     A(row, 7) = -pima1[i][1] * pima2[i][1];
00041     A(row, 8) = -pima1[i][2] * pima2[i][1];
00042     ++row;
00043   }
00044 
00045   A.normalize_rows();
00046   vnl_svd<double> SVD(A);
00047   vnl_vector<double> x=SVD.nullvector();
00048     hmatrix[0][0]=x[0]; hmatrix[0][1]=x[1]; hmatrix[0][2]=x[2];
00049   hmatrix[1][0]=x[3]; hmatrix[1][1]=x[4]; hmatrix[1][2]=x[5];
00050   hmatrix[2][0]=x[6]; hmatrix[2][1]=x[7]; hmatrix[2][2]=x[8];
00051 }