利用 ransac 算法拟合平面

1.前言

最近项目中遇到一个问题， 老板给了一组数据然后要求获取其中处于同一个平面上的数据点的信息， 很明显就是使用ransac 算法进行处理。

2. ransac算法思想

这里我们使用自己的理解来说明下这个算法。
1. 首先我们从给定的数据集中随机挑选几组数据获取一个模型(最好可以保证随机挑选的数据不重复)
2. 将这个拟合方程作用于所有的数据，根据阈值区分出模型的内点和外点信息
3. 重复多次上述操作，挑选出其中包含内点数目最多的模型

具体的关于ransac 的算法说明可以参考 Random sample consensus

3. 具体实现

3.1 实现效果

3.2 具体代码

SimulatePlane.m

%% 利用ransac 方法从一组三维数据点中找到处于同一个平面上的数据点
% @param W3 输入数据点 4 x N
% @return [W, id_best, plane_best, inliner, outliner] 在同一个平面上的物点,
% 挑选的物点在原始数据集中的标号, 最佳拟合平面表示形式（【a, b, c, d】: ax + by + cz + d = 0, 1 X 4),
% 在平面上的内点数目, 不再平面上的外点数目
%

function [W, id_best, plane_best, inliner_best, outliner_best] = SimulatePlane(W3)
    assert(size(W3, 1) == 4, 'please check your input for W3')

    iterations = 0;
    k = 500;
    W3_len = size(W3, 2);
    dis_limit_xishu = 1.3;
    dis_limit = 0.015;
    inliner_best = 0;

    ids_set = [];
    while iterations < k
        %% 随机 4 点获取平面参数
        while 1
            ids = myrand4(1, W3_len);
            % 去重
            if ~ismember(ids, ids_set, 'rows')
                ids_set = [ids_set; ids];
                break;
            end
        end

        A = W3(:, ids)';
        [uu, dd, vv] = svd(A);
        plane = vv(:,end);
        plane = plane ./ plane(end);
        plane = plane';

        %% 计算内点
        dis = abs(plane * W3);
%         dis_limit_best = min(dis_limit_xishu * mean(dis), dis_limit_best);
%        dis_limit_best = min(dis_limit_xishu * mean(dis), dis_limit);
        dis_limit_best = dis_limit;
        inliner = sum(dis < dis_limit_best);
        outliner = sum(dis >= dis_limit_best);

        if inliner > inliner_best
           plane_best = plane;
           inliner_best = inliner;
           outliner_best = outliner;
        end       

        iterations = iterations + 1;
    end

    %% 找到 best_plane 上的内点
    dis = abs(plane_best * W3);
    id_best = find(dis < dis_limit_best);
    W = W3(:, id_best);
end

function [ids] = myrand4(minlimit, maxlimit)
    ids = [];
    count = 0;
    while count < 4
        id = randi([minlimit, maxlimit], 1, 1);
        if ismember(id, ids)
            continue
        end
        ids = [ids, id];
        count = count + 1;
    end
    ids = sort(ids);
end

相应的图形绘制代码：
参考资料：
http://blog.sina.com.cn/s/blog_5cd4cccf0100q90p.html
http://blog.csdn.net/wangcj625/article/details/6287735/

close all
        figure
        plot3(W3(1, :), W3(2, :), W3(3, :), '*');
        hold on        
        plot3(W(1, :), W(2, :), W(3, :), 'o');
        hold on
        X = min(W3(1, :)) : 0.5 : max(W3(1, :));
        Y = min(W3(2, :)) : 0.5 : max(W3(2, :));
        [x,y] = meshgrid(X, Y);
        z = -(plane(1)*x + plane(2) * y + plane(4))/plane(3);
        surf(x, y, z)
        hold off