function [trit Rt pvalMS]=alphashapemaxspace2(X,kdens);

N=size(X,1);
d=size(X,2);

tri=delaunayn(X);
[R,C]=rayonscir(tri,X);
[Rt,It]=sort(R);
trit=tri(It,:);
Nt=size(trit,1);

Vol=zeros(Nt,1);
Vol(1)=volumesimplexe2(X(trit(1,:),:));

for i=2:Nt
    Vol(i)=Vol(i-1)+volumesimplexe2(X(trit(i,:),:));
end

[Idens, Rdens]=knnsearch(X,X,'K',kdens+1);
Ddens=Rdens(:,kdens).^(-1/d);

[XX Ind Pds]=unifsuprnd3(trit,X,N*100);

dd1=Ddens(trit(Ind,:));
DXX=zeros(N*100,1);
for i=1:N*100
    DXX(i)=dd1(i,:)*Pds(i,:)';
end

[ii ra]=knnsearch(X,XX);
va=ra.^d;
MaxSpa=zeros(Nt,1);
pvalMS=2*ones(Nt,1);
alpha=((sqrt(pi)*gamma((d+2)/2)/gamma((d+1)/2))^(d-1))/(gamma(d+1));
omega=pi^(d/2)/gamma((d+2)/2);
for i=1:100*N
    j=Ind(i);
    MaxSpa(i)=0.5*(kdens/N)*(ra(i)/Rdens(ii(i),kdens+1))^d;
    U=N*MaxSpa(i)-log(N)-(d-1)*log(log(N))-log(alpha);
    pv=min(1,1-exp(-exp(-U)));
    %pause
    if pv<pvalMS(j)
        pvalMS(j)=pv;
    end
end

for i=1:Nt
    pvalMS(i)=min(pvalMS(1:i));
end
ii=find(pvalMS>=0.005);
i0=max(ii);
% figure
% subplot(1,2,1);
% plot(Rt,pvalMS,'.')
% hold on 
% plot([Rt(i0) Rt(i0)],[0 1],'-r')
% Rt(i0)
% %axis([0 0.015 0 1])
% axis([0 1 0 1])
% subplot(1,2,2)
% if d==2
% trisurf(trit(1:i0,:),X(:,1),X(:,2),zeros(N,1))
% end
% if d==3
% tetramesh(trit(1:i0,:),X)
% end
