matlab实现ransac算法

RANSAC（Random Sample Consensus）是一种基于统计学原理的参数估计算法，可以用于解决一些具有噪声和异常值的数据拟合问题。在 Matlab 中实现 RANSAC 算法的步骤如下：

1. 定义模型函数

首先需要定义一个模型函数，该函数可以根据数据点集合来计算模型参数。这里以直线拟合为例，模型函数为：

function [k, b] = line_model(data)
% data: 数据点集合，每一行代表一个数据点，第一列为 x 坐标，第二列为 y 坐标
% k, b: 直线方程为 y = kx+b，返回 k 和 b
x = data(:,1);
y = data(:,2);
k = (mean(x.*y) - mean(x)*mean(y)) / (mean(x.^2) - mean(x)^2);
b = mean(y) - k*mean(x);
end

2. 计算模型误差

根据数据点集合和模型参数，计算每个数据点到模型直线的距离作为误差，这里采用欧几里得距离：

function [err] = line_error(data, k, b)
% data: 数据点集合，每一行代表一个数据点，第一列为 x 坐标，第二列为 y 坐标
% k, b: 直线方程为 y = kx+b
% err: 每个数据点到直线的距离
x = data(:,1);
y = data(:,2);
D = abs(y - k*x - b) / sqrt(k^2+1);
err = D;
end

3. RANSAC 迭代

在 RANSAC 迭代过程中，每次从数据点集合中随机选择一定数量的点，计算模型参数并计算误差，然后将误差小于一定阈值的点作为内点，其余点作为外点。迭代过程中记录内点最大的模型参数和内点数量，迭代次数取决于选择的随机点数量、内点数量和置信度。

function [k, b] = ransac(data, iters, n, threshold, confidence)
% data: 数据点集合，每一行代表一个数据点，第一列为 x 坐标，第二列为 y 坐标
% iters: 迭代次数
% n: 随机选择的数据点数量
% threshold: 内点阈值，误差小于该阈值的点为内点
% confidence: 置信度，表示概率上至少有一个模型是正确的
% k, b: 直线方程为 y = kx+b，返回最优的 k 和 b
best_k = 0;
best_b = 0;
max_inliers = 0;
for i = 1:iters
    % 随机选择 n 个数据点
    subset = datasample(data, n, 'Replace', false);
    % 计算模型参数
    [k, b] = line_model(subset);
    % 计算误差
    err = line_error(data, k, b);
    % 标记内点和外点
    inliers = find(err < threshold);
    outliers = find(err >= threshold);
    % 更新最优参数和内点数量
    if length(inliers) > max_inliers
        max_inliers = length(inliers);
        best_k = k;
        best_b = b;
        % 计算迭代次数
        w = max_inliers / size(data, 1);
        p_no_outliers = 1 - w^n;
        iters_required = log(1-confidence) / log(p_no_outliers);
        if i >= iters_required
            break;
        end
    end
end
k = best_k;
b = best_b;
end

4. 测试

使用随机生成的数据测试 RANSAC 算法：

% 生成数据
x = linspace(0, 10, 100);
y = 2*x + 3 + randn(size(x));
data = [x', y'];

% RANSAC 参数
iters = 1000;
n = 2;
threshold = 0.5;
confidence = 0.99;

% RANSAC 迭代
[k, b] = ransac(data, iters, n, threshold, confidence);

% 绘制数据和拟合直线
figure;
plot(x, y, '.');
hold on;
plot(x, k*x+b, 'r');

运行结果如下图所示：