function [meassurf meassurfN0]=croftonDW(X,r,Ntheta,Ny,N0)

L=max(pdist(X));
d=size(X,2);
n=size(X,1);
Betad=gamma(d/2)/(sqrt(pi)*gamma((d+1)/2));

Theta=normrnd(0,1,Ntheta,d);
norTheta=sqrt(sum(Theta.^2,2));
Theta=Theta./(norTheta*ones(1,d));
Theta(:,1)=abs(Theta(:,1));
angle=[0:2*pi/100:2*pi];    


meassurf=0;
meassurfN0=0;
for thet=1:Ntheta
    thet;
    %simulation of Y in Theta perp
    Y=unifrnd(-L,L,Ny,d-1);
    Mat=[Theta(thet,:)' unifrnd(0,1,d,d-1)];
    [Q R]=qr(Mat);
    YY=(Q(:,2:d)*Y')';
    YY*Theta(thet,:)';
    for iy=1:Ny
        %list DWr
        ListPossInter=[];
        Listell=[];
        for i=1:n
            uiy=X(i,:)-YY(iy,:);
            di=sqrt(norm(uiy)^2-(uiy*Theta(thet,:)')^2);
            if di<=r;
                ell1=uiy*Theta(thet,:)'-sqrt(r^2-di^2);
                ell2=uiy*Theta(thet,:)'+sqrt(r^2-di^2);
            
                z1=YY(iy,:)+ell1*Theta(thet,:);
                z2=YY(iy,:)+ell2*Theta(thet,:);
                
                ZZ=[z1;z2];
                LL=[ell1;ell2];
                
                Idx=knnsearch(X,[z1;z2]);
                keep=zeros(2,1);
                for j=1:2
                    if Idx(j)==i
                        keep(j)=1;
                    end
                end
                ii=find(keep==1);
                if isempty(ii)==0
                    ZZ=ZZ(ii,:);
                    LL=LL(ii,:);
                    ListPossInter=[ListPossInter;ZZ];
                    Listell=[Listell;LL];
                end
            end
        end
        [Listellr Ielr]=sort(Listell);
        Extremr=ListPossInter(Ielr,:);
        
        %list DW4r
        ListPossInter=[];
        Listell=[];
        for i=1:n
            uiy=X(i,:)-YY(iy,:);
            di=sqrt(norm(uiy)^2-(uiy*Theta(thet,:)')^2);
            if di<=4*r;
                ell1=uiy*Theta(thet,:)'-sqrt(16*r^2-di^2);
                ell2=uiy*Theta(thet,:)'+sqrt(16*r^2-di^2);
            
                z1=YY(iy,:)+ell1*Theta(thet,:);
                z2=YY(iy,:)+ell2*Theta(thet,:);
                
                ZZ=[z1;z2];
                LL=[ell1;ell2];
                
                Idx=knnsearch(X,[z1;z2]);
                keep=zeros(2,1);
                for j=1:2
                    if Idx(j)==i
                        keep(j)=1;
                    end
                end
                ii=find(keep==1);
                if isempty(ii)==0
                    ZZ=ZZ(ii,:);
                    LL=LL(ii,:);
                    ListPossInter=[ListPossInter;ZZ];
                    Listell=[Listell;LL];
                end
            end
        end        
        Listell4r=sort(Listell);
        
        NCr=size(Listellr,1);
        KCr=NCr/2;
        CCr=zeros(KCr,2);
        cc=0;
        for i=1:KCr
            for j=1:2
                cc=cc+1;
                CCr(i,j)=Listellr(cc);
            end
        end
        
        
        NC4r=size(Listell4r,1);
        KC4r=NC4r/2;
        CC4r=zeros(KC4r,2);
        cc=0;
        for i=1:KC4r
            for j=1:2
                cc=cc+1;
                CC4r(i,j)=Listell4r(cc);
            end
        end        
               
        numcc4r=ones(KCr,1);
        for i=1:KCr
            for j=1:KC4r
                if and(CCr(i,1)>=CC4r(j,1),CCr(i,2)<=CC4r(j,2))
                    numcc4r(i)=j;
                end
            end
        end
        UU=unique(numcc4r);
        
          
        meassurf=meassurf+2*size(unique(numcc4r),1);
        meassurfN0=meassurfN0+2*min(size(unique(numcc4r),1),N0);
        
    end
end

            
meassurf=(meassurf*(2*L)^(d-1))/(Betad*Ny*Ntheta);    
meassurfN0=(meassurfN0*(2*L)^(d-1))/(Betad*Ny*Ntheta);        
    
end