基于matlab的凸包(Convex Hull)算法理解与测试

基于matlab的凸包Convex Hull算法理解与测试

0 引言
1 Graham扫描算法原理
2 Graham扫描算法实现
3 Graham扫描算法关键函数
4 结语

0 引言

💻💻AI一下💻💻

凸包算法是计算给定点集的最小凸包的一种算法。凸包是包含给定点集中所有点的最小凸多边形。根据具体的实现和算法思想，凸包算法可以分为以下几类：

(1) Jarvis算法（也称为包裹算法）：该算法通过找到初始点，然后沿着点集中的边界逐步扩展凸包，直到回到初始点为止。这种算法的时间复杂度为O(nh)，其中n是点集的大小，h是凸包的边数。

(2) Graham扫描算法：该算法首先选取一个最低点作为起始点，然后根据与起始点的极角进行排序。然后依次将每个点加入凸包中，如果加入后出现了凹角，则将凸包中的一些点删除，直到不再出现凹角。

(3) QuickHull算法：该算法通过递归地将点集分成较小的两部分，并构造凸包。具体过程是：选择距离最远的两点作为凸包的两个顶点，然后将点集划分为两个子集，分别在凸包的左侧和右侧。然后对两个子集分别进行递归，直到遍历了所有的点。该算法的时间复杂度取决于点集的分布，平均情况下为O(nlogn)。

(4) Chan算法（也称为增量算法）：该算法结合了Jarvis算法和Graham扫描算法的优点。首先将点集分成多个子集，在每个子集上运用Graham扫描算法构造凸包。然后通过Jarvis算法将所有的子凸包合并成一个凸包。该算法的时间复杂度为O(nlogh)，其中h是凸包的边数。

以上是一些常见的凸包算法分类和介绍，根据具体的应用场景和性能要求，可以选择合适的算法来计算凸包。

本篇结合网上搜集到的资料和程序，介绍下Graham扫描法的原理和实现过程。

1 Graham扫描算法原理

上面讲到了Graham扫描算法的思想是以一点为基准，最好是特殊点（如左下角点、中心点），然后根据与起始点的极角进行排序，排序之后每3个一组遍历所有点集，然后判断3个点的位置关系是否满足凸包原则，把满足的点留下，不满足的点从点集中删除，每删除一个点，都要从新遍历所有点集，直至推出整个遍历过程，那留下的点所构成的图形即为包围所有点集的凸多边形。下面用一个简单案例，详细描述下上述过程：

(1) 为了便于描述，用下面5个点进行原理描述。首先，确定基准点，本过程选取O为基准；

(2) 计算其它点与 $O$ 点的极角 $(α)$ ，并按大小进行排序，下图的排序顺序为 $[D E A BC]$ 的顺序，图中 $XO D$ 即为 $D$ 点对应的极角；

(3) 遍历排序数组，每次拿3个出来，我们以 $A BC$ 为例，以 $A$ 为交点，判断向量 $\vec{AB}$ 和 $\vec{AC}$ 的位置关系，由图可以看出，向量 $\vec{AC}$ 可由向量 $\vec{AB}$ 逆时针旋转得到，或向量 $\vec{AB}$ 在向量 $\vec{AC}$ 右侧，而B点明显是凸点，所以依据该准则能够找到所有满足条件的点;

(4) 再以 $C D E$ 为例，以 $C$ 点为角点，向量 $\vec{CD}$ 到向量 $\vec{CE}$ 需要顺时针旋转，所以 $D$ 点就不是凸点，应该从点集中删除;

(5) 用向量叉乘描述向量 $\vec{AB}$ 和向量 $\vec{AC}$ 的位置关系。具体以下图为例，假设向量 $\vec{OY}=[0,1,0]$ ，向量 $\vec{OA} = [1,1,0]$ ，向量 $\vec{OB} = [-1,1,0]$ ，计算 $\vec{OY} \times \vec{OA = [0,0 ,1]}$ ，计算 $\vec{OY} \times \vec{OB = [0,0 ,-1]}$ ，由叉乘结果的正负就可以判断位置关系，负值为逆时针，正值为顺时针。对应到(3)和(4)就是留下向量叉积小于0的点，删除叉积结果大于0的点，上面5个点中 $D$ 点是要删除的点;

2 Graham扫描算法实现

💦💦💦💦💦

clc;clear
n = 50;
point = rand(n,2);figure
scatter(point(:,1),point(:,2));
konvexHullPoints = GrahamScanAlgorithm(point);
%% Visualise konvex Hull
drawHull(konvexHullPoints(:,1),konvexHullPoints(:,2));

50个随机点构成的凸多边形结果展示。
在这里插入图片描述

3 Graham扫描算法关键函数

💦💦💦💦💦

function [KonvexHullPoints] = GrahamScanAlgorithm(points)
%GRAHAMASCANALGORITHM Summary of this function goes here
xMax = max(points(:,1));
xMin = min(points(:,1));
yMax = max(points(:,2));
yMin = min(points(:,2));Z.x = xMin + (xMax-xMin)/2;
Z.y = yMin + (yMax-yMin)/2;
% scatter(Z.x,Z.y,'*'); %turn on for visualisation of center point%% calculate all arcs between the randpoints and center Z
degHold = []; %holds all the degrees 
for(a=1:1:length(points))degHold(end+1) = calcArc(Z,points(a,:));
end
degHold = transpose(degHold);%% Sort Points 
pointIndice = transpose(1:1:length(points));
degHold = table(degHold, pointIndice,points(:,1),points(:,2));
degHoldSort = sortrows(degHold,'degHold','ascend');
degHoldSort.Properties.VariableNames{3} = 'x';
degHoldSort.Properties.VariableNames{4} = 'y';
KonvexHull = degHoldSort;
%% Calculate Convex Hull
noMoreRemove = false; %condition for staying in while-loop
k = 1;% Increment
cas = 0; % case while(noMoreRemove == false ) %stay in while until no points need to be removed anymoreif(k<=height(KonvexHull)-2)%standard case far away from the vector endP1 = [KonvexHull{k,3},KonvexHull{k,4}];P2 = [KonvexHull{k+1,3},KonvexHull{k+1,4}];P3 = [KonvexHull{k+2,3},KonvexHull{k+2,4}];cas = 0;elseif(k==height(KonvexHull)-1)%k is the penultimate element must now be k+1 = last vector element and k+2 == 1st vector elementP1 = [KonvexHull{k,3},KonvexHull{k,4}];P2 = [KonvexHull{k+1,3},KonvexHull{k+1,4}];P3 = [KonvexHull{1,3},KonvexHull{1,4}];cas = 0;elseif(k==height(KonvexHull))%k is the last element must now be k+1 == 1st vector element and k+2 == 2nd vector element P1 = [KonvexHull{k,3},KonvexHull{k,4}];P2 = [KonvexHull{1,3},KonvexHull{1,4}];P3 = [KonvexHull{2,3},KonvexHull{2,4}];cas = 1;endpolyTest = polyCheck(P1,P2,P3);if(polyTest <= 0 && cas == 0)% concave must be removeKonvexHull(k+1,:) = []; %removek = 0;elseif(polyTest < 0 && cas == 1)KonvexHull(1,:) = []; %removek = 0;elseif(k==height(KonvexHull)&&polyTest>0)     noMoreRemove = true;break;endk = k+1;
endKonvexHullPoints = [KonvexHull{:,3},KonvexHull{:,4}];
end

function [ArcDeg] = calcArc(Z,point)
%CALCARC Summary of this function goes here
toDeg = 180/pi;
dY = point(2) - Z.y;
dX = point(1) - Z.x;
ArcRad  = atan2(dY,dX); if(ArcRad < 0)  ArcDeg = 360 - ((abs(ArcRad)/pi)*180);
elseArcDeg = ArcRad * toDeg;
endend

function  drawHull(HullX,HullY)
%DRAWHULL Summary of this function goes here
%    Plot the wrapping polygon
%% Andreas Bernatzky 19.08.2019
hold(gca,'on');plot(HullX,HullY,'r','LineStyle','--','LineWidth',2);
lastPartX = [HullX(end);HullX(1)];
lastPartY = [HullY(end);HullY(1)];
plot(lastPartX,lastPartY,'r','LineStyle','--','LineWidth',2);
sz = 50;
scatter(HullX,HullY,sz,'*');
%title('Graham-Scan Algorithm for calculating a convex hullpolygon around a pointcloud');
%set(gca,'FontSize',16,'FontWeight','bold');
end