00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #ifndef __IRLS_H
00032 #define __IRLS_H
00033
00034 #include <TooN/wls.h>
00035 #include <cassert>
00036 #include <cmath>
00037
00038 namespace TooN {
00039
00044 template<typename Precision>
00045 struct RobustI {
00046 void set_sd(Precision x){ sd_inlier = x;}
00047 double sd_inlier;
00048 inline Precision reweight(Precision x) {return 1/(sd_inlier+fabs(x));}
00049 inline Precision true_scale(Precision x) {return reweight(x) - fabs(x)*reweight(x)*reweight(x);}
00050 inline Precision objective(Precision x) {return fabs(x) + sd_inlier*log(sd_inlier*reweight(x));}
00051 };
00052
00057 template<typename Precision>
00058 struct RobustII {
00059 void set_sd(Precision x){ sd_inlier = x*x;}
00060 Precision sd_inlier;
00061 inline Precision reweight(Precision d){return 1/(sd_inlier+d*d);}
00062 inline Precision true_scale(Precision d){return d - 2*d*reweight(d);}
00063 inline Precision objective(Precision d){return 0.5 * log(1 + d*d/sd_inlier);}
00064 };
00065
00070 template<typename Precision>
00071 struct ILinear {
00072 void set_sd(Precision){}
00073 inline Precision reweight(Precision d){return 1;}
00074 inline Precision true_scale(Precision d){return 1;}
00075 inline Precision objective(Precision d){return d*d;}
00076 };
00077
00083 template<typename Precision>
00084 struct RobustIII {
00085
00086 void set_sd(Precision x){ sd_inlier = x*x;}
00087 Precision sd_inlier;
00088
00089 Precision reweight(Precision x) const
00090 {
00091 double d = (1 + x*x/sd_inlier);
00092 return 1/(d*d);
00093 }
00095 Precision objective(Precision x) const
00096 {
00097 return x*x / (2*(1 + x*x/sd_inlier));
00098 }
00099 };
00100
00106 template <int Size, typename Precision, template <typename Precision> class Reweight>
00107 class IRLS
00108 : public Reweight<Precision>,
00109 public WLS<Size,Precision>
00110 {
00111 public:
00112 IRLS(int size=Size):
00113 WLS<Size,Precision>(size),
00114 my_true_C_inv(Zeros(size))
00115 {
00116 my_residual=0;
00117 }
00118
00119 template<int Size2, typename Precision2, typename Base2>
00120 inline void add_mJ(Precision m, const Vector<Size2,Precision2,Base2>& J) {
00121 SizeMismatch<Size,Size2>::test(my_true_C_inv.num_rows(), J.size());
00122
00123 Precision scale = Reweight<Precision>::reweight(m);
00124 Precision ts = Reweight<Precision>::true_scale(m);
00125 my_residual += Reweight<Precision>::objective(m);
00126
00127 WLS<Size>::add_mJ(m,J,scale);
00128
00129 Vector<Size,Precision> scaledm(m*ts);
00130
00131 my_true_C_inv += scaledm.as_col() * scaledm.as_row();
00132
00133 }
00134
00135 void operator += (const IRLS& meas){
00136 WLS<Size>::operator+=(meas);
00137 my_true_C_inv += meas.my_true_C_inv;
00138 }
00139
00140
00141 Matrix<Size,Size,Precision>& get_true_C_inv() {return my_true_C_inv;}
00142 const Matrix<Size,Size,Precision>& get_true_C_inv()const {return my_true_C_inv;}
00143
00144 Precision get_residual() {return my_residual;}
00145
00146 void clear(){
00147 WLS<Size,Precision>::clear();
00148 my_residual=0;
00149 my_true_C_inv = Zeros;
00150 }
00151
00152 private:
00153
00154 Precision my_residual;
00155
00156 Matrix<Size,Size,Precision> my_true_C_inv;
00157
00158
00159 IRLS( IRLS& copyof );
00160 int operator = ( IRLS& copyof );
00161 };
00162
00163 }
00164
00165 #endif