凸包算法理解——基于MATLAB代码

个人有个习惯：学习一个算法，要有代码能单步逐行理解。

正好这两天研究了一下凸包算法。留个记录文档吧。

代码如下，有详细备注。

% 凸包算法理解  
clear all;close all;clc;

%% 数据源
data_x=[-2.96,-3.26,-3.26,-2.96,-3.23,-2.93,-2.98,-3.28,-3.83,-4.13,-4.13,-3.83,-3.44,-3.74,-3.74,-3.44,-3.93,-4.23,-4.23,-3.93,-3.96,-4.26];
data_y=[29.21,29.21,28.91,28.91,28.65,28.65,29.46,29.46,28.86,28.86,28.56,28.56,29.42,29.42,29.12,29.12,29.6,29.6,29.3,29.3,29.85,29.85];

data=[data_x;data_y]';
hfig=plot_number(data_x,data_y);  %..绘制数据并标记ID
data_number = length(data_x);
hull = [];  %..记录凸包的ID
inner = [];  %..记录内点的ID
undefined = (1 : data_number)';

[~, lowest_y_id] = min(data_y);  %..查找y最小值点，将该点作为起始点。

% Calculate angles between X axis and each segment, connecting the origin
% with all of the data.
u = data(lowest_y_id,1:2);  %..起始点的数据。
angles = zeros(data_number,1);
for i = 1 : data_number
    v = data(i, 1 : 2);    
    w = v - u;
    angles(i) = atan2d(w(2), w(1));%..所有点到起始点的夹角。
end
% Sort data in according to it's angles with the X axis.
undefined = [undefined,angles,data(:,1)];
undefined = sortrows(undefined, 3); %..先针对x值的大小进行排序
undefined = sortrows(undefined, 2); %..再根据角度的大小进行排序
undefined = undefined(:, 1);  %..排序之后的序号

% First point will be on the convex hull.
hull = undefined(1);  %..将待处理数据集中的第一个点（就是上面的起始点）作为凸包点。
undefined = [undefined; undefined(1)];
undefined(1) = [];%..将起始点放在待处理点集的最后。

% Loop through all data and remove those of them which are not on the border
% of the convex hull.
while length(undefined)>=2
    while true
        test_ids = [hull(end);undefined(1:2)]; %..凸包点集中最后一个点，与待检测点集中前两个点组成test点集。
        test_data = data(test_ids,:); %..待检测点集的数据。
        
        if Direction_determine(test_data(1,:),test_data(2,:),test_data(3,:))  %..根据叉积判断当前3个点集的方式是顺时针还是逆时针。
            break
        end
        inner = cat(1,inner,undefined(1)); %..将待处理点集中第一个元素归到内点点集中。记录内点的ID
        undefined(1)=[]; %..从待处理点集中删除已经归到内点点集的点。
        undefined = cat(1, hull(end), undefined); %..将凸包点集最后一个点重新拉回到待检测点集，并置于首位。
        hull(end,:) = [];%..从凸包点集中删除已经归到待处理点集中的点。
    end
    %..3个点处于逆时针的位置，则可以将待处理点集中第一个元素暂时归到凸包点集中。
    hull = cat(1,hull,undefined(1));  %..归入凸包点集
    undefined(1)=[]; %从待处理点集中删除已经归到凸包点集的点。
end
hull = cat(1,hull,undefined(1)); %..将待处理点集中最后一个点归到凸包点集中。

%..绘制凸包
plot(data(hull,1),data(hull,2),'.-r');
hold on;

function result = Direction_determine(a, b, c)   
% 3个点组成两个向量。然后求两个向量的叉积，并根据叉积的方向确定两个向量的关系。
% 根据右手定则，逆时针向上为正，顺时针向下为负。（需要理解叉集的基本概念）
%
    x1=b(1) - a(1); %..第一个向量是（x1,y1）
    y1=b(2) - a(2);
    x2=c(1) - a(1); %..第二个向量是（x2,y2）
    y2=c(2) - a(2);
    result = ((x1 * y2) - (x2 * y1) >= 0);
end

function [hfig] = plot_number(X,Y)
    %plot_number 以绘制X,Y二维散点图并标出序号
    hfig=figure(100);
    scatter(X,Y,'r.')%绘制散点图
    hold on
    for i=1:max(size(X))
        c = num2str(i);%数字转字符
        text(X(i),Y(i),c);%在图上显示文字
    end
    hold on;
end