Go to the documentation of this file.00001 #include "vgl_h_matrix_1d_compute_optimize.h"
00002 #include <vgl/algo/vgl_h_matrix_1d_compute_linear.h>
00003
00004 #include <vcl_cassert.h>
00005 #include <vcl_iostream.h>
00006 #include <vcl_iomanip.h>
00007
00008 #include <vnl/vnl_least_squares_function.h>
00009 #include <vnl/algo/vnl_levenberg_marquardt.h>
00010 #include <vnl/vnl_double_2x2.h>
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00020 class XXX : public vnl_least_squares_function
00021 {
00022 private:
00023 unsigned N;
00024 const vcl_vector<double> &z1,z2;
00025 public:
00026 XXX(const vcl_vector<double> &z1_,const vcl_vector<double> &z2_)
00027 : vnl_least_squares_function(3, z1_.size(), no_gradient)
00028 , N(z1_.size())
00029 , z1(z1_) , z2(z2_)
00030 {
00031 #ifdef DEBUG
00032 vcl_cerr << "N=" << N << '\n';
00033 #endif
00034 assert(N == z1.size());
00035 assert(N == z2.size());
00036 }
00037 ~XXX() { N=0; }
00038
00039 void boo(const vnl_vector<double> &x) {
00040 assert(x.size()==3);
00041 vcl_cerr << vcl_showpos << vcl_fixed;
00042 double z,y;
00043 for (unsigned i=0;i<N;i++) {
00044 z=z1[i];
00045 y=(z+x[0])/(x[1]*z+1+x[2]);
00046 vcl_cerr << z << ' ' << y << '[' << z2[i] << ']' << vcl_endl;
00047 }
00048 }
00049
00050
00051
00052 void f(const vnl_vector<double>& x, vnl_vector<double>& fx) {
00053 assert(x.size()==3);
00054 assert(fx.size()==N);
00055 double z,y;
00056 for (unsigned k=0;k<N;k++) {
00057 z=z1[k];
00058 y=(z+x[0])/(x[1]*z+1+x[2]);
00059 fx[k]=z2[k]-y;
00060 }
00061 }
00062 };
00063 #endif // DOXYGEN_SHOULD_SKIP_THIS
00064
00065 static
00066 void do_compute(const vcl_vector<double> &z1,const vcl_vector<double> &z2,vgl_h_matrix_1d<double> &M)
00067 {
00068
00069
00070
00071
00072
00073
00074 XXX f(z1,z2);
00075
00076
00077 vnl_vector<double> x(3);
00078 x.fill(0);
00079
00080
00081 vnl_levenberg_marquardt LM(f);
00082
00083
00084
00085 LM.minimize(x);
00086
00087
00088
00089 vnl_double_2x2 T;
00090 T(0,0)=1; T(0,1)=x[0];
00091 T(1,0)=x[1]; T(1,1)=1+x[2];
00092 M.set(T);
00093 }
00094
00095 bool vgl_h_matrix_1d_compute_optimize::
00096 compute_cool_homg(const vcl_vector<vgl_homg_point_1d<double> >&p1,
00097 const vcl_vector<vgl_homg_point_1d<double> >&p2,
00098 vgl_h_matrix_1d<double>& M)
00099 {
00100 unsigned N=p1.size();
00101 assert(N==p2.size());
00102 if (N<3) return false;
00103
00104 vcl_vector<double> z1(N,0.0),z2(N,0.0);
00105 vgl_h_matrix_1d_compute_linear C;
00106 C.compute(p1,p2,M);
00107
00108
00109
00110 for (unsigned i=0;i<N;i++) {
00111 vgl_homg_point_1d<double> v = M(p1[i]);
00112 if (v.w() == 0.0) return false;
00113 z1[i] = v.x()/v.w();
00114 if (p2[i].w()) return false;
00115 z2[i] = p2[i].x()/p2[i].w();
00116 }
00117
00118 vgl_h_matrix_1d<double> K;
00119 do_compute(z1,z2,K);
00120 M= K * M;
00121 return true;
00122 }
00123