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 }