凸包算法 Matlab实现

内容参考 https://blog.csdn.net/viafcccy/article/details/87483567

n_points = 25; % total number of points
points_source = generateRandomPoints(n_points, 500, 0.12);
% n_points = length(x_coor);
% points_source = [y_coor,x_coor];
% test
% a=[1,2,3,4,5,6;2,3,4,5,6,8]';
% n_points = length(a(:,1));
% points_source = [a(:,1),a(:,2)];
points = points_source;
hull = [];
inner = [];
undefined = (1 : n_points)';
% As the origin point of the convex hull we take point with the lowest Y value.
% Origin point is a priori lies on the border of the convex hull.
[~, lowest_point_id] = min(points(:, 2));

% Calculate angles between X axis and each segment, connecting the origin
% with all of the points.
u = points(lowest_point_id,1:2);
angles = zeros(n_points,1);
for i = 1 : n_points
    v = points(i, 1 : 2);    
    w = v - u;
    angles(i) = atan2d(w(2), w(1));
end
% Sort points in according to it's angles with the X axis.
undefined = [undefined,angles,points(:,1)];
undefined = sortrows(undefined, 3);
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 points 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_points = points(test_ids,:);
        
        if isCounterclockwise(test_points(1,:),test_points(2,:),test_points(3,:))
            break
        end
        inner = cat(1,inner,undefined(1));
        undefined(1)=[];
        undefined = cat(1, hull(end), undefined);
        hull(end,:) = [];
    end
    hull = cat(1,hull,undefined(1));
    undefined(1)=[];
end
hull = cat(1,hull,undefined(1));
hfig = figure(1);
plot(points(:,1),points(:,2),'.k');
hold on;
plot(points(hull,1),points(hull,2),'.-r');
hold on;
sysHull = convhull(points(:,1),points(:,2));
plot(points(sysHull,1),points(sysHull,2),'.-b');
function result = isCounterclockwise(a, b, c)
    result = ((b(1) - a(1)) * (c(2) - a(2)) - (b(2) - a(2)) * (c(1) - a(1)) >= 0);
end