00001
00002 #include "vgl_h_matrix_2d_compute_rigid_body.h"
00003
00004
00005
00006 #include <vcl_iostream.h>
00007 #include <vcl_cmath.h>
00008 #include <vcl_cassert.h>
00009 #include <vnl/vnl_inverse.h>
00010 #include <vnl/vnl_transpose.h>
00011 #include <vnl/algo/vnl_svd.h>
00012 #include <vgl/algo/vgl_norm_trans_2d.h>
00013
00014
00015 vgl_h_matrix_2d_compute_rigid_body::vgl_h_matrix_2d_compute_rigid_body()
00016 {
00017 }
00018
00019 const int TM_UNKNOWNS_COUNT = 3;
00020 const double DEGENERACY_THRESHOLD = 0.01;
00021
00022
00023
00024
00025
00026 bool vgl_h_matrix_2d_compute_rigid_body::
00027 solve_rigid_body_problem(int equ_count,
00028 vcl_vector<vgl_homg_point_2d<double> > const& p1,
00029 vcl_vector<vgl_homg_point_2d<double> > const& p2,
00030 vgl_h_matrix_2d<double>& H)
00031 {
00032
00033 vnl_matrix<double> D(equ_count, TM_UNKNOWNS_COUNT+2);
00034 int n = p1.size();
00035 int row = 0;
00036 for (int i = 0; i < n; i++) {
00037 D(row, 0) = p1[i].x();
00038 D(row, 1) = -p1[i].y();
00039 D(row, 2) = 1.0;
00040 D(row, 3) = 0.0;
00041 D(row, 4) = -p2[i].x();
00042 ++row;
00043
00044 D(row, 0) = p1[i].y();
00045 D(row, 1) = p1[i].x();
00046 D(row, 2) = 0;
00047 D(row, 3) = 1.0;
00048 D(row, 4) = -p2[i].y();
00049 ++row;
00050 }
00051
00052 D.normalize_rows();
00053 vnl_svd<double> svd(D);
00054 vnl_vector<double> nullv = svd.nullvector();
00055
00056 if (svd.W(4)>DEGENERACY_THRESHOLD*svd.W(3)) {
00057 vcl_cout << "vgl_h_matrix_2d_compute_rigid_body : inaccurate solution probably due to inconsistent point correspondences\n"
00058 << "W\n" << svd.W() << vcl_endl;
00059 return false;
00060 }
00061
00062
00063
00064 double norm = nullv[4];
00065 nullv /= norm;
00066
00067 double y = nullv[1];
00068 double x = nullv[0];
00069 double angle = vcl_atan2(y,x);
00070 double c = vcl_cos(angle);
00071 double s = vcl_sin(angle);
00072 vnl_matrix_fixed<double, 3,3> M;
00073 M[0][0]=c; M[0][1]= -s; M[0][2] = nullv[2];
00074 M[1][0]=s; M[1][1]= c; M[1][2] = nullv[3];
00075 M[2][0]=0; M[2][1]= 0; M[2][2] = 1;
00076 H.set(M);
00077 return true;
00078 }
00079
00080 bool vgl_h_matrix_2d_compute_rigid_body::
00081 compute_p(vcl_vector<vgl_homg_point_2d<double> > const& points1,
00082 vcl_vector<vgl_homg_point_2d<double> > const& points2,
00083 vgl_h_matrix_2d<double>& H)
00084 {
00085
00086 assert(points1.size() == points2.size());
00087 int n = points1.size();
00088
00089 int equ_count = n * (2);
00090 if (n * 2 < TM_UNKNOWNS_COUNT) {
00091 vcl_cerr << "vgl_h_matrix_2d_compute_rigid_body: Need at least 2 matches.\n";
00092 if (n == 0) vcl_cerr << "Could be vcl_vector setlength idiosyncrasies!\n";
00093 return false;
00094 }
00095
00096 vgl_norm_trans_2d<double> tr1, tr2;
00097 if (!tr1.compute_from_points(points1))
00098 return false;
00099 if (!tr2.compute_from_points(points2))
00100 return false;
00101 vcl_vector<vgl_homg_point_2d<double> > tpoints1, tpoints2;
00102 for (int i = 0; i<n; i++)
00103 {
00104 tpoints1.push_back(tr1(points1[i]));
00105 tpoints2.push_back(tr2(points2[i]));
00106 }
00107 vgl_h_matrix_2d<double> hh;
00108 if (!solve_rigid_body_problem(equ_count, tpoints1, tpoints2, hh))
00109 return false;
00110
00111
00112
00113
00114
00115
00116
00117
00118 vgl_h_matrix_2d<double> tr2_inv = tr2.get_inverse();
00119 H = tr2_inv*hh*tr1;
00120 return true;
00121 }
00122
00123 bool vgl_h_matrix_2d_compute_rigid_body::
00124 compute_l(vcl_vector<vgl_homg_line_2d<double> > const& lines1,
00125 vcl_vector<vgl_homg_line_2d<double> > const& lines2,
00126 vgl_h_matrix_2d<double>& H)
00127 {
00128
00129 assert(lines1.size() == lines2.size());
00130 int n = lines1.size();
00131 int equ_count = 2*n;
00132
00133
00134 vgl_norm_trans_2d<double> tr1, tr2;
00135 if (!tr1.compute_from_lines(lines1))
00136 return false;
00137 if (!tr2.compute_from_lines(lines2))
00138 return false;
00139 vcl_vector<vgl_homg_point_2d<double> > tlines1, tlines2;
00140 for (vcl_vector<vgl_homg_line_2d<double> >::const_iterator
00141 lit = lines1.begin(); lit != lines1.end(); lit++)
00142 {
00143
00144 vgl_homg_line_2d<double> l = tr1(*lit);
00145
00146 vgl_homg_point_2d<double> p(l.a(), l.b(), l.c());
00147 tlines1.push_back(p);
00148 }
00149 for (vcl_vector<vgl_homg_line_2d<double> >::const_iterator
00150 lit = lines2.begin(); lit != lines2.end(); lit++)
00151 {
00152
00153 vgl_homg_line_2d<double> l = tr2(*lit);
00154
00155 vgl_homg_point_2d<double> p(l.a(), l.b(), l.c());
00156 tlines2.push_back(p);
00157 }
00158
00159 vgl_h_matrix_2d<double> hl,hp,tr2inv;
00160 if (!solve_rigid_body_problem(equ_count, tlines1, tlines2, hl))
00161 return false;
00162
00163
00164 vnl_matrix_fixed<double, 3, 3> const & Ml = hl.get_matrix();
00165 vnl_matrix_fixed<double, 3, 3> Mp = vnl_inverse_transpose(Ml);
00166 hp.set(Mp);
00167
00168
00169
00170
00171
00172
00173
00174
00175 tr2inv = tr2.get_inverse();
00176 H = tr2inv*hp*tr1;
00177 return true;
00178 }
00179
00180 bool vgl_h_matrix_2d_compute_rigid_body::
00181 compute_pl(vcl_vector<vgl_homg_point_2d<double> > const& points1,
00182 vcl_vector<vgl_homg_point_2d<double> > const& points2,
00183 vcl_vector<vgl_homg_line_2d<double> > const& lines1,
00184 vcl_vector<vgl_homg_line_2d<double> > const& lines2,
00185 vgl_h_matrix_2d<double>& H)
00186 {
00187
00188 assert(points1.size() == points2.size());
00189 int np = points1.size();
00190
00191 assert(lines1.size() == lines2.size());
00192 int nl = lines1.size();
00193
00194 int equ_count = np * 2 + 2*nl;
00195 if ((np+nl)*2 < TM_UNKNOWNS_COUNT)
00196 {
00197 vcl_cerr << "vgl_h_matrix_2d_compute_rigid_body: Need at least 4 matches.\n";
00198 if (np+nl == 0) vcl_cerr << "Could be vcl_vector setlength idiosyncrasies!\n";
00199 return false;
00200 }
00201
00202 vgl_norm_trans_2d<double> tr1, tr2;
00203 if (!tr1.compute_from_points_and_lines(points1,lines1))
00204 return false;
00205 if (!tr2.compute_from_points_and_lines(points2,lines2))
00206 return false;
00207 vcl_vector<vgl_homg_point_2d<double> > tpoints1, tpoints2;
00208 for (int i = 0; i<np; i++)
00209 {
00210 tpoints1.push_back(tr1(points1[i]));
00211 tpoints2.push_back(tr2(points2[i]));
00212 }
00213 for (int i = 0; i<nl; i++)
00214 {
00215 double a=lines1[i].a(), b=lines1[i].b(), c=lines1[i].c(), d=vcl_sqrt(a*a+b*b);
00216 tpoints1.push_back(tr1(vgl_homg_point_2d<double>(-a*c,-b*c,d)));
00217 a=lines2[i].a(), b=lines2[i].b(), c=lines2[i].c(), d = vcl_sqrt(a*a+b*b);
00218 tpoints2.push_back(tr2(vgl_homg_point_2d<double>(-a*c,-b*c,d)));
00219 }
00220 vgl_h_matrix_2d<double> hh;
00221 if (!solve_rigid_body_problem(equ_count, tpoints1, tpoints2, hh))
00222 return false;
00223
00224 vgl_h_matrix_2d<double> tr2_inv = tr2.get_inverse();
00225 H = tr2_inv*hh*tr1;
00226 return true;
00227 }
00228
00229 bool vgl_h_matrix_2d_compute_rigid_body::
00230 compute_l(vcl_vector<vgl_homg_line_2d<double> > const& ,
00231 vcl_vector<vgl_homg_line_2d<double> > const& ,
00232 vcl_vector<double> const& ,
00233 vgl_h_matrix_2d<double>& )
00234 {
00235 return false;
00236 }
00237
00238 #if 0 // do later
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256 bool vgl_h_matrix_2d_compute_rigid_body::
00257 solve_weighted_least_squares(vcl_vector<vgl_homg_line_2d<double> > const& l1,
00258 vcl_vector<vgl_homg_line_2d<double> > const& l2,
00259 vcl_vector<double> const& w,
00260 vgl_h_matrix_2d<double>& H)
00261 {
00262 int Nc = l1.size();
00263
00264
00265 vnl_vector<double> two_w(2*Nc);
00266 int j =0;
00267 for (int i = 0; i<Nc; i++, j+=2)
00268 {
00269 two_w[j]=w[i];
00270 two_w[j+1]=w[i];
00271 }
00272 vnl_diag_matrix<double> W(two_w);
00273
00274
00275 vnl_matrix<double> D(2*Nc, TM_UNKNOWNS_COUNT);
00276 vnl_matrix<double> M(TM_UNKNOWNS_COUNT, TM_UNKNOWNS_COUNT);
00277
00278 int row = 0;
00279 for (int i = 0; i < Nc; i++)
00280 {
00281 D(row, 0) = l1[i].a() * l2[i].c();
00282 D(row, 1) = l1[i].b() * l2[i].c();
00283 D(row, 2) = l1[i].c() * l2[i].c();
00284 D(row, 3) = 0;
00285 D(row, 4) = 0;
00286 D(row, 5) = 0;
00287 D(row, 6) = -l1[i].a() * l2[i].a();
00288 D(row, 7) = -l1[i].b() * l2[i].a();
00289 D(row, 8) = -l1[i].c() * l2[i].a();
00290 ++row;
00291
00292 D(row, 0) = 0;
00293 D(row, 1) = 0;
00294 D(row, 2) = 0;
00295 D(row, 3) = l1[i].a() * l2[i].c();
00296 D(row, 4) = l1[i].b() * l2[i].c();
00297 D(row, 5) = l1[i].c() * l2[i].c();
00298 D(row, 6) = -l1[i].a() * l2[i].b();
00299 D(row, 7) = -l1[i].b() * l2[i].b();
00300 D(row, 8) = -l1[i].c() * l2[i].b();
00301 ++row;
00302 }
00303 M = vnl_transpose(D)*W*D;
00304 D.normalize_rows();
00305 vnl_svd<double> svd(D);
00306
00307
00308
00309
00310 if (svd.W(7)<DEGENERACY_THRESHOLD*svd.W(8)) {
00311 vcl_cerr << "vgl_h_matrix_2d_compute_rigid_body : design matrix has rank < 8\n"
00312 << "vgl_h_matrix_2d_compute_rigid_body : probably due to degenerate point configuration\n";
00313 return false;
00314 }
00315
00316 H.set(svd.nullvector().data_block());
00317 return true;
00318 }
00319
00320 bool vgl_h_matrix_2d_compute_rigid_body::
00321 compute_l(vcl_vector<vgl_homg_line_2d<double> > const& lines1,
00322 vcl_vector<vgl_homg_line_2d<double> > const& lines2,
00323 vcl_vector<double> const & weights,
00324 vgl_h_matrix_2d<double>& H)
00325 {
00326
00327 assert(lines1.size() == lines2.size());
00328
00329
00330
00331 vgl_norm_trans_2d<double> tr1, tr2;
00332 if (!tr1.compute_from_lines(lines1))
00333 return false;
00334 if (!tr2.compute_from_lines(lines2))
00335 return false;
00336 vcl_vector<vgl_homg_line_2d<double> > tlines1, tlines2;
00337 for (vcl_vector<vgl_homg_line_2d<double> >::const_iterator
00338 lit = lines1.begin(); lit != lines1.end(); lit++)
00339 {
00340
00341 vgl_homg_line_2d<double> l = tr1(*lit);
00342 tlines1.push_back(l);
00343 }
00344 for (vcl_vector<vgl_homg_line_2d<double> >::const_iterator
00345 lit = lines2.begin(); lit != lines2.end(); lit++)
00346 {
00347
00348 vgl_homg_line_2d<double> l = tr2(*lit);
00349 tlines2.push_back(l);
00350 }
00351
00352 vgl_h_matrix_2d<double> hl,hp,tr2inv;
00353 if (!solve_weighted_least_squares(tlines1, tlines2, weights, hl))
00354 return false;
00355
00356
00357 vnl_matrix_fixed<double, 3, 3> const & Ml = hl.get_matrix();
00358 vnl_matrix_fixed<double, 3, 3> Mp = vnl_inverse_transpose(Ml);
00359 hp.set(Mp);
00360
00361
00362
00363
00364
00365
00366
00367
00368 tr2inv = tr2.get_inverse();
00369 H = tr2inv*hp*tr1;
00370 return true;
00371 }
00372 #endif // 0