隨機抽樣一致性算法（RANSAC)示例及源代碼

作者：王先榮

    大約在兩年前翻譯了《隨機抽樣一致性算法RANSAC》，在文章的最后承諾寫該算法的C#示例程序。可惜光陰似箭，轉眼許久才寫出來，實在抱歉。本文將使用隨機抽樣一致性算法來來檢測直線和圓，並提供源代碼下載。

一、RANSAC檢測流程

    在這里復述下RANSAC的檢測流程，詳細的過程見上一篇翻譯文章：

    RANSAC算法的輸入是一組觀測數據，一個可以解釋或者適應於觀測數據的參數化模型，一些可信的參數。

    RANSAC通過反復選擇數據中的一組隨機子集來達成目標。被選取的子集被假設為局內點，並用下述方法進行驗證：
    1.有一個模型適應於假設的局內點，即所有的未知參數都能從假設的局內點計算得出。
    2.用1中得到的模型去測試所有的其它數據，如果某個點適用於估計的模型，認為它也是局內點。
    3.如果有足夠多的點被歸類為假設的局內點，那么估計的模型就足夠合理。
    4.然后，用所有假設的局內點去重新估計模型，因為它僅僅被初始的假設局內點估計過。
    5.最后，通過估計局內點與模型的錯誤率來評估模型。
    這個過程被重復執行固定的次數，每次產生的模型要么因為局內點太少而被舍棄，要么因為比現有的模型更好而被選用。

二、得到觀測數據

    我們沒有實驗（測試）數據，這里用手工輸入的數據來替代——記錄您在PictureBox中的點擊坐標，作為觀測數據。

        /// <summary>
        /// 得到樣本點
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void pbSample_Click(object sender, EventArgs e)
        {
            MouseEventArgs me=(MouseEventArgs)e;
            txtRandomPoints.Text += string.Format("({0},{1}),", me.X, me.Y);
            DrawPoint(new Point(me.X, me.Y));
        }

 

三、檢測直線

    3.1 直線的相關知識

    （1）平面上的任意兩點可以確定一條直線；

    （2）直線的通用數學表達形式為：ax+by+c=0。這種表達形式有三個未知數，需要提供三個點才能解出a,b,c三個參數。由於隨機選擇的三個點不一定在一條直線上，所以程序中放棄這種方式。

    （3）直線可以用y=ax+b及x=c這兩個式子來表示。這兩種形式只有一個或者兩個未知數，只需兩個點就能解出a,b,c三個參數。隨機選擇的兩個點即可得到直線，我們采用這種形式。

    3.2 直線類

    直線類(Line)封裝了跟直線相關的一些屬性及方法，列表如下：

    （1）屬性

    A——y=ax+b中的a

    B——y=ax+b中的b

    C——x=c中的c

    （2）構造函數

    public Line(PointF p1, PointF p2)

    提供兩個點p1及p2，計算出直線的屬性A,B,C。

    （3）方法

    GetDistance——獲取點到直線之間的距離；

    GetY——根據x坐標，獲取直線上點的y坐標；

    ToString——獲取直線的方程式。

using System;
using System.Collections.Generic;
using System.Text;
using System.Drawing;
using System.Drawing.Drawing2D;

namespace Ransac
{
    /// <summary>
    /// 直線類：采用y=a*x+b或者x=c的形式表示直線。
    /// </summary>
    public class Line
    {
        /// <summary>
        /// y=ax+b中的a
        /// </summary>
        public double A { get; set; }
        /// <summary>
        /// y=ax+b中的b
        /// </summary>
        public double B { get; set; }
        /// <summary>
        /// x=c中的c
        /// </summary>
        public double C { get; set; }
        
        /// <summary>
        /// 構造函數（如果直線為y=ax+b形式，則C為Nan；如果直線為x=c形式，則A和B為Nan）
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="c"></param>
        public Line(double a, double b, double c)
        {
            if (!((double.IsNaN(a) && double.IsNaN(b) && !double.IsNaN(c)) || (!double.IsNaN(a) && !double.IsNaN(b) && double.IsNaN(c))))
                throw new ArgumentException("參數錯誤，無效的直線參數。");
            A = a;
            B = b;
            C = c;
        }

        /// <summary>
        /// 構造函數，由兩個點確定直線
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        public Line(PointF p1, PointF p2)
        {
            if (p1.X == p2.X)
            {
                A = double.NaN;
                B = double.NaN;
                C = p1.X;
            }
            else
            {
                A = 1d * (p1.Y - p2.Y) / (p1.X - p2.X);
                B = p1.Y - A * p1.X;
                C = double.NaN;
            }
        }

        /// <summary>
        /// 構造函數，由兩個點確定直線
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        public Line(Point p1, Point p2) :
            this(new PointF(p1.X, p1.Y), new PointF(p2.X, p2.Y))
        {
        }

        /// <summary>
        /// 生成一條隨機的直線
        /// </summary>
        /// <returns></returns>
        public static Line GetRandomLine()
        {
            Random random = new Random();
            int a = random.Next(-10, 10);
            int b = random.Next(-10, 10);
            return new Line(a, b, double.NaN);
        }

        /// <summary>
        /// 獲取點到直線的距離
        /// </summary>
        /// <param name="p">點</param>
        /// <returns>返回點到直線的距離；如果直線通過點，返回0。</returns>
        public double GetDistance(Point p)
        {
            return GetDistance(new PointF(p.X, p.Y));
        }

        /// <summary>
        /// 獲取點到直線的距離
        /// </summary>
        /// <param name="p">點</param>
        /// <returns>返回點到直線的距離；如果直線通過點，返回0。</returns>
        public double GetDistance(PointF p)
        {
            double d = 0d;
            if (double.IsNaN(C))
            {
                //y=ax+b相當於ax-y+b=0
                d = Math.Abs(1d * (A * p.X - p.Y + B) / Math.Sqrt(A * A + 1));
            }
            else
            {
                d = Math.Abs(C - p.X);
            }
            return d;
        }

        /// <summary>
        /// 根據x坐標，得到直線上點的y坐標 
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public double GetY(double x)
        {
            double y;
            if (double.IsNaN(C))
                y = A * x + B;
            else
                y = double.NaN;
            return y;
        }

        /// <summary>
        /// 返回直線方程
        /// </summary>
        /// <returns></returns>
        public override string ToString()
        {
            string formula = "";
            if (double.IsNaN(C))
                formula = string.Format("y={0}{1}", A != 0 ? string.Format("{0:F02}x", A) : "", B != 0 ? (B > 0 ? string.Format("+{0:F02}", B) : string.Format("{0:F02}", B)) : "");
            else
                formula = string.Format("x={0:F02}", C);
            return formula;
        }
    }
}

 

    3.3 檢測直線的過程

    （1）隨機從觀測點中選擇兩個點，得到通過該點的直線；

    （2）用（1）中的直線去測試其他觀測點，由點到直線的距離確定觀測點是否為局內點或者局外點；

    （3）如果局內點足夠多，並且局內點多於原有“最佳”直線的局內點，那么將這次迭代的直線設為“最佳”直線；

    （4）重復（1）~（3）步直到找到最佳直線。

    細心的您估計已經發現我省略了標准RANSAC檢測過程中重新估計模型的步驟，我是故意的，我覺得麻煩且沒什么用處，所以咔嚓了，O(∩_∩)O~。

/// <summary>
        /// 嘗試獲取直線
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void btnGetLine_Click(object sender, EventArgs e)
        {
            //用RANSAC方法獲取最佳直線
            points = GetSamplePoints();
            Line bestLine = null;           //最佳直線
            double bestInliersCount = 0;    //最佳模型的局內點數目
            Random random = new Random();
            for (int idx = 0; idx < nudIterCount.Value; idx++)
            {
                int idx1, idx2;
                GetRandomInliersPoints(random, out idx1, out idx2);
                int inliersCount = 2;
                Line line = new Line(points[idx1], points[idx2]);
                for (int i = 0; i < points.Count; i++)
                {
                    if (i != idx1 && i != idx2)
                    {
                        if (line.GetDistance(points[i]) <= (double)nudMinDistance.Value)
                            inliersCount++;
                    }
                }
                if (inliersCount >= nudMinPointCount.Value)
                {
                    if (inliersCount > bestInliersCount)
                    {
                        bestLine = line;
                        bestInliersCount = inliersCount;
                    }
                }
            }
            //顯示最佳直線
            if (bestLine != null)
            {
                lblFormula.Text = string.Format("方程:{0}\r\nA:{1}\r\nB:{2}\r\nC:{3}\r\n局內點數目:{4}",
                    bestLine.ToString(), bestLine.A, bestLine.B, bestLine.C, bestInliersCount);
                DrawLine(bestLine);
            }
            else
                lblFormula.Text = "沒有獲取到最佳直線。";
        }

 

四、檢測圓

    4.1 圓的相關知識

    （1）平面內不在同一直線上的三個點可以確定一個圓；

    （2）圓的數學表達形式為：(x-a)²+(y-b)²=r²

    其中，(a,b)為圓心，r為半徑。

    4.2 圓類

    圓類(Circle)封裝了跟圓有關的屬性及方法，列表如下：

    （1）屬性

    A——圓心的x坐標

    B——圓心的y坐標

     R——圓的半徑

    （2）構造函數

    public Circle(PointF p1, PointF p2, PointF p3)

    提供三個點p1，p2和p3，計算出圓的屬性A,B,R。

    （3）方法

    GetDistance——獲取點到圓（周）之間的距離，表示點接近或者遠離圓；

    ToString——獲取圓的方程式。

using System;
using System.Collections.Generic;
using System.Text;
using System.Drawing;
using System.Drawing.Drawing2D;

namespace Ransac
{
    /// <summary>
    /// 圓類：用(x-a)**2+(y-b)**2=r**2形式表示
    /// </summary>b
    public class Circle
    {
        /// <summary>
        /// 圓心的X坐標
        /// </summary>
        public double A { get; set; }
        /// <summary>
        /// 圓心的Y坐標
        /// </summary>
        public double B { get; set; }
        /// <summary>
        /// 半徑
        /// </summary>
        public double R { get; set; }

        /// <summary>
        /// 構造函數，提供圓心和半徑。
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="r"></param>
        public Circle(double a, double b, double r)
        {
            A = a;
            B = b;
            if (r < 0)
                throw new ArgumentOutOfRangeException("r", "圓的半徑必須大於0。");
            R = r;
        }

        /// <summary>
        /// 構造函數，提供三個點。
        /// 該算法來自csdn論壇，帖子地址是：http://bbs.csdn.net/topics/50383586，在此感謝5樓的privet。
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="p3"></param>
        public Circle(PointF p1, PointF p2, PointF p3)
        {
            float xMove = p1.X;
            float yMove = p1.Y;
            p1.X = 0;
            p1.Y = 0;
            p2.X = p2.X - xMove;
            p2.Y = p2.Y - yMove;
            p3.X = p3.X - xMove;
            p3.Y = p3.Y - yMove;
            float x1 = p2.X, y1 = p2.Y, x2 = p3.X, y2 = p3.Y;
            double m = 2.0 * (x1 * y2 - y1 * x2);
            if (m == 0)
                throw new ArgumentException("參數錯誤，提供的三點不能構成圓。");
            double x0 = (x1 * x1 * y2 - x2 * x2 * y1 + y1 * y2 * (y1 - y2)) / m;
            double y0 = (x1 * x2 * (x2 - x1) - y1 * y1 * x2 + x1 * y2 * y2) / m;
            R = Math.Sqrt(x0 * x0 + y0 * y0);
            A = x0 + xMove;
            B = y0 + yMove;
        }

        /// <summary>
        /// 構造函數，提供三個點。
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="p3"></param>
        public Circle(Point p1, Point p2, Point p3) :
            this(new PointF(p1.X, p1.Y), new PointF(p2.X, p2.Y), new PointF(p3.X, p3.Y))
        {
        }

        /// <summary>
        /// 獲取點到圓的距離（圓周，不是圓心）
        /// </summary>
        /// <param name="p"></param>
        /// <returns></returns>
        public double GetDistance(PointF p)
        {
            return Math.Abs(R - Math.Sqrt(Math.Pow(p.X - A, 2) + Math.Pow(p.Y - B, 2)));
        }

        /// <summary>
        /// 返回圓方程
        /// </summary>
        /// <returns></returns>
        public override string ToString()
        {
            return string.Format("{0}**2+{1}**2={2}",
                A == 0 ? "x" : (A > 0 ? string.Format("(x-{0:F02})", A) : string.Format("(x+{0:F02})", Math.Abs(A))),
                B == 0 ? "y" : (B > 0 ? string.Format("(y-{0:F02})", B) : string.Format("(y+{0:F02})", Math.Abs(B))),
                R == 0 ? "0" : string.Format("{0:F02}**2", Math.Abs(R)));
        }
    }
}

 

    3.3 檢測圓的過程

    （1）隨機從觀測點中選擇三個點，嘗試得到通過這三個點的圓；

    （2）用（1）中的圓去測試其他觀測點，由點到圓的距離確定觀測點是否為局內點或者局外點；

    （3）如果局內點足夠多，並且局內點多於原有“最佳”圓的局內點，那么將這次迭代的圓設為“最佳”圓；

    （4）重復（1）~（3）步直到找到最佳圓。

/// <summary>
        /// 嘗試獲取圓
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void btnGetCircle_Click(object sender, EventArgs e)
        {
            //用RANSAC方法獲取最佳直線
            points = GetSamplePoints();
            Circle bestCircle = null;       //最佳圓
            double bestInliersCount = 0;    //最佳模型的局內點數目
            Random random = new Random();
            for (int idx = 0; idx < nudIterCount.Value; idx++)
            {
                int idx1, idx2, idx3;
                GetRandomInliersPoints(random, out idx1, out idx2, out idx3);
                int inliersCount = 3;
                Circle circle;
                try
                {
                    circle = new Circle(points[idx1], points[idx2], points[idx3]);
                }
                catch
                {
                    continue;
                }
                for (int i = 0; i < points.Count; i++)
                {
                    if (i != idx1 && i != idx2 && i!=idx3)
                    {
                        if (circle.GetDistance(points[i]) <= (double)nudMinDistance.Value)
                            inliersCount++;
                    }
                }
                if (inliersCount >= nudMinPointCount.Value)
                {
                    if (inliersCount > bestInliersCount)
                    {
                        bestCircle = circle;
                        bestInliersCount = inliersCount;
                    }
                }
            }
            //顯示最佳圓
            if (bestCircle != null)
            {
                lblFormula.Text = string.Format("方程:{0}\r\nA:{1}\r\nB:{2}\r\nR:{3}\r\n局內點數目:{4}",
                    bestCircle.ToString(), bestCircle.A, bestCircle.B, bestCircle.R, bestInliersCount);
                DrawCircle(bestCircle);
            }
            else
                lblFormula.Text = "沒有獲取到最佳圓。";
        }

 

五、本文源代碼

    點擊這里下載本文源代碼。

    感謝您閱讀本文，希望對您有所幫助。