Go to the documentation of this file.00001 #ifndef OPENTISSUE_CORE_MATH_MATH_COVARIANCE_H
00002 #define OPENTISSUE_CORE_MATH_MATH_COVARIANCE_H
00003
00004
00005
00006
00007
00008
00009
00010 #include <OpenTissue/configuration.h>
00011
00012 namespace OpenTissue
00013 {
00014
00015 namespace math
00016 {
00017
00023 template< typename vector3_iterator, typename vector3_type, typename matrix3x3_type>
00024 inline void
00025 covariance(
00026 vector3_iterator begin, vector3_iterator end
00027 , vector3_type & mean, matrix3x3_type & C
00028 )
00029 {
00030 typedef typename vector3_type::value_traits value_traits;
00031
00032 unsigned int N = 0;
00033
00034 mean.clear();
00035
00036 for(vector3_iterator v = begin;v!=end;++v,++N)
00037 mean += (*v);
00038 mean /= N;
00039
00040 C.clear();
00041 for(vector3_iterator v = begin;v!=end;++v,++N)
00042 {
00043 C(0,0) += ((*v)(0) - mean(0))*((*v)(0) - mean(0));
00044 C(1,1) += ((*v)(1) - mean(1))*((*v)(1) - mean(1));
00045 C(2,2) += ((*v)(2) - mean(2))*((*v)(2) - mean(2));
00046 C(0,1) += ((*v)(0) - mean(0))*((*v)(1) - mean(1));
00047 C(0,2) += ((*v)(0) - mean(0))*((*v)(2) - mean(2));
00048 C(1,2) += ((*v)(1) - mean(1))*((*v)(2) - mean(2));
00049 }
00050 C(1,0) = C(0,1);
00051 C(2,0) = C(0,2);
00052 C(2,1) = C(1,2);
00053 C /= N;
00054 }
00055
00062 template<typename vector3_type, typename matrix3x3_type>
00063 inline void
00064 covariance_union(
00065 vector3_type const & mean1, matrix3x3_type const & C1
00066 , vector3_type const & mean2, matrix3x3_type const & C2
00067 , vector3_type & mean, matrix3x3_type & C
00068 )
00069 {
00070 typedef typename vector3_type::value_traits value_traits;
00071
00072 mean = (mean1 + mean2)/value_traits::two();
00073 matrix3x3_type KK = outer_prod<matrix3x3_type,vector3_type>(mean,mean);
00074 matrix3x3_type NN = outer_prod<matrix3x3_type,vector3_type>(mean1,mean1);
00075 matrix3x3_type MM = outer_prod<matrix3x3_type,vector3_type>(mean2,mean2);
00076 C = ((C1 + NN + C2 + MM)/value_traits::two() - KK) ;
00077 }
00078
00079 }
00080
00081 }
00082
00083
00084 #endif