:RANSAC算法)
NSACRandom Sample Consensus是一种鲁棒的参数估计方法特别适用于数据中包含大量离群点outliers的场景。其核心思想是通过反复随机采样一小部分数据来生成候选模型然后评估该模型支持的内点数量最终选择内点最多的模型。下面给出RANSAC的数学推导涵盖迭代次数计算、内点判定准则以及算法收敛性分析。 RANSAC的核心思想是随机选取一小部分点最小子集来生成候选模型然后用整个数据集验证该模型统计符合该模型的点内点数量。重复多次选择内点最多的模型作为最终估计。1. 问题定义与数学符号 给定样本散点数据集D{xi}Ni1xi 可以表示空间点坐标或对应关系等。我们希望估计一个参数为θ∈Rd的模型例如直线、平面、基本矩阵等。模型可由最小的m个数据点唯一确定例如直线模型2点、平面3个点。数据集D包含两类点内点inlier points符合真实模型在允许误差内的数据点假设其噪声服从某种分布如高斯分布。离群点outlier points不符合模型的点可能来自错误测量、其他结构等。内点率ϵ 数据中内点所占比例未知。RANSAC 的目标是找到一组参数θ∗使得内点数量最大或某种代价函数最小。 RANSAC的核心思想是随机选取一小部分点最小子集来生成候选模型然后用整个数据集验证该模型统计符合该模型的点内点数量。重复多次选择内点最多的模型作为最终估计。2. RANSAC算法步骤输 入样本数据集D 最小采样点数m距离阈值τ置信度p (如0.99)最大迭代次数Kmax 。输 出最优化模型参数θfinal 及内点集I∗。初始化令最优点数nbest0, 迭代次数 k0。随机采样从D 中均匀随机抽取m 个不同点构成子集Sk。模型假设由Sk 计算模型参数θk (通常通过最小二乘PCA等其他方法或直接求解优化)。内点计数 对每个xi∈D计算距离did(xi,θk) 。若di≤τ则将xi判为内点。记内点集Ik{xi:diτ}内点数集nk|Ik|.更新最优若nknbest 则令θ∗θk,I∗Iknbestnk并根据当前nbest 更新所需迭代次数K。迭代优化k←k1若kK且未达到其他终止条件返回步骤2;最终优化使用I∗ 中的所有内点重新估计模型得到最终参数θfinal。3.迭代次数K 的概率推导3.1 基本公式 设数据中内点的真实比例为ε。一次采用中m 个点全部为内点的概率为pεm若每次采用独立则经过K次采用后至少有一次全内点采样的概率即成功概率Psuccess为Psuccess1−(1−p)K我们希望Psuccess 不低于某个置信度pconf 如0.99。解得所需最少迭代次数K[log(1−pconf)log(1−εm)]此公式是RANSAC 理论的核心。注意ε 通常未知需估计。3.2 自适用迭代次数 实际中我们不知道ε 但可以在迭代过程中根据当前找到的最佳内点数nbest 来估计内点率ˆεnbestN那么一次采样迭代次数全内点的概率估计为ˆpˆεm。为保证置信度pconf 所需迭代次数为Knew[log(1−pconf)log(1−ˆp)]若当前已迭代次数kKnew则继续否则可提前终止这种方法避免了固定K的盲目性。3.3 概率下界分析当ε 很小K 可能性极大。例如ε0.1,m3则p0.001,K≈log(0.01)/log(0.999)≈4603 。因此RANSAC 适用于内点率不太低的问题。4. 模型质量度量 RANSAC 用内点数量作为模型质量的指标这等价于一个0-1损失函数L(xi,θ){0,d(xi,θ)τ1,d(xi,θ)≥τ总损失为离群点个数最大化内点数量即最小化离群点个数。这种硬阈值损失虽然简单但忽略了距离的具体数值。6.最终模型优化 找到最佳内点集I∗ 后为了获取更精确的参数估计通常使用使用内点重新进行模型拟合如最小二乘。若内点噪声服从独立同分布的高斯分布则最大似然估计等价于最小化平方距离和θfinalargmin这一步显著提高了参数精度因为最初只有m 个点而现在利用了所有内点。7. RANSAC 平面拟合案例 RANSAC通过反复随机采样一小部分点最小子集来生成候选平面然后用整个点云评估该平面的内点数量即到平面距离小于阈值的点数。重复多次后选择内点最多的平面作为最佳模型最后用所有内点重新优化平面参数如使用PCA。对于平面拟合最小子集大小为 m3因为三点确定一个平面。内点判定基于点到平面的垂直距离。7.1 算法步骤输 入散点点云 \mathcal{P}\{\mathbf{p}_i\}_{i1}^N 距离阈值\tau置信度p_{conf}通常0.99。输 出最优平面参数法向量\mathbf{n} 和常数d及内点集的索引。初始化令最优内点数n_{best}0迭代次数k0最大迭代次数K\infty或是一个大数。while (kK) doa.随机采样从\mathcal{P}中随机抽取3个不共线的点\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3 。b.计算平面模型由三点确定平面法平面\mathbf{n}_k(\mathbf{p}_2-\mathbf{p}_1)\times(\mathbf{p}_3-\mathbf{p}_1) 归一化; 平面方程 \mathbf{n}_k^{T}\mathbf{x}d_k0其中d_k-\mathbf{n}_k\mathbf{p}_1。c.内点计数对每个点\mathbf{p}_i ,计算距离d_i|\mathbf{n}_k^T\mathbf{p}_id_k| 因为\mathbf{n}_k 是单位向量此即垂直距离。若d_i\leq{\tau} 则判为内点。记内点集\mathcal{I}_k内点数n_k|\mathcal{I}_k|。d.更新最优若n_kn_{best} ,则更新最优模型(\mathbf{n}^*,d^*)(\mathbf{n}_k,d_k)\mathcal{I}^*\mathcal{I}_k,n_{best}n_k并根据当前内点率\hat\varepsilon\frac{n_{best}}{N} 重新计算所需的迭代次数K\frac{\log(1-p_{conf})}{\log(1-\hat{\varepsilon}^3)}\tag{8}e.kk1最终优化使用最优内点集\mathcal{I}^* 重新拟合平面(例如通过PCA或最小二乘 )、得到更精确的参数(\mathbf{n}_{final},d_{final}) 。返回\mathbf{n}_{final},d_{final}及\mathcal{I}^{*}。注意在实际应用中需要考虑设置最大迭代次数K_{max} 。7.2 RANSAC拟合平面的 matlab 仿真代码 以下代码生成一个带有大量离群点的平面点云用RANSAC拟合平面并与标准PCA拟合结果对比。文件一RANSAC的平面拟合算法函数function [a,b,c,d]ransac_fitplane(data,N,t,T) figure;plot3(data(1,:),data(2,:),data(3,:),o);hold on; % 显示数据 iter N; %抽样次数N number size(data,2); % 总数 maxNum0; %符合拟合模型的数据的个数 for i1:iter %循环次数 sampleidx randperm(number); sampleidx sampleidx(1,1:3); sample data(:,sampleidx); %取样本 [a1,a2,a3,a4]get_nice_plane(sample);%拟合直线方程 zaxbyc plane [-a1/a3,-a2/a3,-1,-a4/a3];%建模型 maskabs(plane*[data; ones(1,size(data,2))]); %求每个数据到拟合平面的距离 totalsum(maskt); %计算数据距离平面小于一定阈值的数据的个数 index maskt; if totalT nsampledata(:,index); [a1,a2,a3,a4]get_nice_plane(nsample); plane [-a1/a3,-a2/a3,-1,-a4/a3]; %zaxbyc maskabs(plane*[data; ones(1,size(data,2))]); totalsum(maskt); %计算数据距离平面小于一定阈值的数据的个数 end; if totalmaxNum %记录最大total maxNumtotal; bestplaneplane;%最好的拟合平面 bestindexindex; bestplane2[a1,a2,a3,a4]; end end %显示符合最佳拟合的数据 maxNum %最大一致集 avgerrorabs(bestplane*[data; ones(1,size(data,2))]); avgerrorsum(avgerror(find(avgerrort)))/maxNum %计算一致集内平均误差 abestplane2(1);bbestplane2(2);cbestplane2(3);dbestplane2(4); % 图形绘制 temp1data(1,bestindex); temp2data(2,bestindex); xfit min(temp1):0.2:max(temp1); yfit min(temp2):0.2:max(temp2); [XFIT,YFIT] meshgrid (xfit,yfit); ZFIT bestplane(1)*XFITbestplane(2)*YFITbestplane(4); mesh(XFIT,YFIT,ZFIT);grid on; xlabel(X); ylabel(Y); end文件二利用PCA拟合平面的函数function [a,b,c,d]get_nice_plane(data) planeDatadata; % 协方差矩阵的SVD变换中最小奇异值对应的奇异向量就是平面的方向 xyz0mean(planeData,1); centeredPlanebsxfun(minus,planeData,xyz0); [~,~,V]svd(centeredPlane); aV(1,3); bV(2,3); cV(3,3); d-dot([a b c],xyz0); end文件三测试主程序mu[0 0 0]; S[2 0 4;0 4 0;4 0 8]; data1mvnrnd(mu,S,400); %外点 mu[2 2 2]; S[8 1 4;1 8 2;4 2 8]; data2mvnrnd(mu,S,500); data[data1,data2];%% 相当于原始数据 [a,b,c,d]ransac_fitplane(data,3000,0.5,300)代码说明数据生成内点来自真实平面加小高斯噪声离群点在大范围内均匀随机分布比例较高。RANSAC主循环每次随机选取抽取3个点计算平面。用所有点计算到平面的距离统计内点。若内点创新点创新高则更新最优模型并根据新内点率动态调整所需迭代次数K。循环直到达到K 或最大迭代次数。