%% Quantifying Puncta Through Intensity Distribution
%Shawn Irgen-Gioro 4/24/2022

% you should only need to change things in this first section

clear all; clc;

%
%this adds some of the functions needed for the script. You will also need
%to have MIJ installed. see http://bigwww.epfl.ch/sage/soft/mij/ for more
%info. Depending on the version of imageJ you have installed, the find
%maxima set parameter might have a different name. Sometimes its called
%"noise", "prominence", etc.

%Tell the program where to find the functions utilized?
folderlocation = 'Z:\Data\20220424\Required Functions';

%adds the functions to your path
addpath(genpath(folderlocation));

%Tell the program where to find the .czi file you want
%make sure you use ' before and after your file paths
dataFolder =' Z:\Manuscripts\Fixation\Raw Imaging Data\Halo - TAF15(IDR)';

%creates a list of files in that folder
dataFiles = dir2(strcat(dataFolder,'\*.czi'));


%Do you want the histogram to have log or linear y axis?
%comment out the one you do not want 

% scale = 'linear'; %add a "%" before this to comment this out. 
scale = 'log'; %add a "%" before this to comment this out. 

%Set the colorbar. What is your output wavelength in nm?
wavelength = 571; %516 is EGFP %571 for JFX549 %583 for DsRed

%For the peak finding algorithm. 
sensitivity = .2;


%%
% for filenum = 1:size(dataFiles,1)
filename = fullfile(dataFiles(filenum).folder, dataFiles(filenum).name);

%load "before"
colormap(emissionWL(wavelength,2^16,.5))

loaddat = ReadImage6D(filename);
data = loaddat{1,1}; %right now the order is [series, time, Z, channel, Y, X], not sure how consistent this will be 
data = permute(data,[6, 5, 3, 2, 4, 1]); %order is now [X, Y, Z, T, C, S], if there's no S it will be a 5D array
metadata = loaddat{1,2};
clear loaddat 
x = 0:metadata.ScaleX:metadata.ScaleX * (metadata.SizeX-1); x = x - median(x);
y = 0:metadata.ScaleY:metadata.ScaleY * (metadata.SizeY-1); y = y - median(y);
z = 0:metadata.ScaleZ:metadata.ScaleZ * (metadata.SizeZ-1); z = z - median(z);

%%  Surface Roughness
clf
maxProj = squeeze(max(data(:,:,:,:,1),[],3));
imagesc(maxProj)
title('Left click to draw outline. Right click to stop drawing')
roi = drawpolyline;
xi = roi.Position(:,1); yi = roi.Position(:,2);

Shape = polyshape(xi,yi);
hold on
plot(Shape)

polymask = poly2mask(xi,yi,size(maxProj,1),size(maxProj,2)); %lets hope that the two images are the same size....

maskedProj = maxProj; maskedProj(~polymask) = 0;
maskvals = maxProj(polymask);

z1masknorm = maskvals./sum(maskvals); 
maskedProj = maskedProj./sum(maskvals);
%% Plot masked Image
subplot(211)
imagesc(x,y,maskedProj)
colorbar
pbaspect([1 1 1])

subplot(212)
numvals = length(unique(z1masknorm));
[a,e] = histcounts(z1masknorm,round(numvals/2));
area(e(1:end-1),a)
ylabel('Number of Pixels')
xlabel('Normalized Intensity')
set(gca,'Yscale',scale)
stdvals = std(z1masknorm);
title(strcat('Normalized Standard Deviation = ',num2str(stdvals,4)))

%% save fig
figname = strrep(dataFiles(filenum).name,'.czi','_histogram.fig');
savepath = fullfile(dataFiles(filenum).folder, figname);
savefig(gcf,savepath)
export_fig(fullfile(dataFiles(filenum).folder,strcat(figname(1:end-4),'.png')),'-native') %as png


%% Find number of puncta using imageJ
ImageJ; %open imageJ
prominence = num2str(sensitivity * max(maskedProj(:)));

MIJ.createImage(maskedProj);
MIJ.run('Find Maxima...',strcat('noise=',prominence,' output=[List]'));
z1peaks = MIJ.getResultsTable;
%
% hold on
% plot(z1peaks(:,1),z1peaks(:,2),'r*')


MIJ.run('Clear Results');
MIJ.closeAllWindows()
MIJ.exit

%
figure(2)
clf; hold on;
imagesc(x,y,maxProj)
p = scatter(x(z1peaks(:,1)+1),y(z1peaks(:,2)+1),'filled',...
    'MarkerFaceColor','red','SizeData',100);
alpha(p,.1)

pbaspect([1 1 1])
axis tight
title(strcat('#_{foci} = ',num2str(size(z1peaks,1))))

%% save fig
figname = strrep(dataFiles(filenum).name,'.czi','_peakfind.fig');
savepath = fullfile(dataFiles(filenum).folder, figname);
savefig(gcf,savepath)
export_fig(fullfile(dataFiles(filenum).folder,strcat(figname(1:end-4),'.png')),'-native') %as png

%% extract FWHM around each point of pair

numslices = 36;
angles = linspace(0, 180, numslices);

L = 10; %how long of a line do we want to draw? currently this is in units of pixels

numpoints = 200;
clim = jet(numslices);
% clf
% imagesc(normz1)
maxFWHM = zeros(1,size(z1peaks,1));
hold on
for ii = 1:size(z1peaks,1)
% for ii = 13

    FWHMs = zeros(1,numslices);
    %find the corresponding pairs to make list in the same order

    for jj = 1:numslices
        %handle the short list first 
        [slice, int_x, int_y] = SliceData(z1peaks(ii,:),angles(jj),L,...
            maxProj,numpoints,true);
        slices(ii,jj,:) = slice;
        int_xs(ii,jj,:) = int_x;
        int_ys(ii,jj,:) = int_y;
        
%         plot(int_x,int_y);
        try
          FWHMs(jj) = fwhm(squeeze(int_xs(ii,1,:)),slice);
        end
    end
    
    [~,FWHM_sizes] = sort(FWHMs,'descend','MissingPlacement','last');
%     FWHM = FWHMs(FWHM_sizes(3)); %picking the 3rd highest peak.
%try picking the value with highest value in histogram
nanidx = isnan(FWHMs); zeroidx = FWHMs==0; removeidx = nanidx|zeroidx;
FWHMs(removeidx) = []; %remove 0's and nans
 
[N, edges, bin] = histcounts(FWHMs,8);
[~, idx] = max(N); 
maxFWHM(ii) = median(FWHMs(bin==idx));

% text(z1peaks(ii,1), z1peaks(ii,2), num2str(mean(diff(x)).*maxFWHM(ii)*1000,3), ...
%         'color','b','FontSize',13,'VerticalAlignment','bottom')%label the circles w. FWHMs
    
%     maxFWHM(ii) = FWHM; % in microns
end
   
%% plot results 

clf
imagesc(maxProj); 

pbaspect([1 1 1])
set(gca,'FontSize',24)

%create scalebar
hold on
plot([size(maxProj,1)*.1; size(maxProj,2)*.1], [size(maxProj,1)*.1; size(maxProj,2)*.1 + (10./mean(diff(x)))], '-w',...
    [size(maxProj,1)*.1 + (10./mean(diff(x))); size(maxProj,2)*.1], [size(maxProj,1)*.1; size(maxProj,2)*.1], '-w', 'LineWidth', 2);
hold off
text(size(maxProj,1)*.1, size(maxProj,2)*.1, '10 \mum', 'color','white','FontSize',18,...
    'VerticalAlignment','bottom') %label scalebar

set(gca, 'Visible', 'off')%turn off labels, etc. around figure
viscircles(z1peaks,maxFWHM./2); %draw circles with width of FWHM
text(z1peaks(:,1), z1peaks(:,2), num2str(mean(diff(x)).*maxFWHM(:)*1000,3), ...
        'color','b','FontSize',13,'VerticalAlignment','bottom')%label the circles w. FWHMs
    

% Get bound fraction from sum image
sumProj = squeeze(sum(data,3));
%     plot(z1peaks(:,1),z1peaks(:,2),'b*')

mask = createCirclesMask(maxProj,z1peaks,maxFWHM./2);
inhub = sum(sumProj(mask));
totfluor = sum(sumProj(polymask));
boundfrac = inhub/totfluor;

text(size(maxProj,1)*.1, size(maxProj,2)*.9, strcat('Bound Percentage = ', num2str(boundfrac*100,2),'%'),...
    'color','white','FontSize',18,...
    'VerticalAlignment','cap')

%% save fig
figname = strrep(dataFiles(filenum).name,'.czi','_boundfrac.fig');
savepath = fullfile(dataFiles(filenum).folder, figname);
savefig(gcf,savepath)
export_fig(fullfile(dataFiles(filenum).folder,strcat(figname(1:end-4),'.png')),'-native') %as png

end
