二维点集求外轮廓Java_从二维点集重建平面形状-浅议凹包算法
问题背景
近期遇到一个计算几何问题,需要从点集中重建一个合理的几何形状。这个问题既有二维的也有三维的,二维的情况相对简单一点,即给出平面区域的一系列散点,求出一定程度上反映这些散点轮廓的平面多边形,给出边的连接方式即可。如从下图的左图散点重建为右图的形状:
二维平面散点
平面多边形
不过这里有一些细节需要注意,必须明确这一系列点的含义,有时给出的点集是表征图形边界的,如左图的情况;有时则是表征图形所覆盖的范围,即在图形的内部也有一定的点分布。这两种不同的含义是和具体的应用场景相关的。
有的场景是求散点描绘的轮廓
有的场景是求散点分布区域的边界轮廓
凹包与Alpha shape
有不少计算几何领域的资料探讨过这一类问题的解决方案,笔者也曾经阅读过部分相关的文献,也在网络上进行过一些搜索。本文不能一一去介绍所有的方法,所以就介绍几种简单直观的方法来解决二维点集重建平面多边形的问题。
首先很多学计算机的人都知道凸包,凸包求取的是覆盖所有点的凸多边形,由于限定死了一定要凸多边形,所以凸包不是本文所讨论问题的解决方案。本文所求的形状肯定不能否定存在凹陷的部分,因而可以联想到这个问题的解是否是求一种“凹包”,或者说一种广义的参数化的凸包。事实上的确有凹包的概念,英文叫做concave hull(凸包叫convex hull)。不过于凸包的情况不同,凹包没有特别明确的定义,给定一个较不规则的点集,可以有各种各样的凹包,但凸包是确定的,如下图所示。
凸包
一种可能的凹包连接方式
这样连也是凹包
进一步查找相关的资料可以发现一个叫做alpha shape的概念。这个概念最早出现于论文"On the shape of a set of points in the plane"中。alpha shape有较为正式的定义,若要简单点解释它,alpha shape其实就是在凸包基础上多了一个可以设定的参数α,因为这个α,在alpha shape重建形状的过程中可能就不会像凸包一样连接相距过远的顶点。参数α若趋于无穷大,则这个alpha shape会无无限接近凸包;而α取小了,则alpha shape会倾向于在一定位置凹陷进去,以更加贴合点集的外形。选择合理的α就能够控制生成的形状让其更加符合点集所描绘的形状。所以给出一个α值,就有点类似于给出一个长度的限制,例如让多边形的任何一边长度不超过R。
为了更好的展示和测试程序,笔者特地开发了一个基于WPF技术的小软件ConcaveGenerator用来展示算法效果,下文展示的很多图都是对这个软件截图的。如有兴趣可以在这里下载。
自制ConcaveGenerator小软件截图
第一种思路--基于Delaunay三角化
有一定计算几何基础的人一定熟悉Delaunay三角化,通过这一个过程可以形成Delaunay三角网,假如我们为想要求取形状的点集上使用Delaunay三角化算法,可以得到如下的三角网。
一个点集的Delaunay三角网
另一个例子
一看到这样的三角网,就不难想到一个好点子:因为我们想要求取的形状就是Delaunay三角网的子集,所以我们只需要从Delaunay三角网的最外层边开始外不断的删去超过长度限制的边,当这个过程结束时,我们就能够得到一个大致符合我们预期的形状。
点集的Delaunay三角网
删掉边上太长的边就能形成预期的形状
所以总结这个思路,输入为点集S和长度限制R的求取凹包的边列表算法的过程如下:
为点集S求取Delaunay三角网M,三角网用标准Mesh形式表示。
为M初始化所有Edge对象,并求取Edge的长度以及邻接三角形集合。其中邻接2个三角形的边为内部边,1个三角形的为边界边,0个三角形的为计算过程中会退化的边。
将所有长度大于R的边界边加入队列,并当队列非空时循环下列过程:从队列中取出一条边E,得到E的唯一邻接三角形T。
找到T中另外两个边E1,E2将他们的邻接三角形集合删去T。
将E1,E2中新形成的长度大于R的边界边加入队列。
将E置无效标记,若E1,E2有退化的,也置无效标记。
收集所有有效的边界边,形成边列表,输出。
以下ConcaveGenerator是在这部分的代码,其中Delaunay三角化的实现有比较多的代码能参考,这里是使用了Stackoverflow上一个snippet:
public classDelaunayTriangulation
{public const int MaxVertices = 500;public const int MaxTriangles = 1000;public dVertex[] Vertex = newdVertex[MaxVertices];public dTriangle[] Triangle = newdTriangle[MaxTriangles];private bool InCircle(long xp, long yp, long x1, long y1, long x2, long y2, long x3, long y3, double xc, double yc, doubler)
{doubleeps;doublem1;doublem2;doublemx1;doublemx2;doublemy1;doublemy2;doubledx;doubledy;doublersqr;doubledrsqr;
eps= 0.000000001;if (Math.Abs(y1 - y2) < eps && Math.Abs(y2 - y3)
{return false;
}if (Math.Abs(y2 - y1)
{
m2= (-(Convert.ToDouble(x3) - Convert.ToDouble(x2)) / (Convert.ToDouble(y3) -Convert.ToDouble(y2)));
mx2= Convert.ToDouble((x2 + x3) / 2.0);
my2= Convert.ToDouble((y2 + y3) / 2.0);
xc= Convert.ToDouble((x2 + x1) / 2.0);
yc= Convert.ToDouble(m2 * (xc - mx2) +my2);
}else if (Math.Abs(y3 - y2)
{
m1= (-(Convert.ToDouble(x2) - Convert.ToDouble(x1)) / (Convert.ToDouble(y2) -Convert.ToDouble(y1)));
mx1= Convert.ToDouble((x1 + x2) / 2.0);
my1= Convert.ToDouble((y1 + y2) / 2.0);
xc= Convert.ToDouble((x3 + x2) / 2.0);
yc= Convert.ToDouble(m1 * (xc - mx1) +my1);
}else{
m1= (-(Convert.ToDouble(x2) - Convert.ToDouble(x1)) / (Convert.ToDouble(y2) -Convert.ToDouble(y1)));
m2= (-(Convert.ToDouble(x3) - Convert.ToDouble(x2)) / (Convert.ToDouble(y3) -Convert.ToDouble(y2)));
mx1= Convert.ToDouble((x1 + x2) / 2.0);
mx2= Convert.ToDouble((x2 + x3) / 2.0);
my1= Convert.ToDouble((y1 + y2) / 2.0);
my2= Convert.ToDouble((y2 + y3) / 2.0);
xc= Convert.ToDouble((m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 -m2));
yc= Convert.ToDouble(m1 * (xc - mx1) +my1);
}
dx= (Convert.ToDouble(x2) -Convert.ToDouble(xc));
dy= (Convert.ToDouble(y2) -Convert.ToDouble(yc));
rsqr= Convert.ToDouble(dx * dx + dy *dy);
r=Convert.ToDouble(Math.Sqrt(rsqr));
dx= Convert.ToDouble(xp -xc);
dy= Convert.ToDouble(yp -yc);
drsqr= Convert.ToDouble(dx * dx + dy *dy);if (drsqr <=rsqr)
{return true;
}return false;
}private int WhichSide(long xp, long yp, long x1, long y1, long x2, longy2)
{doubleequation;
equation= ((Convert.ToDouble(yp) - Convert.ToDouble(y1)) * (Convert.ToDouble(x2) - Convert.ToDouble(x1))) - ((Convert.ToDouble(y2) - Convert.ToDouble(y1)) * (Convert.ToDouble(xp) -Convert.ToDouble(x1)));if (equation > 0)
{return -1;//WhichSide = -1;
}else if (equation == 0)
{return 0;
}else{return 1;
}
}public int Triangulate(intnvert)
{bool[] Complete = new bool[MaxTriangles];long[,] Edges = new long[3, MaxTriangles * 3 + 1];longNedge;longxmin;longxmax;longymin;longymax;longxmid;longymid;doubledx;doubledy;doubledmax;inti;intj;intk;intntri;double xc = 0.0;double yc = 0.0;double r = 0.0;boolinc;
xmin= Vertex[1].x;
ymin= Vertex[1].y;
xmax=xmin;
ymax=ymin;for (i = 2; i <= nvert; i++)
{if (Vertex[i].x
{
xmin=Vertex[i].x;
}if (Vertex[i].x >xmax)
{
xmax=Vertex[i].x;
}if (Vertex[i].y
{
ymin=Vertex[i].y;
}if (Vertex[i].y >ymax)
{
ymax=Vertex[i].y;
}
}
dx= Convert.ToDouble(xmax) -Convert.ToDouble(xmin);
dy= Convert.ToDouble(ymax) -Convert.ToDouble(ymin);if (dx >dy)
{
dmax=dx;
}else{
dmax=dy;
}
xmid= (xmax + xmin) / 2;
ymid= (ymax + ymin) / 2;
Vertex[nvert+ 1].x = Convert.ToInt64(xmid - 2 *dmax);
Vertex[nvert+ 1].y = Convert.ToInt64(ymid -dmax);
Vertex[nvert+ 2].x =xmid;
Vertex[nvert+ 2].y = Convert.ToInt64(ymid + 2 *dmax);
Vertex[nvert+ 3].x = Convert.ToInt64(xmid + 2 *dmax);
Vertex[nvert+ 3].y = Convert.ToInt64(ymid -dmax);
Triangle[1].vv0 = nvert + 1;
Triangle[1].vv1 = nvert + 2;
Triangle[1].vv2 = nvert + 3;
Complete[1] = false;
ntri= 1;for (i = 1; i <= nvert; i++)
{
Nedge= 0;
j= 0;do{
j= j + 1;if (Complete[j] != true)
{
inc=InCircle(Vertex[i].x, Vertex[i].y, Vertex[Triangle[j].vv0].x, Vertex[Triangle[j].vv0].y, Vertex[Triangle[j].vv1].x, Vertex[Triangle[j].vv1].y, Vertex[Triangle[j].vv2].x, Vertex[Triangle[j].vv2].y, xc, yc, r);if(inc)
{
Edges[1, Nedge + 1] =Triangle[j].vv0;
Edges[2, Nedge + 1] =Triangle[j].vv1;
Edges[1, Nedge + 2] =Triangle[j].vv1;
Edges[2, Nedge + 2] =Triangle[j].vv2;
Edges[1, Nedge + 3] =Triangle[j].vv2;
Edges[2, Nedge + 3] =Triangle[j].vv0;
Nedge= Nedge + 3;
Triangle[j].vv0=Triangle[ntri].vv0;
Triangle[j].vv1=Triangle[ntri].vv1;
Triangle[j].vv2=Triangle[ntri].vv2;
Complete[j]=Complete[ntri];
j= j - 1;
ntri= ntri - 1;
}
}
}while (j
{if (Edges[1, j] != 0 && Edges[2, j] != 0)
{for (k = j + 1; k <= Nedge; k++)
{if (Edges[1, k] != 0 && Edges[2, k] != 0)
{if (Edges[1, j] == Edges[2, k])
{if (Edges[2, j] == Edges[1, k])
{
Edges[1, j] = 0;
Edges[2, j] = 0;
Edges[1, k] = 0;
Edges[2, k] = 0;
}
}
}
}
}
}for (j = 1; j <= Nedge; j++)
{if (Edges[1, j] != 0 && Edges[2, j] != 0)
{
ntri= ntri + 1;
Triangle[ntri].vv0= Edges[1, j];
Triangle[ntri].vv1= Edges[2, j];
Triangle[ntri].vv2=i;
Complete[ntri]= false;
}
}
}
i= 0;do{
i= i + 1;if (Triangle[i].vv0 > nvert || Triangle[i].vv1 > nvert || Triangle[i].vv2 >nvert)
{
Triangle[i].vv0=Triangle[ntri].vv0;
Triangle[i].vv1=Triangle[ntri].vv1;
Triangle[i].vv2=Triangle[ntri].vv2;
i= i - 1;
ntri= ntri - 1;
}
}while (i
}public static double Diameter(double Ax, double Ay, double Bx, double By, double Cx, doubleCy)
{doublex, y;double a =Ax;double b =Bx;double c =Cx;double m =Ay;double n =By;double k =Cy;double A = a * b * b + a * n * n + a * a * c - b * b * c + m * m * c - n * n * c - a * c * c - a * k * k - a * a * b + b * c * c - m * m * b + b * k *k;double B = a * n + m * c + k * b - n * c - a * k - b *m;
y= A / B / 2;double AA = b * b * m + m * n * n + a * a * k - b * b * k + m * m * k - n * n * k - c * c * m - m * k * k - a * a * n + c * c * n - m * m * n + k * k *n;double BB = b * m + a * k + c * n - b * k - c * m - a *n;
x= AA / BB / 2;return Math.Sqrt((Ax - x) * (Ax - x) + (Ay - y) * (Ay -y));
}public static double MaxEdge(double Ax, double Ay, double Bx, double By, double Cx, doubleCy)
{double len1 = Math.Sqrt((Ax - Bx) * (Ax - Bx) + (Ay - By) * (Ay -By));double len2 = Math.Sqrt((Cx - Bx) * (Cx - Bx) + (Cy - By) * (Cy -By));double len3 = Math.Sqrt((Ax - Cx) * (Ax - Cx) + (Ay - Cy) * (Ay -Cy));double len = len1 > len2 ?len1 : len2;return len > len3 ?len : len3;
}
}
以下是基于Delaunay的求凹包算法部分核心函数:
public structdVertex
{public longx;public longy;public longz;
}public structdTriangle
{public longvv0;public longvv1;public longvv2;
}public structOrTriangle
{publicPoint2d p0;publicPoint2d p1;publicPoint2d p2;
}public structTriangle
{public intP0Index;public intP1Index;public intP2Index;public intIndex;public Triangle(int p0index, int p1index, intp2index)
{this.P0Index =p0index;this.P1Index =p1index;this.P2Index =p2index;this.Index = -1;
}public Triangle(int p0index, int p1index, int p2index, intindex)
{this.P0Index =p0index;this.P1Index =p1index;this.P2Index =p2index;this.Index =index;
}
}public structEdgeInfo
{public intP0Index;public intP1Index;public ListAdjTriangle;public boolFlag;public doubleLength;public intGetEdgeType()
{returnAdjTriangle.Count;
}public boolIsValid()
{return P0Index != -1;
}public EdgeInfo(intd)
{
P0Index= -1;
P1Index= -1;
Flag= false;
AdjTriangle= new List();
Length= -1;
}
}public classDelaunayMesh2d
{public ListPoints;public ListFaces;publicEdgeInfo[,] Edges;publicDelaunayMesh2d()
{
Points= new List();
Faces= new List();
}public intAddVertex(Point2d p)
{
Points.Add(p);return Points.Count - 1;
}public intAddFace(Triangle t)
{
Faces.Add(t);return Faces.Count - 1;
}public voidInitEdgesInfo()
{
Edges= newEdgeInfo[Points.Count, Points.Count];for (int i = 0; i < Points.Count; i++)
{for (int j = 0; j < Points.Count; j++)
{
Edges[i, j]= new EdgeInfo(0);
}
}for (int i = 0; i < Faces.Count; i++)
{
Triangle t=Faces[i];
SetEdge(t, i);
}
}private void SetEdge(Triangle t, inti)
{if (t.P0Index
{
Edges[t.P0Index, t.P1Index].P0Index=t.P0Index;
Edges[t.P0Index, t.P1Index].P1Index=t.P1Index;
Edges[t.P0Index, t.P1Index].AdjTriangle.Add(i);
Edges[t.P0Index, t.P1Index].Length=BallConcave.GetDistance(Points[t.P0Index], Points[t.P1Index]);
}else{
Edges[t.P1Index, t.P0Index].P0Index=t.P1Index;
Edges[t.P1Index, t.P0Index].P1Index=t.P0Index;
Edges[t.P1Index, t.P0Index].AdjTriangle.Add(i);
Edges[t.P1Index, t.P0Index].Length=BallConcave.GetDistance(Points[t.P0Index], Points[t.P1Index]);
}if (t.P1Index
{
Edges[t.P1Index, t.P2Index].P0Index=t.P1Index;
Edges[t.P1Index, t.P2Index].P1Index=t.P2Index;
Edges[t.P1Index, t.P2Index].AdjTriangle.Add(i);
Edges[t.P1Index, t.P2Index].Length=BallConcave.GetDistance(Points[t.P1Index], Points[t.P2Index]);
}else{
Edges[t.P2Index, t.P1Index].P0Index=t.P2Index;
Edges[t.P2Index, t.P1Index].P1Index=t.P1Index;
Edges[t.P2Index, t.P1Index].AdjTriangle.Add(i);
Edges[t.P2Index, t.P1Index].Length=BallConcave.GetDistance(Points[t.P1Index], Points[t.P2Index]);
}if (t.P0Index
{
Edges[t.P0Index, t.P2Index].P0Index=t.P0Index;
Edges[t.P0Index, t.P2Index].P1Index=t.P2Index;
Edges[t.P0Index, t.P2Index].AdjTriangle.Add(i);
Edges[t.P0Index, t.P2Index].Length=BallConcave.GetDistance(Points[t.P0Index], Points[t.P2Index]);
}else{
Edges[t.P2Index, t.P0Index].P0Index=t.P2Index;
Edges[t.P2Index, t.P0Index].P1Index=t.P0Index;
Edges[t.P2Index, t.P0Index].AdjTriangle.Add(i);
Edges[t.P2Index, t.P0Index].Length=BallConcave.GetDistance(Points[t.P0Index], Points[t.P2Index]);
}
}public void ExecuteEdgeDecimation(doublelength)
{
Queue queue = new Queue();for (int i = 0; i < Points.Count; i++)
{for (int j = 0; j < Points.Count; j++)
{if (i < j &&Edges[i, j].IsValid())
{if (Edges[i, j].GetEdgeType() == 0)
{throw newException();
}if (Edges[i, j].Length > length && Edges[i, j].GetEdgeType() == 1)
{
queue.Enqueue(Edges[i, j]);
}
}
}
}
EdgeInfo[] opp1Temp= new EdgeInfo[2];while (queue.Count != 0)
{
EdgeInfo info=queue.Dequeue();if (info.AdjTriangle.Count != 1)throw newException();int tindex = info.AdjTriangle[0];
Triangle t=Faces[tindex];
InitOppEdge(opp1Temp, t, info);
SetInvalid(info.P0Index, info.P1Index);for (int i = 0; i < 2; i++)
{
EdgeInfo e=opp1Temp[i];
e.AdjTriangle.Remove(tindex);if (e.GetEdgeType() == 0)
{
SetInvalid(e.P0Index, e.P1Index);
}else if (e.GetEdgeType() == 1 && e.Length >length)
{
queue.Enqueue(e);
}
}
}
}public ListGetBoundaryEdges()
{
List list = new List();for (int i = 0; i < Points.Count; i++)
{for (int j = 0; j < Points.Count; j++)
{if (i
{if (Edges[i, j].GetEdgeType() == 1)
{
list.Add(Edges[i, j]);
}
}
}
}returnlist;
}private void SetInvalid(int i, intj)
{
Edges[i, j].AdjTriangle.Clear();
Edges[i, j].Flag= true;
Edges[i, j].P0Index= -1;
Edges[i, j].P1Index= -1;
}private voidInitOppEdge(EdgeInfo[] opp1Temp, Triangle t, EdgeInfo info)
{int vindex = t.P0Index + t.P1Index + t.P2Index - info.P0Index -info.P1Index;if (vindex
{
opp1Temp[0] =Edges[vindex, info.P0Index];
}else{
opp1Temp[0] =Edges[info.P0Index, vindex];
}if (vindex
{
opp1Temp[1] =Edges[vindex, info.P1Index];
}else{
opp1Temp[1] =Edges[info.P1Index, vindex];
}
}
}
下图是用基于Delauney的方法生成的凹包,看上去大致符合我们的预期:
生成结果
显示三角网
第二种思路--滚边法(Edge Pivoting)
这是笔者最初想到的一个从求凸包的Graham Scan算法衍生出来的一个方法。求凸包的Graham Scan算法先找到一个Y最低的点作为起始点,然后使用叉积角度判断的方法去判断点的走向,最后在栈内留下了凸包的点序列。具体的算法讲解与代码,网上一搜各种有,这里就不详细表述。本文要介绍的方法也是和Graham Scan法从同一个出发点出发,通过扫描角度来确定下一个点。具体算法流程从下面的图文说明可以大致看出来:
首先要实现这个算法,需要对随机的一个点查询距离其几何距离在R内的所有点,即求所谓的R邻域。这个R邻域算法是KNN(K-nearest neighbors)算法的一个变形,可以在小于O(n2)的复杂度下求取R领域,本文为了简单起见没有实现这个基于K-d Tree的算法,感兴趣的读者可以查阅相关资料了解这个算法。本文涉及的应用场景因为点数目不大,所以使用的方法没有过多考虑效率,实现R邻域的方式是先建立一个n*n的二维数组存储所有点两两之间的距离,然后遍历一次二维数组,为所有的点建立一个R邻域列表,该列表中存储了对应点R邻域的索引值,这个列表很有用,下面要介绍的滚球法也利用了这个列表。
实际上不难理解,假设AB为凹包的一个边,那么其下一个点C必然是在B的R邻域点中选择。我们实现这个算法的关键,就是为AB边找一个合适的C点。这里笔者设想的寻找C的方法使用如下原则:
假如B的R领域除了A就只有一个点,那么那个点就是C。
以B为圆心从向量BA出发转圈扫描,遇到的第一个点为C。
这里涉及到一个小算法,即所谓的极坐标方向排序问题,这个问题的描述是:给定一个原点P和一个初始方向n,如何为平面上的点集S排序,使得点集中的点P1,P2...PN与P的连接是按从n开始的逆时针排列的。这个问题搜一搜stackoverflow即得到一个很好的解答,这个链接里一个人给出一个用于比较的函数,一旦有了比较函数,排序也不成问题,这个比较函数在后面的方法中会出现。其具体的比较原理如果对向量的点积叉积有所了解也不难理解。这里不妨提一下点集叉积的结果符号的几何含义:
向量a与b的点集结果为一个实数,计算方式是ax*bx+ay*by,满足交换率,为正值代表ab夹角小于90度,为锐角,负值代表夹角大于90度(谈夹角的话是指0-180度范围),为钝角。
向量a与b的叉积结果为一个向量,其数值的计算方法是ax*bx-ay*by,不满足交换率,为正值代表向量b处于向量a的逆时针方向(即a逆时针转一个小于180的角能转到b),负值代表向量b处于向量a的顺时针方向(即a顺时针转一个小于180的角能转到b,非要逆时针转则必然超过180度)。
那么找C点的第二条实现的方式就类似于对一个数组找最小值那样,通过比较找到最小的角度,这个有最小角偏向的就是C点。不过遗憾的是这个思路其实是问题的,测试表明对一些点集采用这个方法会有BUG出现。例如当C点出现在BA向量小于90度的方向时,形成的BC边和AB边具有大致相反的方向,会导致下一步的寻点出现逆向,故这个思路不算是一个成功的思路,不过失败是成功之母,它却引出另一个滚球法的思路,相比这个思路具有更好的鲁棒性。
第三种思路--滚球法(Ball Pivoting)
对于任何点集,我们把这些点想象为钉在平面上的钉子。假如拿一个半径大于一定值的球去从边界接近这个钉群,我们可以用这个球在这些钉子群的边界滚一圈,每次滚动,球都能接触到两个钉子而被卡住。
这个思路要求一个合法的R,R太小就没法形成一个闭合的图形。由于我们讨论问题的初衷是要形成一个合适的多边形而不是0个或多个,这样对R的选择就应该有一个限制,太小的R必然出不了结果,这里姑且假设给的R值是合适的。此过程若形成一个多边形,则多变形的最长的边一定小于球的直径。所以算法输入参数为R,意味着拿一个半径为R/2的圆去滚。如下图显示了这个滚球的过程:
我们不难发现一个性质,对于任何一次翻滚,假设从弦DE翻转到新弦EF,圆不可能扫过任何的其他点,因为假如扫到其他点,必然会被这个点所卡住,那么新弦就不可能是EF了。这样我们只需对极坐标排序后的E点的R邻域依次寻找符合不覆盖其他点的圆即可。
圆在翻滚时候不能扫到其他点
对E的R邻域测试圆覆盖情况寻找合适的F
这个算法的执行过程和滚边法相似,算法结构都是先找起始点,然后循环寻找下一个点,核心问题也是从边DE线段出发找新线段EF,只是不再使用边去扫,而是用圆去扫。这里先给出这个算法的大致步骤:
先像求凸包那样求出一个Y值最小(Y值相同则取X最大)的点,作为初始点A,此点一定在凹包上。
从此点出发,(1,0)为基准向量,先找一个小于R的边作为初始边,这时这个点即为B,此时一个半径为R/2的圆就卡在了AB上,此时第一个弦AB就找到了。
循环寻找接下来的弦,假如上一条弦为DE,则下一条弦必然从E开始,连接到E的一个R领域内的点F。寻找F可以使用如下的原则:先对E的R领域的点,以E为中心ED向量为基准进行极坐标方向排序,之后依次为R领域点F0~FN建立以EFi为弦的圆,然后检查其中是否包含其他F0~FN的点,若不存在,则EFi即为新弦,则跳出循环。
依次找到所有的弦,直到找不到新弦或遇到以前已经作为弦的点为止。
一旦R值给的比较好,这个过程一定能给出一个闭合的图形。下图为几张用ConcaveGenerator生成的图片示例,其中检查参数按钮会自动给一个能让所有点都有至少2个领域的R值。
一个例子,可以看出一些点没有被滚过
上述过程其实感觉还有很大的优化余地,不过在点数不是很多的场合还是能姑且一用的。
在ConcaveGenerator中这部分的代码为下面的C#类,这个类既实现了滚边法也实现了滚球法:
public classBallConcave
{publicMainWindow main;struct Point2dInfo : IComparable{publicPoint2d Point;public intIndex;public doubleDistanceTo;public Point2dInfo(Point2d p,int i,doubledis)
{this.Point =p;this.Index =i;this.DistanceTo =dis;
}public intCompareTo(Point2dInfo other)
{returnDistanceTo.CompareTo(other.DistanceTo);
}public override stringToString()
{return Point+","+Index+","+DistanceTo;
}
}public BallConcave(Listlist)
{this.points =list;
points.Sort();
flags= new bool[points.Count];for (int i = 0; i < flags.Length; i++)
flags[i]= false;
InitDistanceMap();
InitNearestList();
}private bool[] flags;private Listpoints;private double[,] distanceMap;private List[] rNeigbourList;private voidInitNearestList()
{
rNeigbourList= new List[points.Count];for (int i = 0; i < rNeigbourList.Length; i++)
{
rNeigbourList[i]=GetSortedNeighbours(i);
}
}private voidInitDistanceMap()
{
distanceMap= new double[points.Count, points.Count];for (int i = 0; i < points.Count; i++)
{for (int j = 0; j < points.Count; j++)
{
distanceMap[i, j]=GetDistance(points[i], points[j]);
}
}
}public doubleGetRecomandedR()
{double r=double.MinValue;for (int i = 0; i < points.Count; i++)
{if (distanceMap[i, rNeigbourList[i][1]] >r)
r= distanceMap[i, rNeigbourList[i][1]];
}returnr;
}public doubleGetMinEdgeLength()
{double min = double.MaxValue;for (int i = 0; i < points.Count; i++)
{for (int j = 0; j < points.Count; j++)
{if (i
{if (distanceMap[i, j]
min=distanceMap[i, j];
}
}
}returnmin;
}public List GetConcave_Ball(doubleradius)
{
List ret = new List();
List[] adjs = GetInRNeighbourList(2*radius);
ret.Add(points[0]);//flags[0] = true;
int i = 0, j = -1, prev = -1;while (true)
{
j=GetNextPoint_BallPivoting(prev, i, adjs[i],radius);if (j == -1)break;
Point2d p=BallConcave.GetCircleCenter(points[i], points[j], radius);
ret.Add(points[j]);
flags[j]= true;
prev=i;
i=j;
}returnret;
}public List GetConcave_Edge(doubleradius)
{
List ret = new List();
List[] adjs = GetInRNeighbourList(2 *radius);
ret.Add(points[0]);int i = 0, j = -1, prev = -1;while (true)
{
j=GetNextPoint_EdgePivoting(prev, i, adjs[i], radius);if (j == -1)break;//Point2d p = BallConcave.GetCircleCenter(points[i], points[j], radius);
ret.Add(points[j]);
flags[j]= true;
prev=i;
i=j;
}returnret;
}private bool CheckValid(List[] adjs)
{for (int i = 0; i < adjs.Length; i++)
{if (adjs[i].Count < 2)
{return false;
}
}return true;
}public boolCompareAngel(Point2d a, Point2d b, Point2d m_origin, Point2d m_dreference)
{
Point2d da= new Point2d(a.X - m_origin.X, a.Y -m_origin.Y);
Point2d db= new Point2d(b.X - m_origin.X, b.Y -m_origin.Y);double detb =GetCross(m_dreference, db);//nothing is less than zero degrees
if (detb == 0 && db.X * m_dreference.X + db.Y * m_dreference.Y >= 0) return false;double deta =GetCross(m_dreference, da);//zero degrees is less than anything else
if (deta == 0 && da.X * m_dreference.X + da.Y * m_dreference.Y >= 0) return true;if (deta * detb >= 0) {//both on same side of reference, compare to each other
return GetCross(da, db) > 0;
}//vectors "less than" zero degrees are actually large, near 2 pi
return deta > 0;
}public int GetNextPoint_EdgePivoting(int prev,int current,List list,doubleradius)
{if (list.Count == 2 && prev != -1)
{return list[0] + list[1] -prev;
}
Point2d dp;if (prev == -1)
dp= new Point2d(1, 0);elsedp=Point2d.Minus(points[prev], points[current]);int min = -1;for (int j = 0; j < list.Count; j++)
{if (!flags[list[j]])
{if (min == -1)
{
min=list[j];
}else{
Point2d t=points[list[j]];if (CompareAngel(points[min], t, points[current], dp) && GetDistance(t,points[current])
{
min=list[j];
}
}
}
}//main.ShowMessage("seek P" + points[min].Index);
returnmin;
}public int GetNextPoint_BallPivoting(int prev, int current, List list, doubleradius)
{
SortAdjListByAngel(list, prev, current);for (int j = 0; j < list.Count; j++)
{if(flags[list[j]])continue;int adjIndex =list[j];
Point2d xianp=points[adjIndex];
Point2d rightCirleCenter=GetCircleCenter(points[current],xianp,radius);if (!HasPointsInCircle(list,rightCirleCenter,radius,adjIndex))
{
main.DrawCircleWithXian(rightCirleCenter, points[current], points[adjIndex], radius);returnlist[j];
}
}return -1;
}private void SortAdjListByAngel(List list, int prev, intcurrent)
{
Point2d origin=points[current];
Point2d df;if (prev != -1)
df= new Point2d(points[prev].X - origin.X, points[prev].Y -origin.Y);elsedf= new Point2d(1, 0);int temp = 0;for (int i = list.Count; i > 0; i--)
{for (int j = 0; j < i - 1; j++)
{if (CompareAngel(points[list[j]], points[list[j + 1]], origin, df))
{
temp=list[j];
list[j]= list[j + 1];
list[j+ 1] =temp;
}
}
}
}private bool HasPointsInCircle(List adjPoints,Point2d center,double radius,intadjIndex)
{for (int k = 0; k < adjPoints.Count; k++)
{if (adjPoints[k] !=adjIndex)
{int index2 =adjPoints[k];if(IsInCircle(points[index2], center, radius))return true;
}
}return false;
}public static Point2d GetCircleCenter(Point2d a, Point2d b, doubler)
{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 Point2d(-1, -1);
}double sqrt = Math.Sqrt(r * r / (dx * dx + dy * dy) - 0.25);return new Point2d(cx - dy * sqrt, cy + dx *sqrt);
}public static bool IsInCircle(Point2d p,Point2d center,doubler)
{double dis2 = (p.X - center.X) * (p.X - center.X) + (p.Y - center.Y) * (p.Y -center.Y);return dis2 < r *r;
}public List[] GetInRNeighbourList(doubleradius)
{
List[] adjs = new List[points.Count];for (int i = 0; i < points.Count; i++)
{
adjs[i]= new List();
}for (int i = 0; i < points.Count; i++)
{for (int j = 0; j < points.Count; j++)
{if (i < j && distanceMap[i, j]
{
adjs[i].Add(j);
adjs[j].Add(i);
}
}
}returnadjs;
}private List GetSortedNeighbours(intindex)
{
List infos = new List(points.Count);for (int i = 0; i < points.Count; i++)
{
infos.Add(newPoint2dInfo(points[i], i, distanceMap[index, i]));
}
infos.Sort();
List adj = new List();for (int i = 1; i < infos.Count; i++)
{
adj.Add(infos[i].Index);
}returnadj;
}public static doubleGetDistance(Point2d p1, Point2d p2)
{return Math.Sqrt((p1.X - p2.X) * (p1.X - p2.X) + (p1.Y - p2.Y) * (p1.Y -p2.Y));
}public static doubleGetCross(Point2d a, Point2d b)
{return a.X * b.Y - a.Y *b.X;
}
}
求解这个问题理论上还有很多思路,这些思路日后有时间还会继续在这里更新。其实笔者曾经将这个问题联系到旅行商问题和哈密顿回路。因为在有的应用场景,当人们给出一系列点,这些点来表征物体轮廓时,可能会要求所有给出的点都必须在生成的多边形上。假如有这要求的话alpha shape之类的方案就不成了,那或许就是一个旅行商问题,不过貌似旅行商问题不能保证不自交的路径,这方面笔者还没有试过,将来有时间也许会尝试一下。
其实讨论这个问题的一大动机还是为了进一步探讨这一问题扩展到三维的情况,业界在这一方面也有不少成熟的研究,例如三维滚球法,PowerCrust,Cocone,假如能知道点法向,还能使用重建更加光滑局面的隐式曲面的方法如MLS,RBF,POISSON等,对三维的方法的探讨将来还会继续更新。
二维点集求外轮廓Java_从二维点集重建平面形状-浅议凹包算法相关推荐
- 轮廓提取之滚球法、凹包算法
提示:本文章主要介绍轮廓提取算法------滚球法的运用,可用于提取点云等二维.三维数据的轮廓,以及代码框架设计. 文章目录 前言 一.基于二维点数据的模型轮廓提取算法 1.1 适用场景 1.2 问题 ...
- matlab 凹包,二维点集凹包算法介绍
最近遇到一个求二维点集凹包的问题,凹包的叫法不知道是否准确,问题可以描述为:(原文下载在文章末尾) 在二维平面上有一系列的点,求能包围所有点集的二维多边形.(好像搜"离散点边界"或 ...
- 采用CAD刀具插件(飞狼刀模工具箱)的求外轮廓线功能快速提取包装盒刀模图
1.打开需要绘制刀模图的包装盒盒型图,或通过工具箱里的盒型库功能事先绘制盒型. 2.通过工具箱的求外轮廓线快速求出外形,操作步骤为:点击求外轮廓线图标->选择盒型图->按右键. 3.通过板 ...
- Opencv的使用小教程3——利用轮廓检测实现二维码定位
Opencv的使用小教程3--利用轮廓检测实现二维码定位 二维码具有什么特征 实现效果 识别二维码的流程 1.预处理图像 2.寻找轮廓 3.通过寻找到的轮廓确定"回"的位置 4.创 ...
- OpenCV-透视变换及对二维点求透视变换之后的坐标
OpenCV-透视变换及对二维点求透视变换之后的坐标 微信公众号:幼儿园的学霸 个人的学习笔记,关于OpenCV,关于机器学习, -.问题或建议,请公众号留言; 目录 前言 基本原理 OpenCV透视 ...
- 《C#零基础入门之百识百例》(五十二)封装介绍 -- 二维多项式求值
C#零基础入门 面向对象 -- 封装介绍 -- 二维多项式求值 前言 一,封装概念 二,封装属性 三,实例练习 -- 二维多项式求值 3.1 题目描述 3.2 问题分析 3.3 参考代码 前言 本文属 ...
- c语言二维数组对角线输出字符,用C编程一个4*4的矩阵,用二维数组,求对角线元素的和。,c语言如何使二维数组 输出为矩阵的形式...
导航:网站首页 > 用C编程一个4*4的矩阵,用二维数组,求对角线元素的和.,c语言如何使二维数组 输出为矩阵的形式 用C编程一个4*4的矩阵,用二维数组,求对角线元素的和.,c语言如何使二维数 ...
- 二维搜索求函数极值的python程序
二维搜索,也就是优化两个未知变量,通常将二维搜索转化为一维搜索进行求解. 例题 二维搜索转一维--黄金分割 求函数二维函数A在某区间的最大值 A=( 6 - 2×l + l×cos(θ) ) × l× ...
- C语言求二维数组平均数,一道JavaScript的二维数组求平均数的题
JavaScript中只支持一维数组,但是可以在数组中嵌套数组来创建二维以至于多维的数组.今天下午在看书时候,发现一道感觉比较有意思的题,就是js中如何求二维数组的列之和和行之和,现在就给大家分享下, ...
最新文章
- 北京交通大学计算机科学与技术研究生导师,熊轲_北京交通大学研究生导师信息...
- 派对屋3000效果器怎样调试_怎样测试自己的信用卡是否进了小黑屋
- 删除排序链表中的重复元素Python解法
- 安卓逆向_15( 三 ) --- Android NDK 开发【 jni 静态注册、JNI_OnLoad 动态注册】
- 发生身份验证错误_Python MySQL8.0 错误
- 机器学习-吴恩达-笔记-13-大规模机器学习
- hdu 3926 hands in hands
- 快手新作品监控提醒助手,可无水印下载视频和封面
- Cloudera Manager安装之利用parcels方式安装单节点集群(包含最新稳定版本或指定版本的安装)(添加服务)(CentOS6.5)(四)...
- python 实现实时语音对讲
- Android Platform 3.0 SDK和Eclipse ADT安装记录三
- 类型多样的游戏特效网页特效素材,速来收藏
- 数学建模-对策论模型
- 1. Java基础语法
- MFC_改变窗口大小,随着鼠标的拖动改变
- 一文看懂“知识蒸馏”技术
- 怎么把CAD表格转换成Excel表格?CAD表格导出教程
- python双下划线怎么打出来_Python中单、双下划线的区别总结
- 英文字母和汉字的字节长度问题
- 计算机专硕双导师,说说你不知道的双导师制