function H=PHull(X)
%Convex hull algorithm for points in 3D (the columns of X).
%The convex hull is returned as the two field structure H 
%
%   H(k).node  is the column number in X of a convex hull node.
%   H(k).nbr   contains column numbers in X of some of the convex hull 
%              nodes that share an edge with H(k).node.  The remaining 
%              edges connected to H(k).node are given in the nbr field for
%              other elements of H.  Overall there is no redundancy.
%
%This code may crash given theoretically unlikely collections of points 
%such as: 
%
% 1. more than three points that are exactly co-planar.
% 2. more than one point with the minimum x1 coordinate.
%
%These cases occur frequently in practice (e.g., with CMM data to be 
%analyzed for flatness).  The code can still be used if the data is
%perturbed at the level of numerical noise (e.g., 1e-12).
%
%Patch Kessler, August 14, 2008.

n=size(X,2);
free=true(1,n); 
%               Find the first point:
[minx,I1]=min(X(1,:)); 
I1=I1(1);    
H(1).node=I1; 
H(2).nbr=I1; 
free(I1)=false;
%               Find the second point:
a=X-X(:,I1)*ones(1,n);
beta=a(1,:);                   
alpha=sqrt(sum(a([2 3],:).^2));
q1=atan2(beta,alpha);
q2=q1(free);
[minang,K]=min(q2); 
K=find(cumsum(free)==K(1)); 
I2=K(1);
H(2).node=I2; 
H(3).nbr=I2; 
free(I2)=false;
%               Find the third point:
e1=X(:,I2)-X(:,I1);                 
e1=e1/norm(e1);
e2=[1 0 0]'-e1(1)*e1; 
e2=e2/norm(e2);
b=a-e1*(e1'*a);
beta=e2'*b;
alpha=sqrt(sum((b-e2*beta).^2));
q1=atan2(beta,alpha);
q2=q1(free);
[minang,K]=min(q2); 
K=find(cumsum(free)==K(1)); 
I3=K(1);
H(1).nbr=I3;
H(3).node=I3; 
free(I3)=false;
%               Construct the initial ring:
R{1}=[I1 I2 I3;I3 I1 I2];           
if (X(:,I3)-X(:,I1))'*cross(e1,e2)>0
    R{1}=[I1 I3 I2;I2 I1 I3];
end

while 1
   [R,H,free]=growring(R,H,free,X,n);
   if isempty(R) break; end
end   

%---------------------------------------------
function [R,H,free]=growring(R,H,free,X,n)
r=R{end};
free2=free; free2(r(1,3:end))=true; free2(r(2,1))=false;
u=X(:,r(1,2))-X(:,r(1,1));  eu=u/norm(u);
v=X(:,r(2,1))-X(:,r(1,1));  ev=v-eu*(eu'*v);  ev=-ev/norm(ev);
ew=cross(ev,eu);
alpha=ev'*(X-X(:,r(1,1))*ones(1,n));
beta=ew'*(X-X(:,r(1,1))*ones(1,n));
free2(beta<=0)=false;
q1=atan2(beta,alpha); 
q2=q1(free2);
[minang,K]=min(q2); K=find(cumsum(free2)==K(1)); K=K(1);

if ~ismember(K,r(1,:))
    %case 1: min point is not a part of the current ring. 
    r=[[r(1,1);r(1,2)] [K;r(1,1)] r(:,2:end)];
    R{end}=r; 
    free(K)=false;
    H(length(H)+1).node=K;
    H(end).nbr=[r(1,1) r(1,3)];
elseif size(r,2)==3
    %case 2: current (triangular) ring can be closed.
    R=R(1:end-1);
elseif K==r(1,3)
    %case 3: min point is in ring, at node 3.
    R{end}=[[r(1,1);r(1,2)] r(:,3:end)]; 
    t=find([H.node]==K);
    H(t).nbr=[H(t).nbr r(1,1)];
elseif K==r(1,end)
    %case 4: min point is in ring, at end node.
    R{end}=[[r(1,end);r(1,1)] r(:,2:end-1)]; 
    t=find([H.node]==K);
    H(t).nbr=[H(t).nbr r(1,2)];
else
    %case 5: min point is in ring, between node 3 and end node.
    s=find(K==r(1,:));
    R{end}=[[r(1,1);r(2,1)] r(:,s:end)];
    R{end+1}=[r(:,2:s-1) [r(1,s);r(1,1)]];
    t=find([H.node]==K);
    H(t).nbr=[H(t).nbr r(1,1) r(1,2)];
end