%% Parameters and data source

pnseq='C:\Users\gracec\Documents\PICS DATA\YR151026\Sequences\';
pnout='C:\Users\gracec\Documents\PICS DATA\YR151026\auto_threshold_processed';
if ~exist(pnout,'dir')
    mkdir (pnout)
end
%files={'20151026132438.seq'
% '20151026132407.seq'/
% '20151026132330.seq'};

files=dir(fullfile(pnseq,'*.seq'));  %List of sequence files

%Parameters fro local thresdholding and edge detection
GSlim=[0,100]; %grayscale display contrast [0,255]
%local intensity thresholding settings
p.opts_LT.sense=0; %sensitivity of ADAPTTHRESH [0-1]
p.opts_LT.stat='gaussian'; %gaussian,mean,median
%edge detection settings
p.opts_Edge.method='laplacian'; %'laplacian','log','prewitt','sobel'
p.opts_Edge.alpha=0.2;  %Laplacian shape factor [0-1]
p.opts_Edge.logN=[5,5]; %Neighborhood size for Gaussian part of LoG method
p.opts_Edge.sigma=0.5;  %std [>0] LoG method
p.opts_Edge.dir='vert'; %vert or horiz for Prewitt or Sobel methods
%focus thresholds
p.focus.L90=0.02; %90th percentile Laplacian threshold%
p.focus.PoP=10; %percent of perimeter threshold

%sharp.Radius=5; % Default= 1.The size of the region around the edge pixels that is affected by sharpening. A large value sharpens wider regions and a small value sharpens narrower regions around edges.
p.sharp.Radius=2;
p.sharp.Amount=1.28841; %[0 2],Default=0.8; Strength of sharpening effect
%p.sharp.Amount=2;
p.sharp.Thold=0.05; %[0 1], Default=0; Threshold:Minimum contrast requred for a pixel to be considered an edge pixel. 

%Load image infromation
load(fullfile(pnseq,'post.mat')); %get registration info

for k=16:length(files)
%for k=3:5
%for k=4
    tic
    display(['Starting File ', num2str(k), ' out of ',num2str(length(files))])
    [hdr]=load_SEQ(fullfile(pnseq,files(k).name));
    files(k).name
    for J=1:hdr.AllocatedFrames %frame index
        [hdr,I]=load_SEQ(fullfile(pnseq,files(k).name),J);
        
%%Automatic global thresholding for estimate of largest particle the Otsu threshold level
        gt=graythresh(I);
        I1=im_morph(I,gt); %binary image
        [R1,L1]=im_props(I1,post.A,J);
        tmp=cat(1,R1.xSphDia);
        Kmx=find(tmp==max(tmp),1,'first');
        Nmx=R1(Kmx).BoundingBox(3:4);
        Max_d_otsu=R1(Kmx).xSphDia; %esd in mm

%%Sharpen image edges of particles chosen from local threshold
        Isharp=imsharpen(I,'Amount',p.sharp.Amount,'Radius',p.sharp.Radius,'Threshold',p.sharp.Thold);

    %%Local intensity thresholding: This section determines candidate particles for future tracking neighborhood must be odd-integer sized
        Nsize=ceil(1.1*Nmx);
        Nsize=Nsize+~mod(Nsize,2);
 
        if Nsize(1,1) >1000
            Nsize=[9 9];
        end

        %Nsize=[13 13];


    T=adaptthresh(I,p.opts_LT.sense,...
        'ForegroundPolarity','bright',...
        'Statistic',p.opts_LT.stat,...
        'NeighborhoodSize',Nsize);

    Tsharp=adaptthresh(Isharp,p.opts_LT.sense,...
        'ForegroundPolarity','bright',...
        'Statistic',p.opts_LT.stat,...
        'NeighborhoodSize',Nsize);

    I2=imbinarize(I,T);
    I2sharp=imbinarize(Isharp,Tsharp);

    I2=imfill(I2,'holes');
    I2sharp=imfill(I2sharp,'holes');

    [R2,L2]=im_props(I2,post.A,J);
    [R2sharp,L2sharp]=im_props(I2sharp,post.A,J);

    %%Edge detection

   %%form the spatial filter structure
    switch lower(p.opts_Edge.method)
        case 'laplacian'
            h=fspecial('laplacian',p.opts_Edge.alpha);
        case 'log'
            h=fspecial('log',p.opts_Edge.logN,p.opts_Edge.sigma);
        case {'prewitt','sobel'}
            h=fspecial(p.opts_Edge.method);
     switch lower(p.opts_Edge.dir)
        case {'vert','vertical'}
            h=h';
     end
    end
    
%    %filter the image
    I3=imfilter(im2double(I),h,'symmetric');

 %  % Metrics for determining accepted particles
    med=nan(size(R2)); %prealloc
    PixelNum=nan(size(R2)); %prealloc
    Perimeter=nan(size(R2)); %prealloc
    data_I3=cell(size(R2)); %prealloc
    L90=nan(size(R2)); %prealloc
    PoP=nan(size(R2)); %prealloc

for kk=1:length(R2)
    med(kk)=double(median(I(R2(kk).PixelIdxList)));
    PixelNum(kk)=length(R2(kk).PixelIdxList);
    Perimeter(kk)=(R2(kk).Perimeter);
    data_I3{kk}=I3(R2(kk).PixelIdxList(:));
    if PixelNum(kk)>=9
        [cdf,edges]=histcounts(data_I3{kk},'Normalization','cdf');
        L90(kk)=findcumval(cdf,edges(2:end),0.90);
        PoP(kk)=100*sum(data_I3{kk}>p.focus.L90)/Perimeter(kk);
    end
end

  p.L90{J}=cat(2,L90);
  p.PoP{J}=cat(2,PoP);

 %   define accepted partilces 
    in=L90>=p.focus.L90 & PoP>=p.focus.PoP;
    p.in{J}=cat(2,in);
    %%in2=L90>=focus.L90 & PoP>=3 & PoP<focus.PoP; %Extend .L90
    in3=PixelNum<9; %small particles
    in4=PixelNum>9; %
    p.Lg(J)=length(R2(in4));
    p.total(J)=length(R2);
    p.count(J)=length(R2(in));
    p.xSphDia{J}=cat(2,R2.xSphDia);
    p.xSphDiaA{J}=cat(2,R2(in).xSphDia);

    p.PixelList{J,:}=cat(1,R2.PixelList);
    p.PixelListA{J,:}=cat(1,R2(in).PixelList);

     p.R2_accepted{J,:}=R2(in);
     %p.R2{J,:}=R2;

%%Metrics for determining accepted particles (Sharpend Applied)

    med_sharp=nan(size(R2sharp)); %prealloc
    PixelNum_sharp=nan(size(R2sharp)); %prealloc
    Perimeter_sharp=nan(size(R2sharp)); %prealloc
    data_I3_sharp=cell(size(R2sharp)); %prealloc
    L90_sharp=nan(size(R2sharp)); %prealloc
    PoP_sharp=nan(size(R2sharp)); %prealloc


    for kk=1:length(R2sharp)
        med_sharp(kk)=double(median(I(R2sharp(kk).PixelIdxList)));
        PixelNum_sharp(kk)=length(R2sharp(kk).PixelIdxList);
        Perimeter_sharp(kk)=(R2sharp(kk).Perimeter);
        data_I3_sharp{kk}=I3(R2sharp(kk).PixelIdxList(:));
        if PixelNum_sharp(kk)>=9
            [cdf,edges]=histcounts(data_I3_sharp{kk},'Normalization','cdf');
            L90_sharp(kk)=findcumval(cdf,edges(2:end),0.90);
            PoP_sharp(kk)=100*sum(data_I3_sharp{kk}>p.focus.L90)/Perimeter_sharp(kk);
        end
    end

    
  p.L90_sharp{J}=cat(2,L90_sharp);
  p.PoP_sharp{J}=cat(2,PoP_sharp);

%define accepted partilces (compare if I looked at sharp that passed the credentials). 
    in=L90_sharp>=p.focus.L90 & PoP_sharp>=p.focus.PoP;
    p.in_sharp{J}=cat(2,in);
    %%in2=L90_sharp>=focus.L90 & PoP_sharp>=3 & PoP_sharp<focus.PoP; %Extend .L90
    in4=PixelNum_sharp>9; %
    in3=PixelNum_sharp<9; %small particles
    p.Lg_sharp(J)=length(R2sharp(in4));
    p.total_sharp(J)=length(R2sharp);
    p.count_sharp(J)=length(R2sharp(in));

    p.xSphDia_sharp{J}=cat(2,R2sharp.xSphDia);
    p.xSphDia_sharpA{J}=cat(2,R2sharp(in).xSphDia);

    p.PixelList_sharp{J,:}= cat(1,R2sharp.PixelList);
    p.PixelList_sharpA{J,:}= cat(1,R2sharp(in).PixelList);

    p.R2sharp_accepted{J,:}=R2sharp(in);
    p.R2sharp{J,:}=R2sharp;
    
   % [alpha,beta,~,alpha_sharp,beta_sharp]=get_alpha_beta_sobel(I,R2,R2sharp);

    % p.alpha{J}=cat(2,alpha);
     %p.beta{J}=cat(2,beta);

     %p.alpha_sharp{J}=cat(2,alpha_sharp);
     %p.beta_sharp{J}=cat(2,beta_sharp);

    
    end
    
    p.esd=cat(2,p.xSphDia{:});%equivelent spherical diameter (mm)
    p.esd_sharp=cat(2,p.xSphDia_sharp{:});%equivelent spherical diameter (mm)

    p.esdA=cat(2,p.xSphDiaA{:});%equivelent spherical diameter (mm)
    p.esd_sharpA=cat(2,p.xSphDia_sharpA{:});%equivelent spherical diameter (mm)

    
   p.fnseq=files(k).name;
   save(fullfile(pnout,['p_',files(k).name(1:end-4)]),'-struct','p'); %same for output file
toc 
end 