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
00032
00033
00034
00035
00036
00037 #include <Eigen/SVD>
00038 #include <Eigen/LU>
00039 #include <Eigen/Dense>
00040
00041 namespace pcl {
00042
00043 TransformationFromCorrespondences::TransformationFromCorrespondences ()
00044 {
00045 reset();
00046 }
00047
00048 TransformationFromCorrespondences::~TransformationFromCorrespondences ()
00049 {
00050 }
00051
00052 inline void
00053 TransformationFromCorrespondences::reset ()
00054 {
00055 no_of_samples_ = 0;
00056 accumulated_weight_ = 0.0;
00057 mean1_.fill(0);
00058 mean2_.fill(0);
00059 covariance_.fill(0);
00060 }
00061
00062 inline void
00063 TransformationFromCorrespondences::add (const Eigen::Vector3f& point, const Eigen::Vector3f& corresponding_point,
00064 float weight)
00065 {
00066 if (weight==0.0f)
00067 return;
00068
00069 ++no_of_samples_;
00070 accumulated_weight_ += weight;
00071 float alpha = weight/accumulated_weight_;
00072
00073 Eigen::Vector3f diff1 = point - mean1_, diff2 = corresponding_point - mean2_;
00074 covariance_ = (1.0-alpha)*(covariance_ + alpha * (diff2 * diff1.transpose()));
00075
00076 mean1_ += alpha*(diff1);
00077 mean2_ += alpha*(diff2);
00078 }
00079
00080 inline Eigen::Affine3f
00081 TransformationFromCorrespondences::getTransformation ()
00082 {
00083
00084 Eigen::JacobiSVD<Eigen::Matrix<float, 3, 3> > svd (covariance_, Eigen::ComputeFullU | Eigen::ComputeFullV);
00085 const Eigen::Matrix<float, 3, 3>& u = svd.matrixU(),
00086 & v = svd.matrixV();
00087 Eigen::Matrix<float, 3, 3> s;
00088 s.setIdentity();
00089 if (u.determinant()*v.determinant() < 0.0f)
00090 s(2,2) = -1.0f;
00091
00092 Eigen::Matrix<float, 3, 3> r = u * s * v.transpose();
00093 Eigen::Vector3f t = mean2_ - r*mean1_;
00094
00095 Eigen::Affine3f ret;
00096 ret(0,0)=r(0,0); ret(0,1)=r(0,1); ret(0,2)=r(0,2); ret(0,3)=t(0);
00097 ret(1,0)=r(1,0); ret(1,1)=r(1,1); ret(1,2)=r(1,2); ret(1,3)=t(1);
00098 ret(2,0)=r(2,0); ret(2,1)=r(2,1); ret(2,2)=r(2,2); ret(2,3)=t(2);
00099 ret(3,0)=0.0f; ret(3,1)=0.0f; ret(3,2)=0.0f; ret(3,3)=1.0f;
00100
00101 return ret;
00102 }
00103
00104 }