平面點集的凸包計算

平面點集的凸包可理解為包含所有點的最小凸多邊形（點可以在多邊形邊上或在其內）。這里給出一種求解方法。

一、基本思路

先找所有點中 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

四、運行示例