00001
00002 #ifndef vpdt_norm_metric_h_
00003 #define vpdt_norm_metric_h_
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include <vpdl/vpdt/vpdt_field_traits.h>
00017 #include <vpdl/vpdt/vpdt_access.h>
00018 #include <vpdl/vpdt/vpdt_eigen_sym_matrix.h>
00019 #include <vcl_limits.h>
00020 #include <vcl_cmath.h>
00021 #include <vcl_cassert.h>
00022
00023
00024
00025
00026
00027 template<class F, class Tensor, class Disambiguate= void>
00028 struct vpdt_norm_metric;
00029
00030
00031
00032 template<class F>
00033 struct vpdt_norm_metric<F, typename vpdt_eigen_sym_matrix_gen<F>::type,
00034 typename vpdt_field_traits<F>::type_is_vector >
00035 {
00036
00037 typedef typename vpdt_eigen_sym_matrix_gen<F>::type covar_type;
00038
00039 typedef typename vpdt_field_traits<F>::scalar_type T;
00040
00041 typedef typename vpdt_field_traits<F>::vector_type vector;
00042
00043 typedef typename vpdt_field_traits<F>::matrix_type matrix;
00044
00045
00046 static inline T distance(const F& pt1, const F& pt2, const covar_type& c)
00047 {
00048 return vcl_sqrt(sqr_distance(pt1,pt2,c));
00049 }
00050
00051
00052 static inline T sqr_distance(const F& pt1, const F& pt2, const covar_type& c)
00053 {
00054
00055 vector d(pt1-pt2);
00056 return c.inverse_quad_form(d);
00057 }
00058
00059
00060 static inline T sqr_distance_deriv(const F& pt1, const F& pt2,
00061 const covar_type& c, vector& g)
00062 {
00063
00064 vector d(pt1-pt2);
00065 c.inverse_product(d,g);
00066 T sqr_dist = dot_product(d,g);
00067 g *= 2;
00068 return sqr_dist;
00069 }
00070
00071
00072
00073 static inline void compute_covar(matrix& covar, const F& , const covar_type& c)
00074 {
00075 c.form_matrix(covar);
00076 }
00077
00078
00079
00080 static inline T covar_det(const F& , const covar_type& c)
00081 {
00082 return c.determinant();
00083 }
00084 };
00085
00086
00087
00088
00089 template<class F>
00090 struct vpdt_norm_metric<F, typename vpdt_field_traits<F>::vector_type,
00091 typename vpdt_field_traits<F>::type_is_vector >
00092 {
00093
00094 typedef typename vpdt_field_traits<F>::vector_type covar_type;
00095
00096 typedef typename vpdt_field_traits<F>::scalar_type T;
00097
00098 typedef typename vpdt_field_traits<F>::vector_type vector;
00099
00100 typedef typename vpdt_field_traits<F>::matrix_type matrix;
00101
00102
00103 static inline T distance(const F& pt1, const F& pt2, const covar_type& c)
00104 {
00105 return vcl_sqrt(sqr_distance(pt1,pt2,c));
00106 }
00107
00108
00109 static inline T sqr_distance(const F& pt1, const F& pt2, const covar_type& c)
00110 {
00111 const unsigned int d = vpdt_size(c);
00112 assert(vpdt_size(pt1) == d);
00113 assert(vpdt_size(pt2) == d);
00114 T val = T(0);
00115 for (unsigned int i=0; i<d; ++i)
00116 {
00117 T tmp = (vpdt_index(pt1,i)-vpdt_index(pt2,i));
00118 val += tmp*tmp/vpdt_index(c,i);
00119 }
00120 return val;
00121 }
00122
00123
00124 static inline T sqr_distance_deriv(const F& pt1, const F& pt2,
00125 const covar_type& c, vector& g)
00126 {
00127 const unsigned int d = vpdt_size(c);
00128 vpdt_set_size(g,d);
00129 assert(vpdt_size(pt1) == d);
00130 assert(vpdt_size(pt2) == d);
00131 T val = T(0);
00132 for (unsigned int i=0; i<d; ++i)
00133 {
00134 T tmp = (vpdt_index(pt1,i)-vpdt_index(pt2,i));
00135 T& g_i = vpdt_index(g,i) = tmp/vpdt_index(c,i);
00136 val += tmp*g_i;
00137 g_i *= 2;
00138 }
00139 return val;
00140 }
00141
00142
00143
00144 static inline void compute_covar(matrix& covar, const F& , const covar_type& c)
00145 {
00146 const unsigned int d = c.size();
00147 vpdt_set_size(covar,d);
00148 assert(vpdt_size(covar) == d);
00149 for (unsigned int i=0; i<d; ++i)
00150 {
00151 vpdt_index(covar,i,i) = vpdt_index(c,i);
00152 for (unsigned int j=i+1; j<d; ++j)
00153 vpdt_index(covar,i,j) = vpdt_index(covar,j,i) = T(0);
00154 }
00155 }
00156
00157
00158
00159 static inline T covar_det(const F& pt, const covar_type& c)
00160 {
00161 const unsigned int d = c.size();
00162 double det = 1.0;
00163 for (unsigned int i=0; i<d; ++i)
00164 det *= vpdt_index(c,i);
00165 return static_cast<T>(det);
00166 }
00167 };
00168
00169
00170
00171
00172 template<class F>
00173 struct vpdt_norm_metric<F, typename vpdt_field_traits<F>::scalar_type,
00174 typename vpdt_field_traits<F>::type_is_vector >
00175 {
00176
00177 typedef typename vpdt_field_traits<F>::scalar_type covar_type;
00178
00179 typedef typename vpdt_field_traits<F>::scalar_type T;
00180
00181 typedef typename vpdt_field_traits<F>::vector_type vector;
00182
00183 typedef typename vpdt_field_traits<F>::matrix_type matrix;
00184
00185
00186 static inline T distance(const F& pt1, const F& pt2, const covar_type& c)
00187 {
00188 return vcl_sqrt(sqr_distance(pt1,pt2,c));
00189 }
00190
00191
00192 static inline T sqr_distance(const F& pt1, const F& pt2, const covar_type& c)
00193 {
00194 const unsigned int d = vpdt_size(pt1);
00195 assert(vpdt_size(pt2) == d);
00196 T val = T(0);
00197 for (unsigned int i=0; i<d; ++i)
00198 {
00199 T tmp = (vpdt_index(pt1,i)-vpdt_index(pt2,i));
00200 val += tmp*tmp/c;
00201 }
00202 return val;
00203 }
00204
00205
00206 static inline T sqr_distance_deriv(const F& pt1, const F& pt2,
00207 const covar_type& c, vector& g)
00208 {
00209 const unsigned int d = vpdt_size(pt1);
00210 vpdt_set_size(g,d);
00211 assert(vpdt_size(pt2) == d);
00212 T val = T(0);
00213 for (unsigned int i=0; i<d; ++i)
00214 {
00215 T tmp = (vpdt_index(pt1,i)-vpdt_index(pt2,i));
00216 T& g_i = vpdt_index(g,i) = tmp/c;
00217 val += tmp*g_i;
00218 g_i *= 2;
00219 }
00220 return val;
00221 }
00222
00223
00224
00225 static inline void compute_covar(matrix& covar, const F& pt, const covar_type& c)
00226 {
00227 const unsigned int d = vpdt_size(pt);
00228 vpdt_set_size(covar,d);
00229 assert(vpdt_size(covar) == d);
00230 for (unsigned int i=0; i<d; ++i)
00231 {
00232 vpdt_index(covar,i,i) = c;
00233 for (unsigned int j=i+1; j<d; ++j)
00234 vpdt_index(covar,i,j) = vpdt_index(covar,j,i) = T(0);
00235 }
00236 }
00237
00238
00239
00240 static inline T covar_det(const F& pt, const covar_type& c)
00241 {
00242 const unsigned int d = vpdt_size(pt);
00243 double det = 1.0;
00244 for (unsigned int i=0; i<d; ++i)
00245 det *= c;
00246 return static_cast<T>(det);
00247 }
00248 };
00249
00250
00251
00252
00253 template<class F>
00254 struct vpdt_norm_metric<F, typename vpdt_field_traits<F>::scalar_type,
00255 typename vpdt_field_traits<F>::type_is_scalar >
00256 {
00257
00258 typedef typename vpdt_field_traits<F>::scalar_type covar_type;
00259
00260 typedef typename vpdt_field_traits<F>::scalar_type T;
00261
00262 typedef typename vpdt_field_traits<F>::vector_type vector;
00263
00264 typedef typename vpdt_field_traits<F>::matrix_type matrix;
00265
00266
00267 static inline T distance(const F& pt1, const F& pt2, const covar_type& c)
00268 {
00269 return vcl_sqrt(sqr_distance(pt1,pt2,c));
00270 }
00271
00272
00273 static inline T sqr_distance(const F& pt1, const F& pt2, const covar_type& c)
00274 {
00275 T tmp = pt1-pt2;
00276 return tmp*tmp/c;
00277 }
00278
00279
00280 static inline T sqr_distance_deriv(const F& pt1, const F& pt2,
00281 const covar_type& c, vector& g)
00282 {
00283 T tmp = pt1-pt2;
00284 g = 2*tmp/c;
00285 return tmp*tmp/c;
00286 }
00287
00288
00289
00290 static inline void compute_covar(matrix& covar, const F& , const covar_type& c)
00291 {
00292 covar = c;
00293 }
00294
00295
00296
00297 static inline T covar_det(const F& , const covar_type& c)
00298 {
00299 return c;
00300 }
00301 };
00302
00303
00304 #endif // vpdt_norm_metric_h_