多点拟合求平面的RANSAC算法

对于多点求平面,直接的传统方法是三点求平面,通过求解方程组得到平面的法向量和截距。

但仅适用于点集能够恰好拟合到一个平面的情况,如果点集不能恰好拟合到一个平面,则需要使用更复杂的方法来处理,例如使用最小二乘平面拟合或者使用RANSAC算法。
利用所有点到目标平面距离之和最小的条件,筛选出目标平面。

RANSAC是“Random Sample Consensus(随机抽样一致)”的缩写。它可以从一组包含“局外点”的观测数据集中,通过迭代方式估计数学模型的参数。它是一种不确定的算法,有一定的概率得出一个合理的结果。为了提高得出合理结果的概率必须提高迭代次数。

1、基本思想:

RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点,并用下述方法进行验证:

有一个模型适用于假设的局内点,即所有的未知参数都能从假设的局内点计算得出。

用1中得到的模型去测试所有的其它数据,如果某个点适用于估计的模型,认为它也是局内点。

如果有足够多的点被归类为假设的局内点,那么估计的模型就足够合理。

然后,用所有假设的局内点去重新估计模型,因为它仅仅被初始的假设局内点估计过。

最后,通过估计局内点与模型的错误率来评估模型。

2、对上述步骤,进行简单总结如下:

N:样本点个数,K:求解模型需要的最少的点的个数

随机采样K个点
针对该K个点拟合模型
计算其它点到该拟合模型的距离,小于一定阈值当做内点,统计内点个数
重复M次,选择内点数最多的模型
利用所有的内点重新估计模型(可选)
距离代码

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using System.Windows.Media.Media3D;

namespace 平面拟合
{
    public partial class Form1 : Form
    {
        List<Point3D> AllPoints;
        public Form1()
        {
            InitializeComponent();
        }
        /// <summary>
        /// 点到直线的距离
        /// </summary>
        /// <param name="point"></param>
        /// <param name="plane"></param>
        /// <returns></returns>
        double PointToPlaneDis(Point3D point, Plane plane)
        {
            var p = point;
            return Math.Abs(plane.Normal.X * p.X + plane.Normal.Y * p.Y + plane.Normal.Z * p.Z + plane.Intercept) /
                         Math.Sqrt(plane.Normal.X * plane.Normal.X + plane.Normal.Y * plane.Normal.Y + plane.Normal.Z * plane.Normal.Z);
        }

        /// <summary>
        /// 随机抽样法多点拟合平面
        /// </summary>
        /// <param name="points"></param>
        /// <param name="iterations">迭代计算次数一般大于点数量</param>
        /// <param name="threshold">筛选条件点到平面的距离</param>
        /// <returns></returns>
        public Plane FitPlaneRANSAC(List<Point3D> points, int iterations, double threshold)
        {
            // 按照RANSAC算法迭代次数进行循环
            int maxInliers = 0;
            Plane bestPlane = null;
            List<int> recordIdx = new List<int> { 0, 0, 0 };
            List<List<int>> recordIdxList = new List<List<int>>();
            // 随机选择三个点构成一个平面
            Random rand = new Random();
            for (int i = 0; i < iterations; i++)
            {
               
                int index1 = rand.Next(points.Count);
                int index2 = rand.Next(points.Count);
                int index3 = rand.Next(points.Count);
                while (index1 == index2 || index1 == index3 || index2 == index3)
                {
                    index1 = rand.Next(points.Count);
                    index2 = rand.Next(points.Count);
                    index3 = rand.Next(points.Count);
                }

                //检测是否已经计算过该平面
                foreach (var record in recordIdxList)
                {
                    if (record.Contains(index1) && record.Contains(index2) && record.Contains(index3))
                    {                     
                        continue;
                    }
                }
                //记录已经计算的平面
                recordIdxList.Add(new List<int> { index1, index2, index3 });

                Point3D p1 = points[index1];
                Point3D p2 = points[index2];
                Point3D p3 = points[index3];
                Plane plane = new Plane(p1, p2, p3);

                // 计算与该平面距离小于阈值的点数
                int inliers = 0;
                foreach (Point3D p in points)
                {
                    double distance = PointToPlaneDis(p, plane);
                    if (distance <= threshold)
                    {
                        inliers++;
                    }
                }

                // 如果当前平面的支持点数比之前的最好平面要多,则更新最好平面
                if (inliers > maxInliers)
                {
                    maxInliers = inliers;
                    bestPlane = plane;
                }
            }

            return bestPlane;
        }

        private async void button1_Click(object sender, EventArgs e)
        {         
            Plane bestPlane = null;
            Task tk = Task.Factory.StartNew(() => {
                AllPoints = new List<Point3D>();
                int pointcnt = 300;
                Random rand = new Random();

                //假设平面5x+4y+z+1=0;生成距离平面距离10以内的300个点
                //使用这300个点,进行随机平面拟合匹配验证结果和目标平面对比
                while (pointcnt > 0)
                {
                   
                    double X = rand.Next(1, 100);
                    double Y = rand.Next(1, 100);           
                   var mSqrt =  Math.Sqrt(5 * 5 + 4 * 4 + 1);
                    //生成随机距离
                    double randDisZ = (rand.Next(-100, 100)) /10;
                    //利用距离反求出点的坐标Z
                    double Z = randDisZ * mSqrt-5*X-4*Y-1;

                    AllPoints.Add(new Point3D(X, Y, Z));
                    pointcnt--;
                }

                bestPlane = FitPlaneRANSAC(AllPoints, 9900,10);


            });
         
            await tk;
            Vector3D p = bestPlane.Normal;
            var D = bestPlane.Intercept;
            //把Ax+By+Cz+D=0转化两边同时除以C,消去Z的系数
            MessageBox.Show($"求出直线{p.X/ p.Z}X+{p.Y/ p.Z}Y + Z+{D/ p.Z}=0");

        }
    }
    public class Plane
    {
        public Vector3D Normal { get; set; }
        public double Intercept { get; set; }
        public List<Point3D> pointXyzList
        {
            get
            {
                if (Normal.Z != 0 && Normal.X != 0 && Normal.Y != 0)
                    return new List<Point3D>(){new Point3D() { X = 0,Y = 0,Z = -Intercept / Normal.Z},
                    new Point3D() { Z = 0,Y = 0,X = -Intercept / Normal.X},
                    new Point3D() { Z = 0,X = 0,Y = -Intercept / Normal.Y} };
                else
                {
                    return new List<Point3D>() { new Point3D() { X = 0, Y = 0, Z = 0 } };
                }
            }
        }

        public Plane(Point3D p1, Point3D p2, Point3D p3)
        {
            // 计算法向量
            Vector3D v1 = new Vector3D(p2.X - p1.X, p2.Y - p1.Y, p2.Z - p1.Z);
            Vector3D v2 = new Vector3D(p3.X - p1.X, p3.Y - p1.Y, p3.Z - p1.Z);
            Normal = Vector3D.CrossProduct(v1, v2);
            Normal.Normalize();

            // 计算截距
            Intercept = -(Normal.X * p1.X + Normal.Y * p1.Y + Normal.Z * p1.Z);
        }
       
    }
}

该代码首先通过循环执行一定次数的迭代,每次迭代随机选择三个点构成一个平面,计算与该平面距离小于阈值的点数,然后选取支持点数最多的平面作为最好平面。在实际应用中,需要根据具体情况调整迭代次数和阈值的取值。