compute the rigid transformation using SVD with library [Eigen](http://eigen.tuxfamily.org/), usally useful in ICP registration or related works.

computeRigidTransformUsingSVD.cpp

//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::MatrixXd Dt = D.transpose();
	Eigen::Matrix3d H = Dt*S;
	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 Vt = svd.matrixV().transpose();
	Eigen::Matrix3d R = svd.matrixU()*Vt;
	Eigen::Vector3d t = center_dst - R*center_src;	
	
	return std::make_pair(R, t);
}

int main() {
    const int POINT_SIZE = 100;
    srand(time(NULL));
	  while (1) {
  		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();
  	}
    return 0;
}

@iamdpak I ran into a similar problem and I found this paper with some good explanation. To avoid det(R) = -1 , you would negate the third column of V. Note that if you have linear or singular data sets, you can get weird results. The algorithm is designed for planar data sets (more explanation in the paper).

This paper outlines the algorithm above.

nice work! appreciate it!
is there way to figure out accuracy of "RT"?
could i get better accuracy if i do computeTigidTransform more and more?