一、什么是凹包
在计算几何中,凹包(Concave Hull)是一种几何形状,它包围一组点,形状的外部是凹凸不平的。与凸包(Convex Hull)相反,凹包允许形状沿着内部有一些凹陷。
具体来说,凹包是包围一组点的一个多边形,该多边形的边界是点的一个子集,而不是全部。凹包的外形可能会有凹陷的部分,而不是完全凸起。
凹包在计算机图形学、计算机辅助设计(CAD)、地理信息系统(GIS)等领域中都有应用。例如,在路径规划、物体识别、图像处理等任务中,需要找到一组点的凹包以描述它们的整体形状。
二、凹包的应用场景
凹包在计算机科学和工程中有许多应用场景,其中一些包括:
-
图形处理和计算机图形学: 在图形学中,凹包常用于表示和处理复杂的形状。例如,计算机游戏中的地形建模、虚拟现实环境中的场景建模等都可能涉及凹包的计算。
-
路径规划: 在机器人学、自动驾驶等领域,路径规划是一个关键问题。凹包可以用于表示障碍物的形状,帮助规划避障路径。
-
物体识别: 在计算机视觉中,物体识别可能涉及到对物体轮廓的建模,而凹包可以用于更准确地描述物体的形状。
-
地理信息系统(GIS): 在地理信息系统中,地图的边界可以通过凹包来表示,以更精确地捕捉地理区域的形状。
-
CAD(计算机辅助设计): 在设计领域,CAD软件中的形状建模和分析可能需要使用凹包。
-
数据可视化: 在数据可视化中,凹包可以用于突出数据点的集中趋势和形状。
这些应用场景只是冰山一角,凹包的概念和算法在各种计算和工程领域都有广泛的应用。
三、滚球算算法总体思路
3.1 滚球算法思路
滚球法(Rolling Ball Algorithm)是一种凹包算法,其基本思路是从凸包的角度来考虑。该算法的步骤可以概括为以下几点:
-
排序: 首先,对输入的点集按照一定的规则进行排序。通常可以选择按照点的坐标,先按照 y 坐标升序排序,然后在同样 y 坐标的情况下按照 x 坐标升序排序。
-
初始化: 初始化一个空的凹包,以及一个标记数组,用来标记哪些点已经被加入凹包中。初始时,凹包为空,标记数组中所有元素为未标记。
-
选取起始点: 从排序后的点集中选择一个起始点(通常是排序后的第一个点),将其加入凹包,并标记为已加入。
-
滚动球: 从起始点开始,采用滚球的方式寻找下一个凹包点。滚球的过程中,维护一个当前滚球半径。对于每个未加入凹包的点,判断是否在当前滚球半径内,如果是,则将其加入凹包,并将滚球半径设置为该点到凹包中已有点的最大距离。
-
迭代: 重复滚动球的过程,直到无法再加入新的点为止。这通常意味着所有的点都被包含在凹包中。
整体来说,滚球法是一种基于贪心策略的算法,通过不断扩展当前的凹包,直至包含所有的点。该算法相对简单,但需要注意的是滚球的半径的选择,这会影响算法的性能和结果。通常,滚球的半径可以根据数据的特点动态调整。
3.2 滚球法半径R对结果的影响
滚球算法中的半径R对算法的性能和输出结果具有一定的影响。推荐的半径R的计算通常需要根据具体的数据集和应用场景进行调整。以下是一些影响半径R选择的因素以及推荐半径的分析方法:
-
点的分布密度: 如果点的分布比较密集,适当增大半径R可能有助于将更多的点纳入同一个凹包,避免遗漏一些凹包的点。相反,如果点的分布比较稀疏,减小半径R可能更合适,以避免将不相关的点包含在同一个凹包中。
-
数据集的大小: 对于较大的数据集,可能需要更大的半径R,以确保覆盖整个数据集的关键几何特征。对于较小的数据集,较小的半径R可能足以捕捉到数据的主要形状。
-
应用需求: 不同的应用场景可能对凹包的定义有不同的要求。某些情况下,需要较大的凹包,而在其他情况下,可能需要更加精细的凹包。根据应用需求调整半径R可以满足不同的要求。
-
算法性能: 较小的半径R可能导致算法更为敏感,可能会在数据噪声较大或几何形状变化较为复杂的情况下产生不稳定的结果。较大的半径R可能使算法更加稳定,但可能会损失一些细节。
-
可视化效果: 在一些可视化应用中,推荐半径R的选择可能受到视觉效果的影响。合适的半径R应该能够清晰地展示数据的结构,而不至于使可视化失真或过于拥挤。
推荐半径R的计算通常需要在实际应用中进行试验和调整。可以通过尝试不同的半径值,观察凹包的质量和算法的性能,以及与特定应用场景的契合程度,来确定最合适的半径R。
要想顺利实现滚球法 ,需保证R大于最小半径,最小半径即至少能保证,以此为半径的圆能包含住2个点。
总的来说:R越大,则圆球所能滚到的下一个点就越远,凹包边界所包含的点集就越稀疏;R越小,则圆球所能滚到的下一个点就越近,凹包边界所包含的点集就越稠密。
左图为R=20000mm时的凹包边界,右图为R=10000mm时的凹包边界。
四、主要代码段介绍
4.1 首先定义一个二位坐标点类
每个点有 X 和 Y 坐标,Id 是点的编号,Distance 表示到其他点的距离。
public class Point
{
public double X { get; set; }
public double Y { get; set; }
public int Id { get; set; }
public double Distance { get; set; }
public Point(double x, double y)
{
X = x;
Y = y;
}
}
4.2 定义一个实现凹包算法的类
主要包含了初始化邻居列表和点之间距离的方法,计算默认半径的方法,以及使用滚球法获取凹包的方法。在这个类中,也有一些辅助方法,比如根据角度对邻居列表进行排序、检查圆内是否存在其他点等。
public class ConcaveHull
{
//...略
}
4.3 类的构造函数
它接收一个 List<Point>
参数,对这些点进行排序,然后初始化一些数据结构,包括 _signs
数组(用于标记点是否已被使用)、点之间的距离 _distances
数组,以及按距离排序的邻居列表 _neigbours
。
public ConcaveHull(List<Point> list)
{
this._points = list;
_points = _points.OrderBy(p => p.Y).ThenBy(p => p.X).ToList();
_signs = new bool[_points.Count];
for (int i = 0; i < _signs.Length; i++)
{
_signs[i] = false;
}
InitDistance();
InitNeighbours();
}
4.4 核心计算方法
使用滚球法获取凹包,如果不提供半径,会使用 CalDefaultRadius
计算默认半径。在循环中,通过 GetNextPoint
方法找到下一个凹包点,计算该点到当前点的圆心,并将该点加入结果列表。
public List<Point> ***pute(double radius = -1)
{
if (radius == -1)
{
radius = CalDefaultRadius();
}
List<Point> results = new List<Point>();
List<int>[] neighs = GetNeighbourList(2 * radius);
results.Add(_points[0]);
int i = 0, j = -1, pre = -1;
while (true)
{
j = GetNextPoint(pre, i, neighs[i], radius);
if (j == -1)
{
break;
}
Point p = CalCenterByPtsAndRadius(_points[i], _points[j], radius);
results.Add(_points[j]);
_signs[j] = true;
pre = i;
i = j;
}
return results;
}
CalDefaultRadius
方法计算默认的半径,即选择每个点的第二个邻居与其距离最大的点之间的距离最大值。
public double CalDefaultRadius()
{
double r = double.MinValue;
for (int i = 0; i < _points.Count; i++)
{
if (_distances[i, _neigbours[i][1]] > r)
{
r = _distances[i, _neigbours[i][1]];
}
}
return r;
}
五、完整代码
/// <summary>
/// 二维平面上的一个点
/// </summary>
/// </summary>
public class Point
{
/// <summary>
/// x坐标
/// </summary>
public double X { get; set; }
/// <summary>
/// y坐标
/// </summary>
public double Y { get; set; }
/// <summary>
/// 编号
/// </summary>
public int Id { get; set; }
/// <summary>
/// 到其他点的距离
/// </summary>
public double Distance { get; set; }
public Point(double x, double y)
{
X = x;
Y = y;
}
}
/// <summary>
/// 凹包类
/// </summary>
public class ConcaveHull
{
//点点之间距离列表
private double[,] _distances;
//邻居列表
private List<int>[] _neigbours;
private bool[] _signs;
private List<Point> _points;
public ConcaveHull(List<Point> list)
{
this._points = list;
_points = _points.OrderBy(p => p.Y).ThenBy(p => p.X).ToList();
_signs = new bool[_points.Count];
for (int i = 0; i < _signs.Length; i++)
{
_signs[i] = false;
}
InitDistance();
InitNeighbours();
}
/// <summary>
/// 计算默认的半径
/// </summary>
/// <returns></returns>
public double CalDefaultRadius()
{
double r = double.MinValue;
for (int i = 0; i < _points.Count; i++)
{
if (_distances[i, _neigbours[i][1]] > r)
{
r = _distances[i, _neigbours[i][1]];
}
}
return r;
}
/// <summary>
/// 使用滚球法获取凹包
/// 不输入半径时,用默认计算半径代替
/// </summary>
/// <param name="radius"></param>
/// <returns></returns>
public List<Point> ***pute(double radius=-1)
{
//计算默认半径
if (radius == -1)
{
radius = CalDefaultRadius();
}
List<Point> results = new List<Point>();
List<int>[] neighs = GetNeighbourList(2 * radius);
results.Add(_points[0]);
int i = 0, j = -1, pre = -1;
while (true)
{
j = GetNextPoint(pre, i, neighs[i], radius);
if (j == -1)
{
break;
}
Point p = CalCenterByPtsAndRadius(_points[i], _points[j], radius);
results.Add(_points[j]);
_signs[j] = true;
pre = i;
i = j;
}
return results;
}
/// <summary>
/// 初始化每个点的按距离排序的邻居列表
/// </summary>
private void InitNeighbours()
{
_neigbours = new List<int>[_points.Count];
for (int i = 0; i < _neigbours.Length; i++)
{
_neigbours[i] = SortNeighboursByDis(i);
}
}
/// <summary>
/// 初始化点之间的距离
/// </summary>
private void InitDistance()
{
_distances = new double[_points.Count, _points.Count];
for (int i = 0; i < _points.Count; i++)
{
for (int j = 0; j < _points.Count; j++)
{
_distances[i, j] = CalDistanceToPts(_points[i], _points[j]);
}
}
}
/// <summary>
/// 获取下一个凹包点
/// </summary>
/// <param name="pre"></param>
/// <param name="cur"></param>
/// <param name="list"></param>
/// <param name="radius"></param>
/// <returns></returns>
private int GetNextPoint(int pre, int cur, List<int> list, double radius)
{
SortNeighbours(list, pre, cur);
for (int j = 0; j < list.Count; j++)
{
if (_signs[list[j]])
{
continue;
}
int adjIndex = list[j];
Point xianp = _points[adjIndex];
Point rightCirleCenter = CalCenterByPtsAndRadius(_points[cur], xianp, radius);
if (!IsPointsInCircle(list, rightCirleCenter, radius, adjIndex))
{
return list[j];
}
}
return -1;
}
/// <summary>
/// 根据角度对邻居列表进行排序
/// </summary>
/// <param name="list"></param>
/// <param name="pre"></param>
/// <param name="cur"></param>
private void SortNeighbours(List<int> list, int pre, int cur)
{
Point origin = _points[cur];
Point df;
if (pre != -1)
{
df = new Point(_points[pre].X - origin.X, _points[pre].Y - origin.Y);
}
else
{
df = new Point(1, 0);
}
int temp = 0;
for (int i = list.Count; i > 0; i--)
{
for (int j = 0; j < i - 1; j++)
{
if (***pareAngle(_points[list[j]], _points[list[j + 1]], origin, df))
{
temp = list[j];
list[j] = list[j + 1];
list[j + 1] = temp;
}
}
}
}
/// <summary>
/// 检查圆内是否存在其他点
/// </summary>
/// <param name="roundPts"></param>
/// <param name="center"></param>
/// <param name="radius"></param>
/// <param name="roundId"></param>
/// <returns></returns>
private bool IsPointsInCircle(List<int> roundPts, Point center, double radius, int roundId)
{
for (int k = 0; k < roundPts.Count; k++)
{
if (roundPts[k] != roundId)
{
int index2 = roundPts[k];
if (IsPtInCircle(_points[index2], center, radius))
{
return true;
}
}
}
return false;
}
/// <summary>
/// 根据两点及半径计算圆心
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="r"></param>
/// <returns></returns>
private static Point CalCenterByPtsAndRadius(Point a, Point b, double r)
{
double dx = b.X - a.X;
double dy = b.Y - a.Y;
double cx = 0.5 * (b.X + a.X);
double cy = 0.5 * (b.Y + a.Y);
if (r * r / (dx * dx + dy * dy) - 0.25 < 0)
{
return new Point(-1, -1);
}
double sqrt = Math.Sqrt(r * r / (dx * dx + dy * dy) - 0.25);
return new Point(cx - dy * sqrt, cy + dx * sqrt);
}
/// <summary>
/// 检查点是否在圆内
/// </summary>
/// <param name="p"></param>
/// <param name="center"></param>
/// <param name="r"></param>
/// <returns></returns>
private static bool IsPtInCircle(Point p, Point center, double r)
{
double dis2 = (p.X - center.X) * (p.X - center.X) + (p.Y - center.Y) * (p.Y - center.Y);
return dis2 < r * r;
}
/// <summary>
/// 获取半径范围内的邻居列表
/// </summary>
/// <param name="radius"></param>
/// <returns></returns>
private List<int>[] GetNeighbourList(double radius)
{
List<int>[] adjs = new List<int>[_points.Count];
for (int i = 0; i < _points.Count; i++)
{
adjs[i] = new List<int>();
}
for (int i = 0; i < _points.Count; i++)
{
for (int j = 0; j < _points.Count; j++)
{
if (i < j && _distances[i, j] < radius)
{
adjs[i].Add(j);
adjs[j].Add(i);
}
}
}
return adjs;
}
/// <summary>
/// 按距离排序的邻居列表
/// </summary>
/// <param name="index"></param>
/// <returns></returns>
private List<int> SortNeighboursByDis(int index)
{
List<Point> pts = new List<Point>();
for (int i = 0; i < _points.Count; i++)
{
Point pt = _points[i];
pt.Id = i;
pt.Distance = _distances[index,i];
pts.Add(pt);
}
pts= pts.OrderBy(p=>p.Distance).ToList();
List<int> adj = new List<int>();
for (int i = 1; i < pts.Count; i++)
{
adj.Add(pts[i].Id);
}
return adj;
}
/// <summary>
/// 比较两个点相对于参考点的角度
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <param name="origin"></param>
/// <param name="reference"></param>
/// <returns></returns>
private bool ***pareAngle(Point a, Point b, Point origin, Point reference)
{
Point da = new Point(a.X - origin.X, a.Y - origin.Y);
Point db = new Point(b.X - origin.X, b.Y - origin.Y);
//b相对于参考向量的叉积
double detb = CalCrossProduct(reference, db);
// 如果 b 的叉积为零且 b 与参考向量的夹角大于等于零度,则返回 false
if (detb == 0 && db.X * reference.X + db.Y * reference.Y >= 0)
{
return false;
}
//a 相对于参考向量的叉积
double deta = CalCrossProduct(reference, da);
//如果 a 的叉积为零且 a 与参考向量的夹角大于等于零度,则返回 true
if (deta == 0 && da.X * reference.X + da.Y * reference.Y >= 0)
{
return true;
}
// 如果 a 和 b 在参考向量的同一侧,则比较它们之间的叉积大小
if (deta * detb >= 0)
{
//如果叉积大于零,返回 true;否则返回 false
return CalCrossProduct(da, db) > 0;
}
//向量小于零度实际上是很大的,接近 2pi
return deta > 0;
}
/// <summary>
/// 计算两点之间的距离
/// </summary>
/// <param name="p1"></param>
/// <param name="p2"></param>
/// <returns></returns>
private static double CalDistanceToPts(Point p1, Point p2)
{
return Math.Sqrt((p1.X - p2.X) * (p1.X - p2.X) + (p1.Y - p2.Y) * (p1.Y - p2.Y));
}
/// <summary>
/// 计算两个向量的叉积
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <returns></returns>
private static double CalCrossProduct(Point a, Point b)
{
return a.X * b.Y - a.Y * b.X;
}
}