%%% Copyright: T. Hirayama., K. Takabe, Y. Yawata, 2020 %%%
%% Script2D
%  This script and parameters have optimized to analyze "putida2440" data successfully.

[pathname] = uigetdir('','Select a directory');
cd(pathname);
mkdir(pathname,'/Signature/');
savepath = [pathname,'/Signature/'];

CRMlist = dir(['*CRM*','.tif']);
Z = length(CRMlist);

%% Input Parameters
prompt = {'Threshold for image binarization [0-1]', 'Upper threshold for a cell region (in pixels)','Lower threshold for a cell region',...
    'Number of detectors'};
dlg_title = 'Parameters';
num_lines = [1 70];
defaultans = {'0.45','200','10',...
    '32'};
Parameters = inputdlg(prompt,dlg_title,num_lines,defaultans);

BW_th = str2num(Parameters{1});
MAX = str2num(Parameters{2});
MIN = str2num(Parameters{3});
EM_ch = str2num(Parameters{4});

prompt = {'Number of excitation wavelengths used for image acquisition'};
dlg_title = 'Input';
num_lines = [1 50];
defaultans = {'6'};
NumInput = inputdlg(prompt,dlg_title,num_lines,defaultans);

NumEX = str2num(NumInput{1});

Pex = [];
Dex = [];
for i = 1:NumEX
    TempPex = ['Excitation No.',num2str(i)];
    TempDex = [''];
    Pex = [Pex,{TempPex}];
    Dex = [Dex,{TempDex}];
end    

prompt = Pex;
dlg_title = 'Input wavelengths (in a sequence from shortest to longest wavelength)';
num_lines = [1 130];
defaultans = Dex;
EXlist = inputdlg(prompt,dlg_title,num_lines,defaultans);

EX_wl = EXlist;
%% Image processing

[H W] = size(imread(CRMlist(1).name));
CRMstack = uint16(zeros(H,W,Z));
for z=1:Z
    CRMstack(:,:,z)= imread(CRMlist(z).name);
end
IM = max(CRMstack,[],3);

figure(1); 
IM = imadjust(IM);
figure(1); subplot(2,3,1); imshow(IM); title('Original image (imadjust)');

IM_flat = imflatfield(IM,3);
figure(1); subplot(2,3,2); imshow(IM_flat); title('imflatfield');

IM_comp = imadjust(imcomplement(IM_flat));
figure(1); subplot(2,3,3); imshow(IM_comp); title('imcomplement');

IM_bin = imbinarize(IM_comp, BW_th);
figure(1); subplot(2,3,4); imshow(IM_bin); title('imbinarize');

IM_area = bwareafilt(IM_bin,[MIN MAX], 8);
IM_cb = imclearborder(IM_area);
figure(1); subplot(2,3,5); imshow(IM_cb); title('bwareafilt & imclearborder');

IM_ovl = imoverlay(IM, IM_cb);
figure(1); subplot(2,3,6); imshow(IM_ovl); title('Segmented image');

IM_final = IM_cb;
CC = bwconncomp(IM_final, 4);
NUM = CC.NumObjects; %the number of object

%% Extraction of innate fluorescence signature based on segmented regions

EX = length(EX_wl); %The number of excitaions
EM = max(EM_ch); %Max number of detector channels

uiwait(msgbox(('Draw a rectangle to select a background region')));
figure('Name', 'Draw a rectangle to select a background region');
[BG, BG_rect] = imcrop(IM);
close

FL_list = zeros(EX,EM,NUM);
for k = 1:EX
    for l = 1:EM

        FL = dir(['*',EX_wl{k},'c',num2str(l, '%02i'),'z*','.tif']);
        FL_stack=uint16(zeros(H,W,Z));
        for z=1:Z
            FL_stack (:,:,z) = imread (FL(z).name);
        end
        FL_max = max(FL_stack,[],3);
        

        BG_mask = imcrop(FL_max,BG_rect);
        BG_mean = mean(BG_mask(:));
        

        for n = 1:NUM
            mask=zeros(H,W);
            mask(CC.PixelIdxList{1,n})=1;
            ROI = FL_max.*uint16(mask);
            FL_cell = mean(FL_max(mask>0)) - BG_mean;
            
            FL_list(k,l,n) = FL_cell;
        end
    end
end

%% Export signatures as PNG format

FL_list_norm = zeros(EX,EM,NUM);
for n = 1:NUM
    int_lim = stretchlim(uint16(FL_list(:,:,n)),5e-4);
    FL_list_norm(:,:,n) = imadjust(uint16(FL_list(:,:,n)),[int_lim(1) int_lim(2)],[0 1]);
    imwrite(uint16(FL_list_norm(:,:,n)),[savepath,'Signature',num2str(n,'%02d'),'.png']);
end


%% Display cell-map and heatmap

% heatmap
Signature = uint16(zeros(EX*EM,NUM));
for n = 1:NUM
    Signature(:,n) = transpose(reshape(FL_list_norm(:,:,n).',[192,1]));
end
figure(2); 
Heatmap(Signature, pathname);

% cell-map
figure(3); 
imshow(IM_ovl)
Cent = regionprops(CC, 'Centroid');
for l = 1:NUM
    x = Cent(l).Centroid(1);
    y = Cent(l).Centroid(2);
    text(x, y,(num2str(l)));
end

%% Funstion for display of heatmap
function Heatmap(A,path)
    MaxV = max(max(A));
    h = heatmap(A,'Colormap',parula,'ColorLimits',[0 MaxV]);
    h.XLabel = 'Cells';
    h.YLabel = 'Detector Number';
    h.FontSize = 8;
    h.colorbar;
    saveas(gcf,[path,'/Heatmap.png'])
end