clear all;
clc;
close all;
disp('Analysis of STED data.m...');

%% ------Variables to Edit--- %%
savefolderName='Analysis_New4';
threshArea=75; %kernel size for adaptive threshold
threshOffset=-0.1; %more negative = higher threshold
scale=0.0196;%in microns

indexOnLines=0; %Do you want index numbers on the 'line' images: 1 = Yes, 0=No

%% -----Main Programme----- %%
d = dir(pwd);
isub = [d(:).isdir]; %# returns logical vector
nameFolds = {d(isub).name}';
nameFolds(ismember(nameFolds,{'.','..'})) = [];

parentfolder=pwd;

for mainloop=1:length(nameFolds)
output=[];
distanceMaster=[];
areaCh0=[];
area2Ch1=[];

clearvars -except outpit distanceMaster indexOnLines nameFolds parentfolder savefolderName threshArea threshOffset scale mainloop areaCh0 area2Ch1

folder=strcat(pwd,'\',nameFolds(mainloop),'\');
    
cd(char(folder));
ch0 = dir('*C=0.tif');
ch1=dir('*C=1.tif');
for file = 1:length(ch0)
     close all

ch1Name = ch1(file).name;
info = imfinfo(ch1Name);
numberOfImagesch1 = length(info);
[pathstr,namCh1,ext] = fileparts(ch1Name);
name=namCh1(1:end-6);
tempIm=imread(ch1Name);
[sizeX,sizeY]=size(tempIm);
imageStack = zeros(sizeX,sizeY,numberOfImagesch1); 

for k = 1:numberOfImagesch1
    currentImage = imread(ch1Name, k, 'Info', info);
    imageStack(:,:,k) = currentImage;
end 
mipch1Original = max(imageStack, [], 3);
mipch1Original=mat2gray(mipch1Original);

mipch1 = max(imageStack, [], 3);
mipch1=mat2gray(mipch1);
mipch1=imtranslate(mipch1,[0, -1]);%CORRECTION COLOUR
%figure(),imshow(mipch1);
mipch1=wiener2(mipch1,[5,5]);%Wiener filter size
f1=mipch1;
f1=double(f1);
[M,N]=size(f1);

F1=fft2(f1)/sqrt(M*N);
cF1=fftshift(F1); % centering in the origin
cP1=abs(cF1).^2;

T1=0.1; %Fourier threshold
c1=F1.*(abs(F1)>=T1);
fM1=real(ifft2(c1))*sqrt(M*N);
mipch1=fM1;
%figure();imshow(mipch1);


ch0Name = ch0(file).name;
info2 = imfinfo(ch0Name);
numberOfImagesch0 = length(info2);
[pathstr,namech0,ext] = fileparts(ch0Name);
tempIm=imread(ch0Name);
[sizeX2,sizeY2]=size(tempIm);
imageStack2 = zeros(sizeX2,sizeY2,numberOfImagesch0);

for k = 1:numberOfImagesch0
    currentImage2 = imread(ch0Name, k, 'Info', info2);
    imageStack2(:,:,k) = currentImage2;
end 

mipch0Original = max(imageStack2, [], 3);
mipch0Original=mat2gray(mipch0Original);

mipch0 = max(imageStack2, [], 3);
mipch0=mat2gray(mipch0);
%figure(),imshow(mipch0);

mipch0=wiener2(mipch0,[5,5]);
f=mipch0;
f=double(f);
[M,N]=size(f);
F=fft2(f)/sqrt(M*N);
cF=fftshift(F);
cP=abs(cF).^2;

T=0.1;
c=F.*(abs(F)>=T);
fM=real(ifft2(c))*sqrt(M*N);
mipch0=fM;
%figure();imshow(mipch0);


rotation=0;
fractionAssoc=[];
xshiftMaster=[];
yshiftMaster=[];

for masterCount=1:1
    
    if masterCount>1
        xshift=round(rand*512);
        yshift=round(rand*512);
        
        xshiftMaster=[xshiftMaster;xshift];
        yshiftMaster=[yshiftMaster;yshift];
        
    mipch1 = circshift(mipch1, [xshift, yshift]);

    else mipch1=mipch1;
    end 
  


    seFat=strel('disk',1);
    seFat2=strel('disk',1);
    
    mask=adaptivethreshold(mipch1,threshArea,threshOffset,0);%was -0.15
    mask=bwareaopen(mask,20);    
    bw=mask;
    %figure();imshow(bw);
 
    %Ch1 WATERSHED SEPARATION
    D = -bwdist(~bw);
    Ld = watershed(D);
    bw2 = bw;
    bw2(Ld == 0) = 0;
    extendMin = imextendedmin(D,0.35);
    D2 = imimposemin(D,extendMin);
    Ld2 = watershed(D2);
    bw3 = bw;
    bw3(Ld2 == 0) = 0;
    bw3=bwareaopen(bw3,20);
    %figure();imshow(bw3);
   
    mean=[];
    IntInfobw3=regionprops(bw3,mipch1,'MeanIntensity','MaxIntensity','MinIntensity','Centroid','PixelList','Area');
    if isempty(IntInfobw3)
        break
    end
    centroidsch1 = cat(1, IntInfobw3.Centroid);


    %ch0 SEGMENTATION
    
    mask2=adaptivethreshold(mipch0,threshArea,threshOffset,0);
    
    mask2=bwareaopen(mask2,20);%was50

    E = -bwdist(~mask2);
    ELd = watershed(E);
    E2 = mask2;
    E2(ELd == 0) = 0;
    extendMin2 = imextendedmin(E,0.35);
    E3 = imimposemin(E,extendMin2);
    ELd2 = watershed(E3);
    E4 = mask2;
    E4(ELd2 == 0) = 0;
    E4=bwareaopen(E4,20);
   
    IntInfoE4=regionprops(E4,mipch0,'MeanIntensity','MaxIntensity','MinIntensity','Centroid','PixelList','Area');
    if isempty(IntInfoE4)
        break
    end
    centroidsch0 = cat(1, IntInfoE4.Centroid);
        
    %[indexClose,distance]=knnsearch(centroidsch1,centroidsch0,'k',1,'distance','euclidean');
    [indexClose,distance]=knnsearch(centroidsch1,centroidsch0,'k',20,'distance','euclidean');
    indexClose2=indexClose;
    indexClose2(distance>25.51)=0;
    if all(indexClose2==0)
        break
    end
    indexClose3=table(indexClose2);
    distance2=distance;
    distance2(distance>25.51)=0; %25.51 = 500um search radius
    distance3=table(distance2);
    %figure();imshow(imfuse(mipch0Original,mipch1Original));
%     [r,c,d] = size(mipch0Original);              %# Get the image size
%     set(gca,'Units','normalized','Position',[0 0 1 1]);  %# Modify axes size
%     set(gcf,'Units','pixels','Position',[0 0 c r]);  %# Modify figure size
    
     % set(gcf,'PaperPosition',[0,0,40,40]);  

%     rez=1200; %resolution (dpi) of final graphic
%     f=gcf; %f is the handle of the figure you want to export
%     figpos=getpixelposition(f);
%     resolution=get(0,'ScreenPixelsPerInch');
%     set(f,'paperunits','inches','papersize',figpos(3:4)/resolution,'paperposition',[0 0 figpos(3:4)/resolution]);
   
    blank1=zeros(sizeX,sizeY);
    blank2=blank1;
    blank3=blank1;
    fused=imfuse(mipch0Original,mipch1Original);
    %RGB=repmat(mipch0Original,[1,1,3]);
    linesMas=[];
    count=1;
    for i=1:length(IntInfoE4)
        for j=1:length(indexClose2(2,:))
            
            if (indexClose2(i,j)~=0)
            data1=IntInfoE4(i).Centroid;data2=IntInfobw3(indexClose2(i,j)).Centroid;
            %plot(imgca,data1(:,1),data1(:,2),'r*');plot(imgca,data2(:,1),data2(:,2),'b*');
            LineDataX=[data1(1);data2(1)];LineDataY=[data1(2);data2(2)];
            lines=[LineDataX(1),LineDataY(1),LineDataX(2),LineDataY(2)];
            linesMas=[linesMas;lines];
            markerInserter = vision.MarkerInserter('Shape','Plus','BorderColor','white');
            blank1=step(markerInserter, blank1, [lines(1),lines(2)]);

            markerInserter2=vision.MarkerInserter('Shape','Star','BorderColor','white');
            blank3=step(markerInserter2,blank3,[lines(3),lines(4)]);

            %disnumber=count;
            disnumber=i;
            disnumber_string = sprintf('%d',disnumber);
            centroid=IntInfoE4(i).Centroid;
            if indexOnLines==1
            H=vision.TextInserter(disnumber_string);
            H.Color=[1,1,1];
            H.FontSize=12;
            H.Location=[(centroid(1)+2),centroid(2)+2];
            blank2=step(H,blank2);
            end
            
            centroidnumCh0=i;
            centroidnumCh0_string=sprintf('%d',centroidnumCh0);
            ch0Index=vision.TextInserter(centroidnumCh0_string);
            ch0Index.Color=[1,1,1];
            ch0Index.FontSize=12;
            ch0Index.Location=[(centroid(1)+1),centroid(2)+1];
            blank1=step(ch0Index,blank1);

           %centroidnumCh1=indexClose2(i);
            centroidnumCh1=indexClose2(i,j);
            centroidnumCh1_string=sprintf('%d',centroidnumCh1);
            %centroid2=IntInfobw3(indexClose2(i)).Centroid;
            centroid2=IntInfobw3(indexClose2(i,j)).Centroid;
            ch1Index=vision.TextInserter(centroidnumCh1_string);
            ch1Index.Color=[1,1,1];
            ch1Index.FontSize=12;
            ch1Index.Location=[(centroid2(1)+1),centroid2(2)+1];
            blank3=step(ch1Index,blank3);

            count=count+1;
            end
        end
        
    end
 
        indexing=[];
        distance2=distance(indexClose2~=0)*scale;
        nameArray={name};
        nameArray2=repmat(nameArray,length(distance2),1);
        for k=1:length(IntInfoE4)
            if(indexClose2(k)~=0)
                indexing=[indexing;k];
            end
        end
        indexing=[1:length(distance2)];
        indexing=indexing.';
        
        disout=[nameArray2,num2cell(indexing),num2cell(distance2)];
        distanceMaster=[distanceMaster;disout];

        area=[IntInfobw3.Area];
        area=area.';
        area=area.*scale.*scale;
        nameArrayCh0=repmat(nameArray,length(IntInfobw3),1);
        indexingCh0=[1:length(IntInfobw3)];
        indexingCh0=indexingCh0.';
        areaout=[nameArrayCh0,num2cell(indexingCh0),num2cell(area)];
%         
        area2=[IntInfoE4.Area];
        area2=area2.';
        area2=area2.*scale.*scale;
        nameArrayCh1=repmat(nameArray,length(IntInfoE4),1);
        indexingCh1=[1:length(IntInfoE4)];
        indexingCh1=indexingCh1.';
        area2out=[nameArrayCh1,num2cell(indexingCh1),num2cell(area2)];
        
        areaCh0=[areaCh0;areaout];
        area2Ch1=[area2Ch1;area2out];
    
        matrix=[]
        matrix=zeros(length(area2)+1,length(area)+1);
        for index1=1:length(area2);
            for index2=1:length(area);
                matrix(index1+1,1)=index1;
                matrix(1,index2+1)=index2;
            end
        end
        
        for k=1:length(indexClose);
            for h=1:length(indexClose(1,:))
               
                if indexClose2(k,h)~=0;
                    matrix(k+1,indexClose2(k,h)+1)=distance(k,h);
                end
            end
        end
        
        
        area2Mat=[0,area2.'];
        area2Mat=area2Mat.';
        matrix2=[matrix,area2Mat];
        areaMat=[0,area.',0];
        matrix3=[matrix2;areaMat];
    
        matrix4=table(matrix3);
        
        white = uint8([255 255 255]);
        shapeInserter = vision.ShapeInserter('Shape','Lines','BorderColor','White');%,'CustomBorderColor',white);
        blank2=step(shapeInserter, blank2, linesMas);
        RGB=cat(3,mipch1Original,mipch0Original,blank2);
        mkdir(char(folder),char(savefolderName));
        savedir=strcat(char(folder),char(savefolderName),'\');
        imwrite(RGB,strcat(savedir,name,'_originalMerge.tif'))
        RGBfiltered=cat(3,mipch1,mipch0,blank2);
        imwrite(RGBfiltered,strcat(savedir,name,'_filteredMerge.tif'));
        imwrite(mipch0Original,strcat(savedir,name,'_layered.tif'));   
        imwrite(mipch1Original,strcat(savedir,name,'_layered.tif'),'WriteMode','append');
        imwrite(blank2,strcat(savedir,name,'_layered.tif'),'WriteMode','append');
        imwrite(blank1,strcat(savedir,name,'_layered.tif'),'WriteMode','append');
        imwrite(blank3,strcat(savedir,name,'_layered.tif'),'WriteMode','append');
        

        
        
%CENTROID AND LINE IMAGES        
%         figure();imshow(imfuse(mipch0Original,mipch1Original),[])
%         hold(imgca,'on')
%         plot(imgca,centroidsch0(:,1), centroidsch0(:,2), 'r*')
%         plot(imgca,centroidsch1(:,1), centroidsch1(:,2), 'b+');
%         hold(imgca,'off')
%         
%         figure();imshow(imfuse(mipch0Original,mipch1Original));
%         hold(imgca,'on')
%         for i=1:length(IntInfoE4)
%             if indexClose2(i)~=0
%             data1=IntInfoE4(i).Centroid;data2=IntInfobw3(indexClose2(i)).Centroid;
%             plot(imgca,data1(:,1),data1(:,2),'r*');plot(imgca,data2(:,1),data2(:,2),'b+');
%             LineDataX=[data1(1);data2(1)];LineDataY=[data1(2);data2(2)];
%             line(LineDataX,LineDataY,'Color','w');
%             end
%         end
%         hold off
    
end

end

areaCh0=table(areaCh0);
area2Ch1=table(area2Ch1);
disTable=table(distanceMaster);
writetable(areaCh0,strcat(savedir,char(nameFolds(mainloop)),'.xlsx'),'Sheet','Area','Range','A1');
writetable(area2Ch1,strcat(savedir,char(nameFolds(mainloop)),'.xlsx'),'Sheet','Area','Range','D1');
writetable(disTable,strcat(savedir,char(nameFolds(mainloop)),'.xlsx'),'Sheet','Distance');
writetable(matrix4,strcat(savedir,char(nameFolds(mainloop)),'.xlsx'),'Sheet','Matrix');
% writetable(indexClose3,strcat(savedir,char(nameFolds(mainloop)),'.xlsx'),'Sheet','index');
% writetable(distance3,strcat(savedir,char(nameFolds(mainloop)),'.xlsx'),'Sheet','Distance');
cd(char(parentfolder));
end
close all