作者:王先荣
大约在两年前翻译了《随机抽样一致性算法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(-, );
int b = random.Next(-, );
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 + ));
}
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 != ? string.Format("{0:F02}x", A) : "", B != ? (B > ? string.Format("+{0:F02}", B) : string.Format("{0:F02}", B)) : "");
else
formula = string.Format("x={0:F02}", C);
return formula;
}
}
}
Line类
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 = ; //最佳模型的局内点数目
Random random = new Random();
for (int idx = ; idx < nudIterCount.Value; idx++)
{
int idx1, idx2;
GetRandomInliersPoints(random, out idx1, out idx2);
int inliersCount = ;
Line line = new Line(points[idx1], points[idx2]);
for (int i = ; 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)2+(y-b)2=r2
其中,(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 < )
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 = ;
p1.Y = ;
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 == )
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, ) + Math.Pow(p.Y - B, )));
} /// <summary>
/// 返回圆方程
/// </summary>
/// <returns></returns>
public override string ToString()
{
return string.Format("{0}**2+{1}**2={2}",
A == ? "x" : (A > ? string.Format("(x-{0:F02})", A) : string.Format("(x+{0:F02})", Math.Abs(A))),
B == ? "y" : (B > ? string.Format("(y-{0:F02})", B) : string.Format("(y+{0:F02})", Math.Abs(B))),
R == ? "" : string.Format("{0:F02}**2", Math.Abs(R)));
}
}
}
Circle类
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 = ; //最佳模型的局内点数目
Random random = new Random();
for (int idx = ; idx < nudIterCount.Value; idx++)
{
int idx1, idx2, idx3;
GetRandomInliersPoints(random, out idx1, out idx2, out idx3);
int inliersCount = ;
Circle circle;
try
{
circle = new Circle(points[idx1], points[idx2], points[idx3]);
}
catch
{
continue;
}
for (int i = ; 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 = "没有获取到最佳圆。";
}
获取圆
五、本文源代码
感谢您阅读本文,希望对您有所帮助。