%% Merge all the divided STORM movies in the same layer
%  Plot the point in figures

%% Input
%  *.dat file from STdata5_2S (Redata.m)
%  *.dat file from Insight_to_Matlab.m
%
%    Column 01 : X coordinate in pixel
%    Column 02 : Y coordinate in pixel
%    Column 03 : Intensity (height of the Gaussian peak)
%    Column 04 : Ellipticity (Elp=a-b for Rahul Roy)
%    Column 05 : Sx (Standard deviation of the Gaussian fit in X axis)
%    Column 06 : Sy (Standard deviation of the Gaussian fit in Y axis)
%    Column 07 : Background
%    Column 08 : fitter
%    Column 09 : photoncounts
%    Column 10 : Elliptical ratio (ElpR=a/b)
%    Column 11 : Frame No.

%% Output
%  *.dat file for DriftCorrect.m & Precisionstat.m
%
%    column  1  X in um   MICRO-Meter   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    column  2  Y in um   MICRO-Meter   %%%%% Careful about the unit %%%%%
%    column  3  Z in nm   NANO -Meter   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    column  4  Intensity
%    column  5  Ellipticity
%    column  6  Sx in pixel
%    column  7  Sy in pixel
%    column  8  Background
%    column  9  filterr
%    column 10  photoncounts
%    column 11  Elp ratio
%    column 12  frame no.

%% Modified on Jan 05, 2011 
%  to include the frametime of data acquisition
%  that will be used to correct drift in the sample

%% Rahul: Aug 15, 2011
%  modified Merge_and_plot to now select a few user defined files in the
%  current folder

%% Rahul: Aug 17, 2011
%  modified to print an info file with all the limits and located points

%% Revised by Q.Peter Su Nov 09 2012 
%  1. Use different Z to Ellipticity Calibration Value for different color
%  2. Ellipticity histogram limitation
%  3. Figure background color
%  4. delete some codes

%% Modified by Q.Peter Su on Nov 12-14 2012
%  Different color drift correction -- Red, Green, Rahul Roy, or without Z

%% Modified by Q. Peter Su on Nov 19-21 2012
%  Use strcmp(colorjudge,'n') to judge the string
%  Use strcmp(terminate,'zz') to judge the string

%% Revised by Q. Peter Su on Dec 28 2012
%  Update Six histograms, and input limit due to six histograms
%  1. Intensity(Gaussian peak)   2. Ellipticity=Sx-Sy   
%  3. X position in pixel        4. Y position in pixel
%  5. X Std of width in nm       6. Y Std of width in nm

%% Revised by Q. Peter Su on Dec 28 2012
%  Input the X and Y pixel size in nm

%% Revised by Q.Peter Su on Jan 08 2013
%  Revise the data import for a single file
%  if strcmp(class(FileName),'cell')    numberfiles=length(FileName); 
%  else                                 numberfiles=1;                 end

%% Revised by Q.Peter Su on Apr 16 2013
%    to import multiple files
%    Data = importdata (Afname);
%         Data(:,11) = Data(:,11) + framenumber* (i-1);
%         DataAll = [DataAll; Data];




%function ST02_Merge_and_plotfiles(Xpixel,Ypixel)
function ST02_Merge_and_plotfiles()
clc
clear all;
close all;
DataAll = [];

%% load *.DAT files in the same layer
[FileName,PathName] = uigetfile(... 
    { '*.dat',  'All Data Files (*.dat)'; ...
    %'*RLS*a.dat',' Sheet data set a (*a.dat)'; ...
    %'*RLS*b.dat',' Sheet data set b (*b.dat)'; ...
    '*RLS*.dat','All RLS data(*RLS*.dat)'; ...
    %'*Epi*a.dat',' Epi data set a (RLS*a.dat)'; ...
    %'*Epi*b.dat',' Epi data set b (Epi*b.dat)'; ...
    '*Epi*.dat','All Epi data(*epi*.dat)'}, ...
    'Select *.dat files to use for STORM image generation', 'MultiSelect', 'on');
  %revised on sep 17 2012 by qpsu %revised on sep18 2012 by qpsu

framenumber = input('Input the Frame Number in Each Movie. - ');  
%number of frames per STORM movie
cd(PathName);
if strcmp(class(FileName),'cell')
    numberfiles=length(FileName); 
else     % Revised by Q.Peter Su on Jan 08 2013
    numberfiles=1;
end
Xlim = [];Ylim = [];

%% Figure out
scrsz = get(0,'ScreenSize');
figure('OuterPosition',[1 0.45*scrsz(4) 700 600]);
                      %[left bottom width height]
                      %Marked by Q.Peter Su
whitebg([0 0 0]);
hdl1=gcf;
clf;
figure('OuterPosition',[1 1 500 500]);
whitebg([0 0 0]);
hdl2=gcf;
clf;

%%
sel = 'n'; 
%sel = input('Do you want to analyze only original movie files?   ','s');
inc = 1;
i= 1;
if strcmp(class(FileName),'cell')  % Revised by Q.Peter Su on Jan 08 2013
while (i<=numberfiles)
    Afname=char(FileName(i));
    if lower(sel) == 'y'
        TT = (Afname(end-4)~='a') & (Afname(end-4)~='b');
    else
        TT = (Afname(end-2:end)=='dat');
    end 
    if TT
        fname=Afname(1:end-4);
        str = ['Filename:' fname];
        disp(str);
        TT  = (Afname(end-4)=='a') | (Afname(end-4)=='b')| (Afname(end-4)=='m');
        if (TT)
            filenumber = str2num(Afname(end-8:end-5));
        else
            filenumber = str2num(Afname(end-7:end-4));
        end      
        Data = importdata (Afname);
        Data(:,11) = Data(:,11) + framenumber* (i-1);
        DataAll = [DataAll; Data];
    end  
    i=i+inc;
end
else   % Revised by Q.Peter Su on Jan 08 2013
    Afname=FileName
    Data = importdata (Afname);
    Data(:,11) = Data(:,11);
    DataAll = [DataAll; Data];  
end

%% X and Y Pixel Size
%if nargin<2
   Xpixel      =157;%input('Please input X pixel size (in nm)   - ');
   Ypixel      =157;%input('Please input Y pixel size (in nm)   - ');
%end

%% Preview Figure
figure(hdl1); 
scatter(DataAll(:,1),DataAll(:,2), 1, DataAll(:,4), 'filled'); 
xlabel('X position in pixel');
ylabel('Y position in pixel');
 
terminate = 'zz';

figure(hdl2);
      subplot(3,2,1);
      hist(DataAll(DataAll(:,3)>1000 & DataAll(:,3)<2*10^5,3),100);
      xlabel('Intensity');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,2);
      hist(DataAll(DataAll(:,4)>-2 & DataAll(:,4)<2,4),100);
      xlabel('Ellipticity');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,3);
      hist(DataAll(:,1));
      xlabel('X (pixel)');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,4);
      hist(DataAll(:,2));
      xlabel('Y (pixel)');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,5);
      hist(DataAll(:,5)*Xpixel);          % Standard deviation of the
      xlabel('X width Std (nm)');      % Gaussian fit in X in nm
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,6);           
      hist(DataAll(:,6)*Ypixel);           % Standard deviation of the
      xlabel('Y width Std (nm)');      % Gaussian fit in Y in nm
      ylabel('Occurrence');
      grid on;      zoom on;
  
      mcounter = 0;
%%
while (terminate ~= 'y')   
    if (mcounter == 0)    
        XSel = (0 < DataAll(:,1));
        YSel = (0 < DataAll(:,2));
    end   
%% Intensity limit
      Intensitylim=input('Intensity limits? [IntlimL IntlimH] - ');
%% Elp limit low and high
      Elplim      =input('Elp limits? [ElplimL ElplimH]       - ');
%% judgement
      if strcmp(terminate,'xx')
       Xlim       =input('X limits? [XlimL XlimH]             - ');
       XSel = (Xlim(1)<DataAll(:,1)) & (DataAll(:,1) <Xlim(2));
       mcounter =1;
      end
       
      if strcmp(terminate,'xt')
       Ylim       =input('Y limits? [YlimL YlimH]             - ');
       YSel = (Ylim(1)<DataAll(:,2)) & (DataAll(:,2) <Ylim(2));
       mcounter = 1;
      end
      
       if strcmp(terminate,'zz')
       Xlim      =input('X limits in pixel? [XlimL XlimH]    - ');
       XSel = (Xlim(1)<DataAll(:,1)) & (DataAll(:,1) <Xlim(2));
       
       Ylim      =input('Y limits in pixel? [YlimL YlimH]    - ');
       YSel = (Ylim(1)<DataAll(:,2)) & (DataAll(:,2) <Ylim(2));
       
       Sxlim     =input('Std limit in X in nm [XstdL XstdH]  - ');
       SxSel = (Sxlim(1)<(DataAll(:,5)*Xpixel)) & ...
               ((DataAll(:,5)*Xpixel)<Sxlim(2));
       
       Sylim     =input('Std limit in Y in nm [YstdL YstdH]  - ');
       SySel = (Sylim(1)<(DataAll(:,6)*Ypixel)) & ...
               ((DataAll(:,6)*Ypixel)<Sylim(2));
       
       mcounter = 1;
       
       end
      
       IntSel = (Intensitylim(1)<DataAll(:,3)) &...
                (DataAll(:,3) <Intensitylim(2));
       ElpSel = (Elplim(1)<DataAll(:,4)) & (DataAll(:,4) <Elplim(2));
       DataSel = IntSel & ElpSel & XSel & YSel & SxSel & SySel;
       tSel = find(DataSel);
      
%% Figure 1 Preview dots   
      figure(hdl1);
      scatter(DataAll(tSel,1), DataAll(tSel,2), 1, DataAll(tSel,4), 'filled');
      xlabel('X position in pixel');
      ylabel('Y position in pixel');
      %axis equal
      colorbar;
      
      MeanInt = mean(DataAll(tSel,3));
      StdInt = std(DataAll(tSel,3));
      nInt = 10*ceil((max(DataAll(tSel,3))-min(DataAll(tSel,3)))/std(DataAll(tSel,3)));
    
      MeanElp = mean(DataAll(tSel,4));
      StdElp = std(DataAll(tSel,4));
      nElp = 4*ceil((max(DataAll(tSel,4))-min(DataAll(tSel,4)))/std(DataAll(tSel,4)));
    
      CalX = Xpixel/1000*DataAll(tSel,1);
      MeanX = mean(CalX);
      StdX = std(CalX);
      nX = 8*ceil((max(CalX)-min(CalX))/std(CalX));
    
      CalY = Ypixel/1000*DataAll(tSel,2);
      MeanY = mean(CalY);
      StdY = std(CalY);
      nY = 8*ceil((max(CalY)-min(CalY))/std(CalY));
%% Figure 2 Histogram      
    figure(hdl2);
      subplot(3,2,1);
      hist(DataAll(tSel,3),nInt);
      xlabel('Intensity');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,2);
      hist(DataAll(tSel,4),nElp);
      xlabel('Ellipticity');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,3);
      hist(DataAll(tSel,1),nX);
      xlabel('X (pixel)');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,4);
      hist(DataAll(tSel,2),nY);
      xlabel('Y (pixel)');
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,5);
      hist(DataAll(tSel,5)*Ypixel);      % Standard deviation of the
      xlabel('X width Std (nm)');     % Gaussian fit in X in nm
      ylabel('Occurrence');
      grid on;      zoom on;
      
      subplot(3,2,6);           
      hist(DataAll(tSel,6)*Xpixel);      % Standard deviation of the
      xlabel('Y width Std (nm)');     % Gaussian fit in Y in nm
      ylabel('Occurrence');
      grid on;      zoom on;
      
      totalpoints = length(DataAll);
      selpoints = length(tSel);
      fprintf('Total STORM mapped points is %u out of %u fitted points \n', selpoints, totalpoints);

   figure(hdl1);
   disp('Hi, Are you happy with the data???');
   disp('Please press(y) for yes,(xx) for xlims,');
   disp('(xt) for y limits and (zz) for x & y limits.');
   terminate = input('', 's');
    
end

%% Z Calculation and Information for different color
% Revised on Nov 09 2012 by Q Peter Su
disp('Press r for Red Color, g for Green Color.');
colorjudge=input('rr for Rahul Roy, n without Z value.  ','s');
if strcmp(colorjudge,'r')
    Zcalfile = 'F:\Q. Peter Su\0-Matlab codes for STORM\Calibration Values Red.dat';
end
if strcmp(colorjudge,'g')
    Zcalfile = 'F:\Q. Peter Su\0-Matlab codes for STORM\Calibration Values Green.dat';    
end
if strcmp(colorjudge,'rr')
    Zcalfile = 'F:\Q. Peter Su\0-Matlab codes for STORM\Calibration Values.dat';
end
ZCalibration=importdata(Zcalfile);
p=polyfit(ZCalibration(:,2),ZCalibration(:,1),3);
CalZ = p(1)*(DataAll(tSel,4)).^3+...
       p(2)*(DataAll(tSel,4)).^2+...
       p(3)*(DataAll(tSel,4)).^1+...
       p(4);                          % in nm
if strcmp(colorjudge,'n')
    CalZ=0*DataAll(tSel,4);
end

%% Figure
figure(hdl1);
scatter(CalX, CalY, 1, CalZ, 'filled');
xlabel('X position in um');
ylabel('Y position in um');
%axis equal
whitebg([0 0 0]);
colorbar;

%% all the information for each dots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XYZ = [CalX             ... % column  1  X in um   MICRO-Meter
       CalY             ... % column  2  Y in um   MICRO-Meter
       CalZ             ... % column  3  Z in nm   NANO -Meter
       DataAll(tSel,3)  ... % column  4  Intensity
       DataAll(tSel,4)  ... % column  5  Ellipticity
       DataAll(tSel,5)  ... % column  6  Sx in pixel
       DataAll(tSel,6)  ... % column  7  Sy in pixel
       DataAll(tSel,7)  ... % column  8  Background
       DataAll(tSel,8)  ... % column  9  filterr
       DataAll(tSel,9)  ... % column 10  photoncounts
       DataAll(tSel,10) ... % column 11  Elp ratio
       DataAll(tSel,11)];   % column 12  frame no.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%outputfilename = input('Enter output file name \n', 's');

%cd ..

outputfilename = Afname(1:end-15);
outfilename  = strcat(outputfilename,'ST', '.dat')
save(outfilename, 'XYZ', '-ascii', '-tabs');
figurefilename = strcat(outputfilename,'ST', '.fig')
saveas(gcf,figurefilename);
%% the log file *.inf
%infofile = strcat(PathName, 'MergePlotlog.inf');
infofile = strcat(outputfilename,'ST', '.inf')
fid = fopen(infofile, 'w+');
fprintf(fid, 'Filepath %s\n', PathName);
for i= 1:numberfiles 
    FF =char(FileName(i));
    fprintf(fid, 'Filename %s\n', FF);
end
fprintf(fid, 'ZCal Filename %s\n', Zcalfile);

if isempty(Ylim)
    Ylim(2) = max(DataAll(:,2));
    Ylim(1) = min(DataAll(:,2));
end

if isempty(Xlim)
    Xlim(2) = max(DataAll(:,1));
    Xlim(1) = min(DataAll(:,1));
end

fprintf(fid, ' X limits %s %s\n Y limits %s %s\n Intensity limits %s %s\n Ellipticity limits %s %s\n',...
    num2str(Xlim(1), '%8.2f'), num2str(Xlim(2), '%8.2f'), num2str(Ylim(1), '%8.2f'), num2str(Ylim(2), '%8.2f'), ...
    num2str(Intensitylim(1), '%8.2f'), num2str(Intensitylim(2), '%8.2f'), num2str(Elplim(1), '%8.2f'), num2str(Elplim(2), '%8.2f')); 

fprintf(fid, 'Total STORM mapped points is %u out of %u fitted points\n', selpoints, totalpoints);
fclose(fid);

end