两组3D点结算两坐标系的R和T

最新推荐文章于 2023-08-14 23:36:41 发布

原创最新推荐文章于 2023-08-14 23:36:41 发布 · 1k 阅读

3 ·

CC 4.0 BY-SA版权

文章标签：

#3d #c++

GSP 专栏收录该内容

5 篇文章

订阅专栏

每组最少三个点解两坐标系的R和T，采用论文《Least-Squares Fitting of Two 3-D Point Sets》的方法，参照nghiaho的博客：http://nghiaho.com/?page_id=671，及gihub的MATLAB代码，github地址为：https://github.com/nghiaho12/rigid_transform_3D，用C++及Eigen库实现了此功能（考虑了两组点尺度不同问题），记录一下。

#include <pangolin/pangolin.h>
#include <Eigen/Core>
#include <unistd.h>
#include<iostream>
#include<fstream>
#include<algorithm>
#include<chrono>
#include<opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <math.h>


using namespace std;
using namespace cv;
using namespace Eigen;

  

int main(int argc, char **argv) {

  //vector<Isometry3d, Eigen::aligned_allocator<Isometry3d>> poses;
   //Eigen::Quaterniond q;
    //MatrixXd t = MatrixXd::Random(3,1);
    Vector3d t = Vector3d::Random(),error = Vector3d::Random();
    Matrix3d R, R_estimate,V,U;


    R << -0.6208, 0.7763,-0.1091,
         -0.2820,-0.3509,-0.8929,
         -0.7315,-0.5236, 0.4368;
    //cout<<R.determinant()<<endl;
         
    int n = 10;
    //A: 10*3,  B: 3*10
    MatrixXd A = 100*MatrixXd::Random(n,3),A_move,B_move,A_norm,B_norm,lam2;
    MatrixXd B = 100*R*A.transpose()+(100*t).replicate(1,n)+error.replicate(1,n);
    
    Vector3d centroid_A = A.colwise().sum()/n;
    Vector3d centroid_B = B.rowwise().sum()/n;
    //cout<<"centroid_A = \n"<<centroid_A<<endl;
    //cout<<"centroid_B = \n"<<centroid_B<<endl;

    
    A_move = A.transpose() - centroid_A.replicate(1,n);
    B_move = B - centroid_B.replicate(1,n);
    A_norm = A_move.cwiseProduct(A_move).colwise().sum();
    B_norm = B_move.cwiseProduct(B_move).colwise().sum();
    //cout<<"A_norm.array():\n"<<A_norm.array()<<endl;
    lam2 = A_norm.cwiseQuotient(B_norm);
    double lam_2 = sqrt(lam2.mean());
    cout<<"lam_2:\n"<<1/lam_2<<endl;
    
    

    Matrix<double,3,3> H = (A.transpose()-centroid_A.replicate(1,n))*((lam_2*B-(lam_2*centroid_B).replicate(1,n)).transpose());
    //cout<<"H = "<<H<<endl;
    //SVD decomposition
    JacobiSVD<MatrixXd> svd(H, ComputeThinU | ComputeThinV);
    //cout<<"svd.matrixU():\n"<<svd.matrixU()<<endl;
    //cout<<"svd.matrixV():\n"<<svd.matrixV()<<endl;
    //cout<<"svd.singularValues():\n"<<svd.singularValues()<<endl;

    U = svd.matrixU();
    V = svd.matrixV();
    R_estimate = V*U.transpose();
    if(R_estimate.determinant()<0)
    {
        cout<<"Reflection detected\n"<<endl;
        V.col(2) = -1*V.col(2);
        R_estimate = V*U.transpose();
    }
    t = -R_estimate*centroid_A + lam_2 * centroid_B;
    
    Matrix<double,3,10> B2,err;
    double err_sum,rmse;
    
    B2 = R_estimate*A.transpose()+t.replicate(1,n);
    err = B2 -lam_2 *B;
    err =  err.cwiseProduct(err);
    err_sum = err.sum();
    rmse = sqrt(err_sum);
    cout<<"rmse:\n"<<rmse<<endl;
    return 0;
}

运行结果：