計算幾何-凸包算法 Python實現與Matlab動畫演示

凸包算法是計算幾何中的最經典問題之一了。給定一個點集，計算其凸包。凸包是什么就不羅嗦了

本文給出了《計算幾何——算法與應用》中一書所列凸包算法的Python實現和Matlab實現，並給出了一個Matlab動畫演示程序。

啊，實現誰都會實現啦╮(╯▽╰)╭，但是演示就不一定那么好做了。

算法CONVEXHULL(P) 
輸入：平面點集P 
輸出：由CH(P)的所有頂點沿順時針方向組成的一個列表
1.   根據x-坐標，對所有點進行排序，得到序列p1, …, pn
2.   在Lupper中加入p1和p2（p1在前）
3. for(i←3 ton) 
4.   do 在Lupper中加入pi
5.   while(Lupper中至少還有三個點，而且最末尾的三個點所構成的不是一個右拐) 
6.   do 將倒數第二個頂點從Lupper中刪去
7.   在Llower 中加入pn和pn-1（pn在前）
8. for(i←n-2 downto1) 
9.   do 在Llower 中加入pi
10.    while(Llower 中至少還有三個點，而且最末尾的三個點所構成的不是一個右拐) 
11.    do 將倒數第二個頂點從Llower 中刪去
12.  將第一個和最后一個點從Llower 中刪去
（以免在上凸包與下凸包聯接之后，出現重復頂點）
13.  將Llower 聯接到Lupper后面（將由此得到的列表記為L）
14.  return(L) 

看看，不是很多的樣子是吧。
這里面需要說明的地方只有一點，那就是方向的判定問題。

設有三個點P，Q，R，現需求R位於向量PQ的左側還是右側(或者R在PQ上)。

計算PR與PQ的叉積，也就是外積。

如果將向量PR以P為旋轉中心旋轉只向量PQ的方向走過一個正角(逆時針)，意味着R在PQ的右側，此時外積為正。

另外需要注意的是，如果在凸包上有三點共線的情況，在本例中三點是均位於凸包邊界點集中的。如果想避免這一點，可以通過微量抖動數據的方式解決。

廢話不多說，Python實現如下：

#!/usr/bin/env python
# coding: gbk
 
########################################################################
#Author: Feng Ruohang
#Create: 2014/10/02 13:39
#Digest: Computate the convex hull of a given point list
########################################################################
 
 
direction = lambda m: (m[2][0] - m[0][0]) * (m[1][1] - m[0][1]) - (m[1][0] - m[0][0]) * (m[2][1] - m[0][1])
'''
    A Quick Side_check version Using Lambda expression
    Input:  Given a list of three point : m should like [(p_x,p_y), (q_x,q_y), (r_x,r_y)]
    Output: Return a Number to indicate whether r on the right side of vector(PQ).
    Positive means r is on the right side of vector(PQ).
    This is negative of cross product of PQ and PR: Defined by:(Qx-Px)(Ry-Py)-(Rx-Px)(Qy-Py)
    Which 'negative' indicate PR is clockwise to PQ, equivalent to R is on the right side of PQ
'''
 
 
def convex_hull(point_list):
    '''
    Input:  Given a point List: A List of Truple (x,y)
    Output: Return a point list: A List of Truple (x,y) which is CONVEX HULL of input
    For the sake of effeciency, There is no error check mechanism here. Please catch outside
    '''
    n = len(point_list)  #Total Length
    point_list.sort()
 
    #Valid Check:
    if n < 3:
        return point_list
 
    #Building Upper Hull: Initialized with first two point
    upper_hull = point_list[0:1]
    for i in range(2, n):
        upper_hull.append(point_list[i])
        while len(upper_hull) >= 3 and not direction(upper_hull[-3:]):
            del upper_hull[-2]
 
    #Building Lower Hull: Initialized with last two point
    lower_hull = [point_list[-1], point_list[-2]]
    for i in range(n - 3, -1, -1):  #From the i-3th to the first point
        lower_hull.append(point_list[i])
        while len(lower_hull) >= 3 and not direction(lower_hull[-3:]):
            del lower_hull[-2]
    upper_hull.extend(lower_hull[1:-1])
    return upper_hull
 
 
#========Unit Test:
if __name__ == '__main__':
    test_data = [(i, i ** 2) for i in range(1, 100)]
    result = convex_hull(test_data)
    print result

2015年1月23日

使用很簡單，看DocString就行。

下面順便給出了Matlab 的實現，以及可視化的算法演示：

效果就是個小動畫，像這樣吧。

Matlab的凸包算法有三個文件：

side_check2：檢查三個點構成的彎折的方向

Convex_hull: 凸包算法Matlab實現

Convex_hull_demo：凸包算法的演示。

拷在一個目錄里

運行convex_hull_demo( randn(200,2)*100); 就可以看到可視化演示了

這個是輔助函數

%filename: side_check2.m
%Input:     Matrix of three point: (2x3 or 3x2)
%           P(p_x,p_y),Q(q_x,q_y),R(r_x,r_y)
%Output:    如果P Q R三點構成一個右拐，返回True
%           右拐意味着點R在PQ向量的右側.此時

function result = side_check2(D)
if all(size(D) ~= [3,2])
    if all(size(D)==[2,3])
        D = D';
    else
        error('error dimension')
    end
end
result = (det([[1;1;1], D]) < 0 );

這個是純算法實現。

%filename:     convex_hull.m
%CONVEX_HULL
%INPUT:     Point Set:(n x 2)
%OUPUT:     HULL Point List: (x x 2)
function L=t(P)
[num,dimension] = size(P);
if dimension ~= 2
    error('dimension error')
end

P = sortrows(P,[1,2]);
%if there is only one or two point remain,return it
if num < 3
     L = P;
     return 
end


%STEP ONE: Upper Hull:
L_upper = P([1,2],:); %Take first two points
for i = 3:num
    L_upper = [L_upper;P(i,:)]; %add the point into list 
    while size(L_upper,1) >= 3
        l_size = size(L_upper,1);
        if det([ones(3,1),L_upper(l_size-2:l_size,:)])= 3
        l_size = size(L_lower,1);
       if det([ones(3,1),L_lower(l_size-2:l_size,:)])

這個是演示：

%CONVEX_HULL
%INPUT:     Point Set:(n x 2)
%OUPUT:     HULL Point List: (x x 2)
%Samples:   convex_hull_demo( randn(200,2)*100)

function L=convex_hull_demo(P)

%Test Data
%data_size = data_size
%P = randi([-50,50],[data_size,2]);
[num,dimension] = size(P);
if dimension ~= 2
    error('dimension error')
end


P = sortrows(P,[1,2]);

%====Visual Lization
board_left = min(P(:,1));
board_right = max(P(:,1));
board_bottom = min(P(:,2));
board_up = max(P(:,2));
x_padding = (board_right- board_left)*0.1;
y_padding = (board_up- board_bottom)*0.1;
plot_range= [board_left - x_padding,board_right + x_padding,board_bottom-y_padding,board_up+y_padding];

clf;
scatter(P(:,1),P(:,2),'b.');
axis(plot_range);
hold on
%====VisualLization

%if there is only one or two point remain,return it
if num < 3
     L = P;
end

%STEP ONE: Upper Hull:
L_upper = P([1,2],:); %Take first two points
hull_handle = plot(L_upper(:,1),L_upper(:,2),'ob-');
for i = 3:num
    L_upper = [L_upper;P(i,:)]; %add the point into list
    
    while size(L_upper,1) >= 3
        l_size = size(L_upper,1);
        if side_check2(L_upper(l_size-2:l_size,:)) %Check if it is valid
            break;  %Quit if Valid
        else
            L_upper(l_size-1,:) = []; %Remove the inner point and continue if not
        end
        set(hull_handle,'XData',L_upper(:,1),'YData',L_upper(:,2));drawnow;
        
    end
    set(hull_handle,'XData',L_upper(:,1),'YData',L_upper(:,2));drawnow;
end
 
 
 %Visualization
 plot(L_upper(:,1),L_upper(:,2),'bo-');
 %Visualization

  
%STEP Two: Build the lower hull
L_lower = [P([num,num-1],:)]; % Add P(n) and P(n-1)
set(hull_handle,'XData',L_lower(:,1),'YData',L_lower(:,2));drawnow;


for i = num-2:-1:1
    L_lower = [L_lower;P(i,:)];
    while size(L_lower,1) >= 3
        l_size = size(L_lower,1);
       if side_check2(L_lower(l_size-2:l_size,:)) %Check if it is valid
            break;  %Quit if Valid
        else
            L_lower(l_size-1,:) = []; %Remove the inner point and continue if not
       end    
       set(hull_handle,'XData',L_lower(:,1),'YData',L_lower(:,2));drawnow;
    end
    set(hull_handle,'XData',L_lower(:,1),'YData',L_lower(:,2));drawnow;
end

L_lower([1,size(L_lower,1)],:) = [];
if isempty(L_lower)
    L = L_upper;
else
    L = [L_upper;L_lower(2:size(L_lower,1)-1,:)];
end
hold off;
return