Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save wongkaiweng/0339aa069f17525f04e651a35fb52bd5 to your computer and use it in GitHub Desktop.
Save wongkaiweng/0339aa069f17525f04e651a35fb52bd5 to your computer and use it in GitHub Desktop.
compute the rigid transformation using SVD with library [Eigen](http://eigen.tuxfamily.org/), usally useful in ICP registration or related works.
//using Eigen's SVD to fastly compute the rigid transformation between two point clouds.
#include <iostream>
#include <ctime>
#include <Eigen/SVD>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/Geometry>
using namespace Eigen;
using namespace std;
typedef std::pair<Eigen::Matrix3d, Eigen::Vector3d> TransformType;
typedef std::vector<Eigen::Vector3d> PointsType;
TransformType computeRigidTransform(const PointsType& src, const PointsType& dst)
{
assert(src.size() == dst.size());
int pairSize = src.size();
Eigen::Vector3d center_src(0, 0, 0), center_dst(0, 0, 0);
for (int i=0; i<pairSize; ++i)
{
center_src += src[i];
center_dst += dst[i];
}
center_src /= (double)pairSize;
center_dst /= (double)pairSize;
Eigen::MatrixXd S(pairSize, 3), D(pairSize, 3);
for (int i=0; i<pairSize; ++i)
{
for (int j=0; j<3; ++j)
S(i, j) = src[i][j] - center_src[j];
for (int j=0; j<3; ++j)
D(i, j) = dst[i][j] - center_dst[j];
}
Eigen::Matrix3d H = S.transpose()* D;
Eigen::Matrix3d W, U, V;
JacobiSVD<Eigen::MatrixXd> svd;
Eigen::MatrixXd H_(3, 3);
for (int i=0; i<3; ++i) for (int j=0; j<3; ++j) H_(i, j) = H(i, j);
svd.compute(H_, Eigen::ComputeThinU | Eigen::ComputeThinV );
if (!svd.computeU() || !svd.computeV()) {
std::cerr << "decomposition error" << endl;
return std::make_pair(Eigen::Matrix3d::Identity(), Eigen::Vector3d::Zero());
}
Eigen::Matrix3d Ut = svd.matrixU().transpose();
Eigen::DiagonalMatrix<double, 3> M_diag(1,1,(svd.matrixV()*Ut).determinant());
Eigen::Matrix3d R = svd.matrixV()*M_diag*Ut;
Eigen::Vector3d t = center_dst - R*center_src;
return std::make_pair(R, t);
}
int main() {
const int POINT_SIZE = 100;
srand(time(NULL));
int count = 0;
while (count < 10) {
PointsType p1s, p2s;
p1s.resize(POINT_SIZE);
for (int i=0; i<POINT_SIZE; ++i) {
p1s[i][0] = rand()%256*1.0 / 512.0;
p1s[i][1] = rand()%256*1.0 / 512.0;
p1s[i][2] = rand()%256*1.0 / 512.0;
}
TransformType RT;
RT.first = AngleAxisd(rand()%180*1.0, Vector3d::UnitZ())
* AngleAxisd(rand()%180*1.0, Vector3d::UnitY())
* AngleAxisd(rand()%180*1.0, Vector3d::UnitZ());
RT.second = Eigen::Vector3d(3, 4, 1);
std::cout << RT.first << std::endl;
std::cout << (RT.second)[0] << " " << (RT.second)[1] << " " << (RT.second)[2] << endl;
for (int i=0; i<POINT_SIZE; ++i) {
p2s.push_back(RT.first*p1s[i] + RT.second);
}
cout << "computing the rigid transformations...\n";
RT = computeRigidTransform(p1s, p2s);
std::cout << RT.first << endl;
std::cout << (RT.second)[0] << " " << (RT.second)[1] << " " << (RT.second)[2] << endl;
cout << endl;
getchar();
count++;
}
return 0;
}
@wongkaiweng
Copy link
Author

See this paper for the algorithm and this paper on why you need planar data sets.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment