平面点集的凸包计算

平面点集的凸包可理解为包含所有点的最小凸多边形（点可以在多边形边上或在其内）。这里给出一种求解方法。

一、基本思路

先找所有点中 y 坐标最大最小的点P_max、P_min，所找点必定是凸包上的点；

找距离直线P_maxP_min两侧最远的点P₁,P₀，构成初始三角形, ；

再对每个三角形新生成的边（、和、）继续找与改变对应顶点（）不在同一侧的最远点。

二、算法流程

1 找所有点中 y 坐标最大和最小的点

     1.1 若找到的点少于两个,return,输出（无凸包结构）

     1.2 若y坐标最大最小点各只有一个记为Pmax,Pmin,找直线PmaxPmin两侧最远的点P₁,P₀,将构成的三角形, 放入堆栈TriStack

     1.3 若找到的点大于两个,把这些点能组成的三角形放入堆栈TriStack

2 若TriStack不为空

     2.1 三角形出栈，找三角形前两个顶点的对边与该点异侧的最远点

     2.2 若点存在,边与点组成三角形放入TriStack

     2.3 若点不存在,该边存入Boundary,返回2

3 返回 Boundary

 

三、实现代码

该算法由matlab实现：

clc;
clear;
N = 74;
DataPoints = [(1:N)', rand(N, 2).*100];
plot(DataPoints(:, 2), DataPoints(:, 3), '.');
% grid on;
X = DataPoints(:,2);
Y = DataPoints(:,3);

%% 
% 根据 y 方向最值点确立初始三角形
% 找 y 最大和最小值的点 --------
[Ymax, Ymax_i] = max(Y);
[Ymin, Ymin_i] = min(Y);
PymaxPymin = [X(Ymin_i) - X(Ymax_i), Y(Ymin_i) - Y(Ymax_i)];
PtNum = N;
% 找距离直线 PymaxPymin 两侧最远点 --------
PtNum_p = 1; PtNum_n = 1;
Dis_p = 0; Dis_n = 0;
for j = 1 : PtNum
    PjPymax = [X(Ymax_i) - X(j), Y(Ymax_i) - Y(j)];
    PjPymin = [X(Ymin_i) - X(j), Y(Ymin_i) - Y(j)];
    Tri_erea = det([PjPymax; PjPymin]);
    if Tri_erea < 0
        if Dis_n < abs(Tri_erea)
            Dis_n = abs(Tri_erea);
            PtNum_n = j;
        end
    else
        if Dis_p < Tri_erea
            Dis_p = Tri_erea;
            PtNum_p = j;
        end
    end
end

% 计算凸包边界 ----------
TriStack = [];
TriStack(1, :) = [Ymax_i, Ymin_i, PtNum_p];
TriStack(2, :) = [Ymax_i, Ymin_i, PtNum_n];
Boundary = FindBoundary(TriStack, DataPoints);
hold on;
for i = 1 : size(Boundary, 1)
    plot([X(Boundary(i,1)), X(Boundary(i,2))], ...
    [Y(Boundary(i,1)), Y(Boundary(i,2))], '-');
end

function TriPtNum_new = TriFindPt(TriPtNum, DataPoints)
% 扩展三角形的两边
% TriPtNum [i,j,k]是三角形 PiPjP 的顶点序号, 过顶点 P 的两边要扩展
% DataPoints 是所有点的坐标
% TriPtNum_new [i,j]是两条扩展边最外侧顶点序号, 若不存在取0

X = DataPoints(:,2);
Y = DataPoints(:,3);
% 点 Pi, Pj的对边, 扩展 -------------------
PjP = [X(TriPtNum(3)) - X(TriPtNum(2)), Y(TriPtNum(3)) - Y(TriPtNum(2))];
PiP = [X(TriPtNum(3)) - X(TriPtNum(1)), Y(TriPtNum(3)) - Y(TriPtNum(1))];
PiPj = [X(TriPtNum(2)) - X(TriPtNum(1)), Y(TriPtNum(2)) - Y(TriPtNum(1))];
% PiPj x PiP, PjPi x PjP, 向量叉积
PiPjxPiP = det([PiPj; PiP]);
PjPixPjP = det([-PiPj; PjP]);

% 找点 Pi, Pj 的对边距离最远点 ------------------ 
PtNum = length(X);
PtNum_pi = 0; PtNum_pj = 0;
Dis_pi = 0; Dis_pj = 0;
for k = 1 : PtNum
    % 点 Pi 的对边找最远点
    PkPj = [X(TriPtNum(2)) - X(k), Y(TriPtNum(2)) - Y(k)];
    PkP = [X(TriPtNum(3)) - X(k), Y(TriPtNum(3)) - Y(k)];
    PkPjxPkP = det([PkPj; PkP]);
    % 点 Pj 的对边找最远点
    PkPi = [X(TriPtNum(1)) - X(k), Y(TriPtNum(1)) - Y(k)];
    PkPixPkP = det([PkPi; PkP]);
    
    if(PiPjxPiP * PkPjxPkP < 0)
        Tri_ereai = abs(PkPjxPkP);
        if(Dis_pi < Tri_ereai)
            Dis_pi = Tri_ereai;
            PtNum_pi = k;
        end
    end
    
    if(PjPixPjP * PkPixPkP < 0)
        Tri_ereaj = abs(PkPixPkP);
        if(Dis_pj < Tri_ereaj)
            Dis_pj = Tri_ereaj;
            PtNum_pj = k;
        end
    end
    TriPtNum_new = [PtNum_pi, PtNum_pj];
end

function Boundary = FindBoundary(TriStack, DataPoints)
% 从初始三角形开始查找边界

Boundary = [];
while size(TriStack, 1)
    TriPtNum_new = TriFindPt(TriStack(end, :), DataPoints);
    TriPt = TriStack(end, :);
    TriStack(end, :) = [];
    if TriPtNum_new(1)
        TriStack(end+1, :) = [TriPt(:, 2:3), TriPtNum_new(1)];
    else
        Boundary(end+1, :) = TriPt(:, 2:3);
    end
    if TriPtNum_new(2)
        TriStack(end+1, :) = [TriPt(1, 1), TriPt(1, 3), TriPtNum_new(2)];
    else
        Boundary(end+1, :) = [TriPt(1, 1), TriPt(1, 3)];
    end
end

四、运行示例