计算几何-凸包算法 Python实现与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) 

 1 #!/usr/bin/env python
 2 # coding: gbk
 3  
 4 ########################################################################
 5 #Author: Feng Ruohang
 6 #Create: 2014/10/02 13:39
 7 #Digest: Computate the convex hull of a given point list
 8 ########################################################################
 9  
10  
11 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])
12 ‘‘‘
13     A Quick Side_check version Using Lambda expression
14     Input:  Given a list of three point : m should like [(p_x,p_y), (q_x,q_y), (r_x,r_y)]
15     Output: Return a Number to indicate whether r on the right side of vector(PQ).
16     Positive means r is on the right side of vector(PQ).
17     This is negative of cross product of PQ and PR: Defined by:(Qx-Px)(Ry-Py)-(Rx-Px)(Qy-Py)
18     Which ‘negative‘ indicate PR is clockwise to PQ, equivalent to R is on the right side of PQ
19 ‘‘‘
20  
21  
22 def convex_hull(point_list):
23     ‘‘‘
24     Input:  Given a point List: A List of Truple (x,y)
25     Output: Return a point list: A List of Truple (x,y) which is CONVEX HULL of input
26     For the sake of effeciency, There is no error check mechanism here. Please catch outside
27     ‘‘‘
28     n = len(point_list)  #Total Length
29     point_list.sort()
30  
31     #Valid Check:
32     if n < 3:
33         return point_list
34  
35     #Building Upper Hull: Initialized with first two point
36     upper_hull = point_list[0:1]
37     for i in range(2, n):
38         upper_hull.append(point_list[i])
39         while len(upper_hull) >= 3 and not direction(upper_hull[-3:]):
40             del upper_hull[-2]
41  
42     #Building Lower Hull: Initialized with last two point
43     lower_hull = [point_list[-1], point_list[-2]]
44     for i in range(n - 3, -1, -1):  #From the i-3th to the first point
45         lower_hull.append(point_list[i])
46         while len(lower_hull) >= 3 and not direction(lower_hull[-3:]):
47             del lower_hull[-2]
48     upper_hull.extend(lower_hull[1:-1])
49     return upper_hull
50  
51  
52 #========Unit Test:
53 if __name__ == ‘__main__‘:
54     test_data = [(i, i ** 2) for i in range(1, 100)]
55     result = convex_hull(test_data)
56     print result
57 
58 2015年1月23日

%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,:)])