ICP point cloud registration

对原文（ICP配准MATLAB实现）稍作整理，方便自己观看。

        相机由一个视角采集到的每帧点云数据通常仅仅包含所测物体的部分点云信息。因此，为了获取所测物体完整的点云数据，往往需要对物体进行多站点多视角测量，并通过配准把不同视角测得的点云数据变换到同一坐标系下。 

创建rotate.m函数计算平移旋转

function [data_q,T] = rotate(data,x,y,z,t)

%欧拉角转旋转矩阵

x = x/180*pi;
y = y/180*pi;
z = z/180*pi;

Rx = [1      0      0;
    0 cos(x) -sin(x);
    0 sin(x) cos(x)];
Ry = [cos(y)  0 sin(y);
    0       1      0;
    -sin(y) 0 cos(y)];
Rz = [cos(z) -sin(z) 0;
    sin(z) cos(z)  0;
    0      0       1];
T = Rz*Ry*Rx;        %旋转矩阵

T = [T(1,1),T(1,2),T(1,3),t(1);
    T(2,1),T(2,2),T(2,3),t(2);
    T(3,1),T(3,2),T(3,3),t(3);
    0 0 0 1];  % 复合

rows=size(data,2);
rows_one=ones(1,rows);
data=[data;rows_one];    %化为齐次坐标

data_q=T*data;
data_q=data_q(1:3,:);    %返回三维坐标

ICP配准代码

% ICP 算法
clear;
close all;
clc;

data_source=load('could.txt');
data_source=data_source';
%旋转角度
alpha = 20;
theta = 5;
grama = 5;

t=[1,1,30];   %平移向量

[data_target,T0]=rotate(data_source,alpha ,theta,grama,t);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%绘制两幅原始图像
x1=data_source(1,:);
y1=data_source(2,:);
z1=data_source(3,:);
x2=data_target(1,:);
y2=data_target(2,:);
z2=data_target(3,:);
figure(1);
title('初始位置');
scatter3(x1,y1,z1,'b*');
hold on;
scatter3(x2,y2,z2,'r*');
hold off;


T_final=eye(4,4);   %旋转矩阵初始值
iteration=0;
Rf=T_final(1:3,1:3);
Tf=T_final(1:3,4);
data_target=Rf*data_target+Tf*ones(1,size(data_target,2));    %初次更新点集（代表粗配准结果）
err=1;
while(err>0.0001)
    iteration=iteration+1;    %迭代次数
    disp(['迭代次数ieration=',num2str(iteration)]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %利用欧式距离找出对应点集 搜索每个点
    k=size(data_target,2);
    for i = 1:k
        data_q1(1,:) = data_source(1,:) - data_target(1,i);    % 两个点集中的点x坐标之差
        data_q1(2,:) = data_source(2,:) - data_target(2,i);    % 两个点集中的点y坐标之差
        data_q1(3,:) = data_source(3,:) - data_target(3,i);    % 两个点集中的点z坐标之差
        distance = data_q1(1,:).^2 + data_q1(2,:).^2 + data_q1(3,:).^2;  % 欧氏距离
        [min_dis, min_index] = min(distance);   % 找到距离最小的那个点
        data_mid(:,i) = data_source(:,min_index);   % 将那个点保存为对应点
        error(i) = min_dis;     % 保存距离差值
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %去中心化
    % 质心
    data_target_mean=mean(data_target,2);
    data_mid_mean=mean(data_mid,2);
    
    data_target_c=data_target-data_target_mean*ones(1,size(data_target,2));
    data_mid_c=data_mid-data_mid_mean*ones(1,size(data_mid,2));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %SVD分解
    W=zeros(3,3);
    % 计算协方差矩阵
    for j=1:size(data_target_c,2)
        W=W+data_mid_c(:,j)*data_target_c(:,j)';
    end
    [U,S,V]=svd(W);
    Rf=U*V';
    % 计算质心平移矢量
    Tf=data_mid_mean-Rf*data_target_mean;
    
    err=mean(error);
    
    T_t=[Rf,Tf];
    T_t=[T_t;0,0,0,1];
    T_final=T_t*T_final;   %更新变换矩阵
    disp(['误差err=',num2str(err)]);
    disp('变换矩阵T=');
    disp(inv(T_final));
    
    data_target=Rf*data_target+Tf*ones(1,size(data_target,2));    %更新点集
    if iteration >= 200
        break
    end
end

disp('真值');
disp(T0);  %旋转矩阵真值,对矩阵求逆

x1=data_source(1,:);
y1=data_source(2,:);
z1=data_source(3,:);
x2=data_target(1,:);
y2=data_target(2,:);
z2=data_target(3,:);
figure(2);
title('配准后');
%scatter3(x1,y1,z1,'r*');
scatter3(x1,y1,z1,'r.');
hold on;
%scatter3(x2,y2,z2,'b*');
scatter3(x2,y2,z2,'b.');
%grid off;  %关闭网格
xlabel('X(m)');
ylabel('Y(m)');
zlabel('Z(m)');
hold off;

其中，could.txt为自己的需要配准的两组点云数据。