Structural Bioinformatics Library
Template C++ / Python API for developping structural bioinformatics applications.
User Manual

Point_cloud_rigid_registration_3

Authors: F. Cazals and T. Dreyfus

Introduction

The point set registration consists on finding a spatial transformation aligning a set of points onto another one. The registration is said rigid if the distance between all pairs of points is conserved by the transformation. This package defines a rigid registration between two set of 3D points where the spatial transformation is a composition of translations and rotations.

Implementation

The rigid registration is done using a linear algebra algorithm that computes the optimal rotation matrix between the two input set of points using the Singular Value Decomposition (SVD) of their covariance matrix. Basically, once centered on the origin, the covariance matrix $ C $ of the two input sets of points $ X $ and $ Y $ is computed. The SVD of $ C = V * S W^t $ then produces orthogonal matrices representing rotations ( $ V, W $) and deformations ( $ S $) from $ X $ to $ Y $. By "removing" the deformation, it is then possible to compute the optimal rotation for aligning both set of points.

The algorithm then consists :

  • 1) computes the centroids of $ X $ and $ Y $ to center the input set of points, and fills the matrices $ M_X $ and $ M_Y $ of size $ (3,N) $ and $ (3,M) $ with the new centered points;
  • 2) computes the covariance matrix $ C = M_X * M_Y^t $;
  • 3) computes the SVD of $ C = V * S * W^t $;
  • 4) computes the matrix $ E $ that is the identity matrix $ (3, 3) $, except that $ E(2, 2) = \mid C\mid $, and then computes the optimal rotation matrix $ U = W * E * V^t $;
  • 5) rotates the matrix $ M_X $ using $ U $ and translates the resulting points using the centroid of $ Y $;

Design

The algorithm is coded in the class SBL::GT::T_Point_cloud_rigid_registration_3< FT >, where FT is the number type used for representing the coordinates of the points (by default, the double type). The class is a functor that takes as input the two set of points and fills an output container of 3D points. The type of the output 3D points has to be constructible from the type of the input 3D points. In other words, if the input and output points have the same type, it should be copy constructible. There are three way to call the functor :

  • by providing two random-access containers of points (e.g arrays or vectors),
  • by providing two pairs of (begin, end) iterators over the containers of points,
  • by providing two triples (size, begin, end), corresponding to the size and iterators over containers of the coordinates of the input points.

The last one is useful when the input sets of points are represented by D-dimensional points, for exemple when computing the distance between molecular conformations (see package Molecular_distances)

Examples

Registration with rotated and translated set of points

The following example makes a registration of a set of points with the same set of points, but translated and rotated. It then computes the distance between pairs of corresponding points, computing the so called l-RMSD (see package Molecular_distances).

#include <iostream>
#include <fstream>
#include <SBL/GT/Point_cloud_rigid_registration_3.hpp>
#include <CGAL/Cartesian.h>
typedef CGAL::Cartesian<double> K;
typedef K::Point_3 Point_3;
int main(int argc, char *argv[])
{
//Read 3D points from an input file
if(argc < 2) return -1;
std::vector<Point_3> points;
std::ifstream in(argv[1]);
while(!in.eof())
{
Point_3 p; in >> p;
points.push_back(p);
}
in.close();
Point_cloud_rigid_registration_3 registration;
std::vector<K::FT> points_registered;
std::vector<Point_3> points_transformed;
//Rotates and translates the input points
for(unsigned i = 0; i < points.size(); i++)
{
K::FT x = std::cos(CGAL_PI/4)*points[i].x() - std::sin(CGAL_PI/4)*points[i].y();
K::FT y = std::sin(CGAL_PI/4)*points[i].x() + std::cos(CGAL_PI/4)*points[i].y();
points_transformed.push_back(Point_3(x + 1, y, points[i].z()));
}
//Makes the registration
registration(points, points_transformed, std::back_inserter(points_registered));
//Comptes the distance between pairs of 3D points
K::FT d = 0;
for(unsigned i = 0;i < points.size(); i++)
d += (points[i] - Point_3(points_registered[3*i], points_registered[3*i + 1], points_registered[3*i + 2])).squared_length()/(K::FT)points.size();
std::cout << "Rotated and Translated Registration: " << d << std::endl;
return 0;
}

Applications

This package offers alors a program for performing the rigid registration between a reference set of 3D points and a collection of sets of 3D points : sbl-align-point-clouds-3.exe . An input set of 3D points is represented as a D dimensional point, where D is divisable by 3. A file listing D dimensional points is a text file where each point is represented by its dimension followed by its coordinates. Note that breaklines are not considered, so that having one D dimensional point per line, or one 3D point per line are both acceptable :

6 0 0 0 1 1 1

or

6
0 0 0
1 1 1

The program sbl-align-point-clouds-3.exe takes two files as input : the first one containing one D dimensional point representing the reference set of points, and the second containing any number of D dimensional points.