function [tf,DE,DF]=flatness(X)
%compute flatness of the columns of X.

%Compute the convex hull.  
H=PHull(X+1e-10*(rand(size(X))-0.5));   %Perturb input to numerical precision.

%Edges
EN=numel([H.nbr]);
E=zeros(2,EN); 
index=1;
for k=1:length(H)
    n=length(H(k).nbr);
    E(1,index:index+n-1)=H(k).node;
    E(2,index:index+n-1)=H(k).nbr;
    index=index+n;
end
DE=inf;
optE=zeros(2);  %columns give indices of closest edges
for k1=1:EN
    for k2=k1+1:EN
        if E(1,k1)==E(1,k2) | E(1,k1)==E(2,k2) | E(2,k1)==E(1,k2) | E(2,k1)==E(2,k2) 
            continue
        end
        x1=X(:,E(1,k1));
        e1=X(:,E(2,k1))-x1;  e1=e1/norm(e1); 
        x2=X(:,E(1,k2));
        e2=X(:,E(2,k2))-x2;  e2=e2/norm(e2);
        e3=cross(e1,e2);
        if norm(e3)<1e-10, disp('warning!! e3 is too small.'), continue, end
        e3=e3/norm(e3);  
        s1=sign(e3'*(X-x1*ones(1,size(X,2)))); s1(E(:,k1))=0; 
        s2=sign(e3'*(X-x2*ones(1,size(X,2)))); s2(E(:,k2))=0;
        if any(s1.*s2>0), continue, end
        d=abs(e3'*(x1-x2));
        if d<DE, DE=d; optE=E(:,[k1 k2]); end
    end
end

%Faces
DF=inf;
optF=zeros(1,4);    %first three indecies identify triangle, last one identifies point.
QF=[];   %find all triangular rings of edges (faces and more).
for k=1:length(E)
    T=[E(1,E(2,:)==E(1,k)) E(2,E(1,:)==E(1,k))]; T=T(T~=E(2,k));
    B=[E(1,E(2,:)==E(2,k)) E(2,E(1,:)==E(2,k))]; B=B(B~=E(1,k));
    I=intersect(T,B);
    QF=[QF;[ones(length(I),1)*E(:,k)' I']];
end
for k=1:length(QF)
    a=X(:,QF(k,1))-X(:,QF(k,2));
    b=X(:,QF(k,3))-X(:,QF(k,2));
    e=cross(a,b);
    e=e/norm(e);
    d=e'*(X-X(:,QF(k,2))*ones(1,length(X)));
    if any(d>0)&any(d<0), continue, end
    [d,Id]=max(abs(d));
    if d<DF, DF=d; optF=[QF(k,:) Id]; end
end

tf=min(DE,DF);