%% Calculate the STORM image colocalizaion ratio in 3D
% based on Rahul's STcoloc.m
% using filtered precision.mat file output by Yujie's precisionstat.m
%             or    *-Corr.dat file output by DriftCorrect_3D.m
%% Input 
%  1. load *.mat file from precisionstat.m for cluster information
%  2. input distance between the two color in 3D
%  3. Nearest Neighbour.m
%% Output
%  1. plot colocalization figure in 3D with tile
%  2. save two figures in 3D
%% Revised by Q. Peter Su on Oct 25 2012
%  3D STcolocalization calculation
%% Revised by Q. Peter Su on Nov 12 2012
%  coloc definition : X Y within 3 times resolution in X Y (~100 nm)
%                      Z  within 3 times resolution in  Z  (~200 nm)
%% Revised by Q. Peter Su on Nov 14 2012
% add 2D 3D colocalization input choice
%% Revised by Q.Peter Su on Jan 09 2013
%  Input global variable 'distance' and 'dimensionjudge'



function ST09_1_STcoloc_3D(distance,dimensionjudge)
%% clear
clc
%clear all;
close all
%% 1. Data import
% Red dat file
[fileR pathR] = uigetfile({'*prec*.mat','Precision *.mat file';...
                          '*.mat','All mat file';...
                          '*filter*.dat','Filtered dots *.dat file';...
                          '*.dat','All dat files'},...
               'Choose the Filtered-precision Red Cluster *.mat file.',...
               'Multiselect', 'off');
disp(pathR); 
disp('Red DAT File : ');
disp(fileR);
cd(pathR);

% Green dat file
[fileG pathG] = uigetfile({'*prec*.mat','Precision *.mat file';...
                          '*.mat','All mat file';...
                          '*filter*.dat','Filtered dots *.dat file';...
                          '*.dat','All dat files'},...
             'Choose the Filtered-precision Green Cluster *.mat file.',...
             'Multiselect', 'off');
disp('Green DAT File : ');
disp(fileG);
cd(pathG);
red_images = importdata([pathR fileR]);
green_images = importdata([pathG fileG]);
%% 1.1 change matrix
% *.mat file output by Yujie's precisionstat.m
if strcmp(fileR(end-3:end),'.mat');
    % change the red_image to Centers Red
    centersR=[]; 
    centersR(:,1)=red_images(:,1);   % red cluster center X    in um
    centersR(:,2)=red_images(:,3);   % red cluster center Y    in um
    centersR(:,3)=red_images(:,5);   % red cluster center Z    in um
    % change the green_image to Centers Green
    centersG=[];
    centersG(:,1)=green_images(:,1); % green cluster center X  in um 
    centersG(:,2)=green_images(:,3); % green cluster center Y  in um
    centersG(:,3)=green_images(:,5); % green cluster center Z  in um
    
% *.dat file output by Rahul's Merge_and_plotfile.m and Driftcorrct.m
elseif strcmp(fileR(end-3:end),'.dat');
    % change the red_image to Centers Red
    centersR=[]; 
    centersR(:,1)=red_images(:,1);   % red cluster center X    in um
    centersR(:,2)=red_images(:,2);   % red cluster center Y    in um
    centersR(:,3)=red_images(:,3);   % red cluster center Z    in um
    % change the green_image to Centers Green
    centersG=[];
    centersG(:,1)=green_images(:,1); % green cluster center X  in um 
    centersG(:,2)=green_images(:,2); % green cluster center Y  in um
    centersG(:,3)=green_images(:,3); % green cluster center Z  in um
end
    
%% 2. input distance between the two color
if nargin<1
    distance      =100%input('Enter the distance between the two color in nm. - ');
end
distanceRandG=distance/1000;
if nargin<1
    dimensionjudge=2;%input('You Want 2D or 3D Colocalization?  Press 2 or 3 - ');
end

%% 3.1 Nearest neighbour 
%% 3.1.1 Calculate the ratio of (Red cluster coloc to Green)
%        Finds the nearest neighbour by Euclidean distance to 
%        each point in CentersR from CentersG
IDXR2G=nearestneighbour(centersR(:,1:2)',centersG(:,1:2)');
disR2Gxy=sqrt((centersG(IDXR2G,1)-centersR(:,1)).^2 ...
             +(centersG(IDXR2G,2)-centersR(:,2)).^2);
disR2Gz =sqrt((centersG(IDXR2G,3)-centersR(:,3)).^2);
coloc1not0=zeros(length(centersR),1);   % coloc =1 not coloc =0
xR2G=[];yR2G=[];zR2G=[];
% Generate coloc=1,not=0 vector for Red clusters
if dimensionjudge==2
    for i=1:length(centersR);
        if disR2Gxy(i)<distanceRandG;
           coloc1not0(i,1)=1;
           xR2G=[xR2G centersR(i,1)];
           yR2G=[yR2G centersR(i,2)];
           zR2G=[zR2G centersR(i,3)];
        end
    end
elseif dimensionjudge==3
    for i=1:length(centersR);
        if disR2Gxy(i)<distanceRandG && disR2Gz(i)<(2*distanceRandG);
           coloc1not0(i,1)=1;
           xR2G=[xR2G centersR(i,1)];
           yR2G=[yR2G centersR(i,2)];
           zR2G=[zR2G centersR(i,3)];
        end
    end
end    
    centersR(:,4)=coloc1not0(:,1);
    disR2G_ct=sum(centersR(:,4));

%% 3.1.2 Calculate the ratio of (Green cluster coloc to Red)
%        Finds the nearest neighbour by Euclidean distance to
%        each point in CentersG from CentersR
IDXG2R=nearestneighbour(centersG(:,1:2)',centersR(:,1:2)');
disG2Rxy=sqrt((centersR(IDXG2R,1)-centersG(:,1)).^2 ...
             +(centersR(IDXG2R,2)-centersG(:,2)).^2);
disG2Rz =sqrt((centersR(IDXG2R,3)-centersG(:,3)).^2);
coloc1not0=zeros(length(centersG),1);   % coloc =1 not coloc =0
xG2R=[];yG2R=[];zG2R=[];
% Generate coloc=1,not=0 vector for Red clusters
if dimensionjudge==2   
    for i=1:length(centersG);
        if disG2Rxy(i)<distanceRandG
           coloc1not0(i,1)=1;
           xG2R=[xG2R centersG(i,1)];
           yG2R=[yG2R centersG(i,2)];
           zG2R=[zG2R centersG(i,3)];
        end
    end
elseif dimensionjudge==3
    for i=1:length(centersG);
        if disG2Rxy(i)<distanceRandG && disG2Rz(i)<(2*distanceRandG)
           coloc1not0(i,1)=1;
           xG2R=[xG2R centersG(i,1)];
           yG2R=[yG2R centersG(i,2)];
           zG2R=[zG2R centersG(i,3)];
        end
    end
end
    centersG(:,4)=coloc1not0(:,1);
    disG2R_ct=sum(centersG(:,4));

%% 3.2 Display the ratio
OverlapG2R = 100*disG2R_ct/length(centersG);
dispGreenCluster=strcat('Percentage of Green Clusters which are ',...
                       num2str(distance),'nm from a Red Cluster =',...
                       num2str(OverlapG2R),'%');

totalG=strcat('Total Green Clusters = ', num2str(length(centersG)));
disp(totalG);
disp(dispGreenCluster);
disp(' ');

OverlapR2G = 100*disR2G_ct/length(centersR);
dispRedCluster=strcat('Percentage of Red Clusters which are',...
                       num2str(distance),'nm from Green Clusters =',...
                       num2str(OverlapR2G),'%');
totalR=strcat('Total Red Clusters = ', num2str(length(centersR)));
disp(totalR);
disp(dispRedCluster);
disp(' ');
        
%% 4.1 Plot the dots
figure1 = figure('Color',...
    [0.313725501298904 0.313725501298904 0.313725501298904]);
%% 4.1.1 Create axes 1
axes1 = axes('Parent',figure1,'ZColor',[1 1 1],...
    'YDir','reverse','YColor',[1 1 1],...
    'XColor',[1 1 1],...
    'Position',[-0.2 0.09 0.9 0.66],...
    'DataAspectRatio',[1 1 1],...
    'Color',[0 0 0]);
xlim(axes1,[0 40]);
ylim(axes1,[0 40]);
hold(axes1,'all');

xGall=centersG(:,1);
yGall=centersG(:,2);
zGall=centersG(:,3);
% Create plot
plot3(xGall,yGall,zGall,'MarkerFaceColor',[0 1 0],... %[0 1 0] means green
                        'MarkerEdgeColor',[0 1 0],'Marker','.',...
                        'LineStyle','none','Color',[0 1 0]);
% Create plot
plot3(xG2R,yG2R,zG2R,'MarkerFaceColor',[1 1 1],... %[1 1 1] means white
                     'MarkerEdgeColor',[1 1 1],'Marker','*',...
                     'LineStyle','none','Color',[0 1 0]);

title({[fileG(1:end-4)];...
       ['Total Green Clusters = ' num2str(length(green_images))];...
       ['Distance between Two Color =' num2str(distance) 'nm'];...
       ['% of Green to Red =' num2str(OverlapG2R) '%']},...
       'FontWeight','bold','FontSize',10,'FontName','Arial',...
       'Color',[1 1 1]);

%% 4.1.2 Create axes 2
axes2 = axes('Parent',figure1,'ZColor',[1 1 1],...
    'YDir','reverse','YColor',[1 1 1],...
    'XColor',[1 1 1],...
    'Position',[0.3 0.09 0.9 0.66],...
    'DataAspectRatio',[1 1 1],...
    'Color',[0 0 0]);
xlim(axes2,[0 40]);
ylim(axes2,[0 40]);
hold(axes2,'all');

xRall=centersR(:,1);
yRall=centersR(:,2);
zRall=centersR(:,3);
% Create plot
plot3(xRall,yRall,zRall,'MarkerFaceColor',[1 0 0],... %[1 0 0] means red
                       'MarkerEdgeColor',[1 0 0],'Marker','.',...
                       'LineStyle','none','Color',[1 0 0]);
% Create plot
plot3(xR2G,yR2G,zR2G,'MarkerFaceColor',[1 1 1],... %[1 1 1] means white
                     'MarkerEdgeColor',[1 1 1],'Marker','*',...
                     'LineStyle','none','Color',[0 1 0]);
                 
title({[fileR(1:end-4)];...
       ['Total Red Clusters = ' num2str(length(red_images))];...
       ['Distance between Two Color =' num2str(distance) 'nm'];...
       ['% of Red to Green =' num2str(OverlapR2G) '%']},...
       'FontWeight','bold','FontSize',10,'FontName','Arial',...
       'Color',[1 1 1]);
axis equal

if dimensionjudge==2
    saveas(gcf, ['0-STcoloc2D-' fileR(1:21) '-R&G-within'...
              num2str(distance) 'nm-dualcolor.fig'], 'fig');
    saveas(gcf, ['0-STcoloc2D-' fileR(1:21) '-R&G-within'...
              num2str(distance) 'nm-dualcolor.bmp'], 'bmp');
end
if dimensionjudge==3
saveas(gcf, ['0-STcoloc3D-' fileR(1:21) '-R&G-within'...
              num2str(distance) 'nm-dualcolor.fig'], 'fig');
saveas(gcf, ['0-STcoloc3D-' fileR(1:21) '-R&G-within'...
              num2str(distance) 'nm-dualcolor.bmp'], 'bmp');
end

%% 4.2 Plot the Coloc
figure2 = figure('Color',...
    [0.313725501298904 0.313725501298904 0.313725501298904]);
axes3 = axes('Parent',figure2,'ZColor',[1 1 1],...
    'YDir','reverse','YColor',[1 1 1],...
    'XColor',[1 1 1],...
    'Color',[0 0 0]);
box(axes3,'on');
hold(axes3,'all');

% Create plot
plot3(xG2R,yG2R,zG2R,'Marker','.','LineStyle','none','Color',[0 1 0]);
% [0 1 0] means Green

% Create plot
plot3(xR2G,yR2G,zR2G,'Marker','.','LineStyle','none','Color',[1 0 0]);
% [1 0 0] means Red
if dimensionjudge==2
    title({['0-STcoloc2D-',fileR(1:11),'-R&G-within',...
        num2str(distance),'nm-merge.fig']},...
        'FontWeight','bold','FontSize',14,'FontName','Arial',...
        'Color',[1 1 1]);
    saveas(gcf, ['0-STcoloc2D-',fileR(1:21),'-R&G-within',...
        num2str(distance),'nm-merge.fig'], 'fig');
    saveas(gcf, ['0-STcoloc2D-',fileR(1:21),'-R&G-within',...
        num2str(distance),'nm-merge.bmp'], 'bmp');
end

if dimensionjudge==3
    title({['0-STcoloc3D-',fileR(1:11),'-R&G-within',...
        num2str(distance),'nm-merge.fig']},...
        'FontWeight','bold','FontSize',14,'FontName','Arial',...
        'Color',[1 1 1]);
    saveas(gcf, ['0-STcoloc3D-',fileR(1:21),'-R&G-within',...
        num2str(distance),'nm-merge.fig'], 'fig');
    saveas(gcf, ['0-STcoloc3D-',fileR(1:21),'-R&G-within',...
        num2str(distance),'nm-merge.bmp'], 'bmp');
end

end
