%% Write by Yujie Sun & Q. Peter Su on May 13-15 2013
%  for transfer *list.bin file to matrix with can be called in Matlab
%
%h
%  Part O  - Calculate the Z Calibration Value
%            X width & Y width output by Insight3, done by Yujie Sun
%            using box size (ROI) of 13.
%
%  Part I  - Preview ~2000 Frame Molecule (Intensity to Width)
%            This part is optional when you analyze a new sample.
%
%  Part II - Convert Mol.List *.bin to *.dat file
%            Suitble for ST05_DriftCorrect_3D.m
%                        ST07_precisionstat_3D.m
%       II-1-Read the *list.bin file, X Y Ax W Int BG & frame no. 
%            with ReadMasterMoleculeList.m from Xiaowei Zhuang's Lab.
%       II-2-Calculate the Z position with Ax & W according to 20130412
%            Z Calibration Value for Red(Alexa647) & Green(Atto565) Channel
%       II-3-For precise image, delete the dots which z position is out of
%            [ -300 nm , 250 nm ], then apply the coefficient number 0.72.
%       II-4-Delete Dots in BF frames
%       II-5-Output the dots/localizations as *list_UnCorr.dat
%                                          or *list_AutoCorr.dat
%
%  PartIII - DriftCorrect for single-color if necessary
%            DriftCorrect for two-color    if necessary
%
%  Part IV - Delete UnExpected Points if necessary
%
%  Part V  - precisionstat_3D
%
%  Part VI - ST08_1_fastmap(Xpixel, Ypixel,binsize);
%
%  sunmoon286@gmail.com
%  petersq001@gmail.com


%% Input
%      The molecule list *.bin file auto-output by insight3
%      after analyzing dots, choose 'KeepMList'
%                                   'Fix azimuth angle'
%       Column  1 : Cas### - Number of All the HR molecules
%       Column  2 : X      - X Coordinate in Pixel
%       Column  3 : Y      - Y Coordinate in Pixel
%       Column  4 : Xc     - X in Pixel after Auto Drift Correction
%       Column  5 : Yc     - Y in Pixel after Auto Drift Correction
%       Column  6 : Height - Intensity/Height of Gaussian Fitting
%       Column  7 : Area   - Area of the Bottom of Gaussian Fitting
%       Column  8 : Width  - sqrt(Wx * Wy)
%       Column  9 : Phi    - 
%       Column 10 : Ax     - Wx/Wy
%       Column 11 : BG     - Background
%       Column 12 : I
%       Column 13 : Frame  - Frame No.
%       Column 14 : Length - Molecule length in Frame
%       Column 15 : Link
%       Column 16 : Valid
%       Column 17 : Z      - Z in nm
%       Column 18 : Zc     - Z in nm after Auto Correction

%% Output
%    *.dat file suitable for DriftCorrect.m & precisionstat.m
%       Column  1 : X in Micrometer 
%       Column  2 : Y in Micrometer
%       Column  3 : Z in Nanometer
%       Column  4 : Intensity/Height
%       Column  5 : Ellipticity = Wx - Wy in nm
%       Column  6 : X Width in pixel
%       Column  7 : Y Width in pixel
%       Column  8 : Background
%       Column  9 : Dots length in fames (06th July 2017 by Q. Peter SU)
%       Column 10 : 
%       Column 11 : Ellipticity ratio = Wx / Wy
%       Column 12 : Frame No.

%% Revised by Q.Peter Su on May 20 2013
%  Connect Part II III IV V with 'outDATfilename'
%  Avoid select the *.dat file repeatedly

%% Revised by Q.Peter Su on May 24 2013
%  1. Show the AutoDrift Corr Curve generated by Insight3
%     DriftFile=[Afname(1:end-9), '_drift.txt'];
%     while exist(DriftFile)~=2
%         disp('Cannot Find the Drift.txt File in This Folder.');
%         judge=input('Do STORM Auto Drift Corr in Insight3, then Press Enter','s');
%     end
%     disp(DriftFile);Drift=importdata(DriftFile);
%     close all
%     figure1=figure;
%  2. Diary log file
%     diaryFile=strcat(Afname(1:end-4),datestr(now, 1),'-diary.txt');
%     diary(diaryFile);

%% Revised by Q.Peter Su on May 24 2013
%  Hist 1 - figure of Height to Width 
%  Hist 2 - hist of intensity
%  Hist 3 - hist of Ellipticity
%  Hist 4 - hist of width
%  Hist 5 - hist of dots counts to frame
%  Hist 6 - hist of dots length

%% Revised by Q.Peter Su on Nov 13 2013
%  display and calculate error for some pretty big file

%% Revised by Q.Peter Su on Nov 18 2013
%   binsize=5;
%   resolution=25;
%   ST08_1_fastmap(outDATfilename,Xpixel,Ypixel,binsize,resolution);

%% Revised by Q.Peter Su on Nov 20 2013
%  function ST21_Insight_to_Matlab_V3_6(Xpixel,Ypixel)
%  input Z coordinate from Insight3

%% Revised by Q.Peter Su on Nov 20 2013
%  avoid wrong color info...
%     judge_color='z';
%     color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
%     while ~color_continue
%         judge_color=input('Type -r- for Red Channel, -g- for Green, -b- for Blue - ','s');
%         color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
%     end

%% Revised by Q.Peter Su on Nov 21 2013
%  plot3k function is much much faster than scatter3...
%
%     figure(hdl1)
%     if Dimension=='2'
%         plot3k([XYZt(1:MolLim,1)*1000,XYZt(1:MolLim,2)*1000,XYZt(1:MolLim,3)],...
%                 XYZt(1:MolLim,3),[-100 100],{'o',1});
%     elseif Dimension=='3'
%         plot3k([XYZt(1:MolLim,1)*1000,XYZt(1:MolLim,2)*1000,XYZt(1:MolLim,3)],...
%             XYZt(1:MolLim,3),Zrange,{'o',1});
%     end
%     %scatter3(XYZt(1:MolLim,1),XYZt(1:MolLim,2),XYZt(1:MolLim,3)/1000,1,XYZt(1:MolLim,3)/1000,'filled');
%     xlabel('X Coordinate (nm)');ylabel('Y Coordinate (nm)');
%     zlabel('Z Coordinate (nm)');
%     axis equal
%     view([0 -90]);

%% Revised by Q.Peter Su on Nov 21 2013 
%    diaryFile=strcat(Afname(1:end-4),'-', datestr(now,'mmm_dd_yyyy-HH_MM_SS'),'-diary.txt');
%    diary(diaryFile);

%% Revised by Q.Peter Su on Nov 22 2013
%  add some time point display
%     disp('Date and Time Now ')
%     disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
%     disp('------------------------------------------------------------------')

%% Revised by Q.Peter Su on Dec 30 2013
%  fix some problems in histogram
%  

%% Revised by Q.Peter Su on Feb 19 2014
%  optimize the DELETE BF parts
%
%     for j=2:(BFIndexNo-1)
%         indextemp=[];
%         for i=1:rangeDel
%             indextemp=[indextemp BFIndex(BFIndexCycle(j-1))-i];
%             indextemp=[indextemp BFIndex(BFIndexCycle(j)-1)+i];
%         end
%         BFIndexDel=[BFIndexDel indextemp];
%     end
%
%% Revised by Q.Peter Su on May 20 2014
%  judge_input23=input('Type in 2 for input 2D data only, 3 for input 3D data - ');

%% Revised by Q.Peter SU on Nov 10 2015
%  some bugs for 3D data are modified
%  renmae to ST21_Insight_to_Matlab_V3_7(Xpixel,Ypixel)

%% Revised by Q.Peter SU on 06th July 2017
%  for Wei REN's UCNPs sample, the STORM data need read channel information
%  from Insight3, actually, to get rid of the channel 0 or 9???, which are the
%  repeatative localized dots, just delete the dots continuously appeared
%  in more than 3 (or more frames)





%%
function ST21_Insight_to_Matlab_V3_8(Xpixel,Ypixel)
clc
close all;

disp('*******************************************************************');
disp('Thanks for using ST21_Insight_to_Matlab_V3_8, written by Q.Peter SU');
disp('This is V3_8, which is updated on 06th July 2017 in HRI,USYD,SYD,AU');
disp('*******************************************************************');
%% Part O - Z Calibration Value
% Z calibration value for Alexa647 calculated by Sunmoon
%   Data from Apr 12 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% As the Objective move downward, the distance between %%%%
%%%% dye molecules and focal plan changes from negative Z %%%%
%%%% value to positive Z value, with the X width increase %%%%
%%%% and Y width decrease (Gaussian fit of the fluo spot) %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% InsightMeanX_R=[362.4348798	349.4594216	326.2824944	332.0297645	314.7796171...
%     308.4924923	314.2696309	310.471812	296.741819	295.3795521	285.1561133...
%     297.1898356	303.987819	317.8996355	335.446684	350.0655242	369.2775761...
%     408.4805658	424.5484789	489.0491904	464.9190438	490.3707009	517.2138633...
%     535.203943	584.6960731	672.315563	661.2182837	682.4663547];
InsightMeanX_R=[310.471812	296.741819	295.3795521	285.1561133...
    297.1898356	303.987819	317.8996355	335.446684	350.0655242	369.2775761...
    408.4805658	424.5484789	489.0491904	464.9190438	490.3707009	517.2138633...
    535.203943	584.6960731	672.315563	661.2182837	682.4663547];
InsightMeanX_R=InsightMeanX_R';


% InsightMeanY_R=[590.7396688	561.6270423	589.870637	585.4556092	617.4934034...
%     584.3981228	606.6216049	599.625734	574.0457175	547.9374236	544.8792301...
%     519.1246196	507.8802412	508.5791658	486.9015503	472.0951934	432.4498765...
%     430.4685334	395.422518	375.8295222	372.5964183	362.2698378	342.5236975...
%     361.5479024	360.1438998	335.3342759	362.3436225	382.6974581];
InsightMeanY_R=[599.625734	574.0457175	547.9374236	544.8792301...
    519.1246196	507.8802412	508.5791658	486.9015503	472.0951934	432.4498765...
    430.4685334	395.422518	375.8295222	372.5964183	362.2698378	342.5236975...
    361.5479024	360.1438998	335.3342759	362.3436225	382.6974581];
InsightMeanY_R=InsightMeanY_R';

InsightMeanZ_R=[-300:25:200];
InsightMeanZ_R=InsightMeanZ_R';
Ellip_R=InsightMeanX_R-InsightMeanY_R;

p_R=polyfit(Ellip_R,InsightMeanZ_R,3);
Ellip_R_plot=[-400:1:400]';
Z_cal_R=p_R(1) * Ellip_R_plot.^3 + ...
        p_R(2) * Ellip_R_plot.^2 + ...
        p_R(3) * Ellip_R_plot.^1 + ...
        p_R(4) ;

figure1 = figure;
whitebg([0 0 0])
subplot(1,2,1)
plot(Ellip_R,InsightMeanZ_R,'ro');
xlabel('Red Channel - Ellipticity=Wx-Wy (nm)');
ylabel('Red Channel - Z position relative to focal plane');
hold on
plot(Ellip_R_plot,Z_cal_R,'r')
title('Red Channel Ellip to Z');
xlim([-600 600]);
ylim([-600 600]);
hold on
plot(0,0,'MarkerFaceColor',[1 1 1],...
    'MarkerEdgeColor',[1 1 1],...
    'MarkerSize',10,...
    'Marker','*',...
    'LineStyle','none',...
    'Color',[1 1 1]);
hold off

% Z calibration value for Atto565 calculated by Sunmoon
%   Data from Apr 12 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% As the Objective move downward, the distance between %%%%
%%%% dye molecules and focal plan changes from negative Z %%%%
%%%% value to positive Z value, with the X width increase %%%%
%%%% and Y width decrease (Gaussian fit of the fluo spot) %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% InsightMeanX_G=[370.5023146	269.1347544	293.8403971	276.1902599	281.1003878...
%     293.989937	274.916144	275.2354091	279.575947	277.4642244	269.4089540...
%     287.1993199	291.7407696	293.8339509	301.3398031	316.7378146	310.2278026...
%     343.5733048	349.9468974	359.7114469	370.8807685	410.7851087	407.5608816...
%     440.5770014	462.3683643	488.2584649	521.0760543	542.9720148	587.7812541...
%     611.9224658	628.720666];
InsightMeanX_G=[275.2354091	279.575947	277.4642244	269.4089540...
    287.1993199	291.7407696	293.8339509	301.3398031	316.7378146	310.2278026...
    343.5733048	349.9468974	359.7114469	370.8807685	410.7851087	407.5608816...
    440.5770014	462.3683643	488.2584649	521.0760543	542.9720148	587.7812541...
    611.9224658	628.720666];
InsightMeanX_G=InsightMeanX_G';

% InsightMeanY_G=[487.5577709	548.829634	581.505102	566.8532811	561.4851037...
%     568.559944	547.5692956	557.3492212	532.7592617	530.5241671	518.9658021...
%     521.8337922	546.4427402	505.5871527	480.0016993	488.3404226	424.4859898...
%     424.7169258	413.156248	404.6029457	388.4410475	380.2735909	361.9215875...
%     357.5210161	336.6801145	335.2794242	327.2307519	335.3867448	349.3669540...
%     355.0399337	354.3676616];
InsightMeanY_G=[557.3492212	532.7592617	530.5241671	518.9658021...
    521.8337922	546.4427402	505.5871527	480.0016993	488.3404226	424.4859898...
    424.7169258	413.156248	404.6029457	388.4410475	380.2735909	361.9215875...
    357.5210161	336.6801145	335.2794242	327.2307519	335.3867448	349.3669540...
    355.0399337	354.3676616];
InsightMeanY_G=InsightMeanY_G';

InsightMeanZ_G=[-300:25:275];
InsightMeanZ_G=InsightMeanZ_G';
Ellip_G=InsightMeanX_G-InsightMeanY_G;

p_G=polyfit(Ellip_G,InsightMeanZ_G,3);
Ellip_G_plot=[-400:1:400]';
Z_cal_G=p_G(1) * Ellip_G_plot.^3 + ...
        p_G(2) * Ellip_G_plot.^2 + ...
        p_G(3) * Ellip_G_plot.^1 + ...
        p_G(4) ;

subplot(1,2,2)
plot(Ellip_G,InsightMeanZ_G,'go');
xlabel('Green Channel - Ellipticity=Wx-Wy (nm)');
ylabel('Green Channel - Z position relative to focal plane');
hold on
plot(Ellip_G_plot,Z_cal_G,'g')
title('Green Channel Ellip to Z')
xlim([-600 600]);
ylim([-600 600]);
hold on
plot(0,0,'MarkerFaceColor',[1 1 1],...
    'MarkerEdgeColor',[1 1 1],...
    'MarkerSize',10,...
    'Marker','*',...
    'LineStyle','none',...
    'Color',[1 1 1]);
hold off


%% Part I - Preview ~2000 Frame Molecule
disp('-------------------------------------------------------------')
disp('If You are Just Starting the Analysis of a New Sample, Type in -new-');
judge_pre=input('Or, If You Already Got the Right Patameters for Insight3, press Enter - ','s');
if  strcmp(judge_pre,'new')
disp('-------------------------------------------------------------');
disp('Analyze about 10k Dots(1k~2k frames) in Insight3 with Default Settings.')
disp('1. View -> Preferences -> ');
disp('   Parallel threads - 8 - ');
disp('2. Analysis -> Parameters -> ');
disp('   Minimum height - 1000;  Maximum height  - 20000      ');
disp('   Fit ROI        - 7   ;  CCD pixel size  - 157 or 105 ');
disp('   CCD baseline   - ~180 ; Photons/counts  - 0.41       ');
disp('   Model Botton 0 -     Simple Gaussian for 2D Analysis ');
disp('   Model Botton 1 - Elliptical Gaussian for 3D Analysis ');
disp('   Minimum width  - 100 ;  Maximum width   - 800        ');
disp('   Initial width  - 300 ;  Max axial ratio - 2.5        ');
disp('   * Choose the Fix azimuth angle ');
disp('   * Disable Elliptical PSF for 2D Analysis');
disp('   * Choose  Elliptical PSF for 3D Analysis');
disp('3. STORM -> STORM analysis -> ');
disp('   Input Start frame and End frame ; 0 for continuous activation');
disp('   * Choose Keep frame MList and Drift Corr');

disp('-------------------------------------------------------------');
judge=input('Press Enter Key to Continue. - ','s');
close all

% Import molecule *list.bin file with Zhuang's ReadMasterMoleculeList();
[preMList, memoryMap1] = ReadMasterMoleculeList();
preX      = preMList.x;
preY      = preMList.y;
preHeight = preMList.h;              % Counts : Height of Gaussian Fit
preWidth  = preMList.w;              % nm     : SQRT(Wx * Wy)
preAx     = preMList.ax;             % A.U    : Wx / Wy
preWidthX = preWidth .* sqrt(preAx); % nm     : Wx = sqrt(Wx * Wy) * sqrt(Wx / Wy)
preWidthY = preWidth ./ sqrt(preAx); % nm     : Wy = sqrt(Wx * Wy) / sqrt(Wx / Wy)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% In Insight3 the Wx and Wy are reversed, so, change them %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
temp      = preWidthX;preWidthX = preWidthY;preWidthY = temp; %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
preEllip  = preWidthX - preWidthY;   % nm     : Ellipticity = Wx - Wy
preFrame  = preMList.frame;  preFrame  = double(preFrame);
preLength = preMList.length; preLength = double(preLength);

Afname1=memoryMap1.Filename;
Afname='';
NameLength=length(Afname1);
i=1;
while Afname1(NameLength-i+1)~='\'
    Afname=[Afname1(NameLength-i+1) Afname];
    i=i+1;
end

close all
scrsz = get(0,'ScreenSize');
figure('OuterPosition',[1 0.5*scrsz(4) 0.75*scrsz(4) 0.5*scrsz(4)]);
                        %[left bottom width height]
whitebg([0 0 0]);hdl1=gcf;clf;
figure('OuterPosition',[1 1  0.75*scrsz(4) 0.5*scrsz(4)]);
whitebg([0 0 0]);hdl2=gcf;clf;

figure(hdl1)
scatter(preX,preY,1,preEllip,'filled');
xlabel('X Coordiante in Pixel');ylabel('Y Coordinate in Pixel');axis equal
title('Pre-view the HR Image, Colorbar Means Ellipticity');
set(gca,'YDir','reverse');
colorbar;

figure(hdl2)
whitebg([0 0 0]);

subplot(2,3,1);
%scatter(preWidth,preHeight,2,'filled');
plot(preWidth,preHeight,'MarkerSize',1,'Marker','.',...
    'LineStyle','none','Color',[1 1 0]);
hold on
plot([300 250 650 300],[16000 1000 1000 16000],'Marker','o','LineWidth',2,...
    'LineStyle','--','Color',[1 0 0]);
hold off
xlim([100 800]);ylim([0 20000]);
xlabel('Width of Gaussian Fit');ylabel('Height of the Gaussian Fit');
title('Find a Triangle ROI in this Figure');grid on;zoom on;

subplot(2,3,2);
hist(preHeight,100);xlim([0 20000]);xlabel('Intensity (A.U.)');ylabel('Occurrence');view([90 -90]);
title('Intensity of Every Dots');grid on;zoom on;

subplot(2,3,3)
hist(preEllip,100);xlim([-500 500]);xlabel('Ellipticity(nm)');
ylabel('Occurrence');title('Ellipticity of Every Dots');grid on;zoom on;

subplot(2,3,4);
hist(preWidth,100);xlim([100 800]);xlabel('Width of Gaussian Fit');
ylabel('Occurrence');title('Width of Every Dots');grid on;zoom on;

subplot(2,3,5)
hist(preFrame,max(preFrame));xlabel('Frame No.');ylabel('Dots Occurrence');
title('Dots Counts in Every Frame');grid on;zoom on;

subplot(2,3,6)
hist(preLength,max(preLength));xlabel('Dots Length');ylabel('Occurrence');
title('Dots Length');grid on;zoom on;


disp('-------------------------------------------------------------');
disp('Before the Final Analysis, Re-name the former *list.bin file to *list_pre.bin.')
disp('Now, Use the Width Range to Finish the Analysis of All the Frames in Insight3.')
judge=input('Press Enter Key to Continue. - ','s');
else
    disp('-------------------------------------------------------------');
end

%% Part II - Convert Mol.List *list.bin to *.dat file
% Import molecule *list.bin file with Zhuang's ReadMasterMoleculeList();
if isempty(judge_pre)
    disp('-------------------------------------------------------------');
    disp('Now, Use the Right Width Range and Intensity Threshold to');
    disp('Finish the Analysis of All the Frames in Insight3.');
    judge=input('Press Enter Key to Continue. - ','s');
end
[dataMList, datamemoryMap] = ReadMasterMoleculeList();
Afname1=datamemoryMap.Filename;
Afname='';
NameLength=length(Afname1);
i=1;
while Afname1(NameLength-i+1)~='\'
    Afname=[Afname1(NameLength-i+1) Afname];
    i=i+1;
end

diaryFile=strcat(Afname(1:end-4),'-diary','-', datestr(now,'mmm_dd_yyyy-HH_MM_SS'),'.txt');
diary(diaryFile);
disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')



 if nargin<1 
%     if Xpixel==0 % for GUI input not a number
        Xpixel      =input('Please input X pixel size (in nm)   - ');
        Ypixel      =input('Please input Y pixel size (in nm)   - ');
%     end
 else
    Xpixelstrout=strcat('X Pixel Size = ',num2str(Xpixel));
    disp(Xpixelstrout);
    Ypixelstrout=strcat('Y Pixel Size = ',num2str(Ypixel));
    disp(Ypixelstrout);
 end

disp('-------------------------------------------------------------');
disp('You Wanna Import the Auto-Drift-Corrected HR Molecules from Insight3?');
judge_file=input('Press c for Corrected Data, u for Un-Corrected Data. - ','s');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Matrix XYZt for all dots/localizations %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% XYZt=[]; % empty matrix got a pretty low speed
XYZt=zeros([size(dataMList.x,1),12]);



    %%%%%%%%%%%% Column 1 : X coordinate in um %%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%% Column 2 : Y coordinate in um %%%%%%%%%%%%%%%%%%
if judge_file=='u'
    XYZt(:,1) = dataMList.x * Xpixel/1000;   % X coordinate (Micrometer)
    XYZt(:,2) = dataMList.y * Ypixel/1000;   % Y coordinate (Micrometer)
elseif judge_file=='c'
    % in the bin file, column of xc is not the corrected x coordinate 
    DriftFile=[Afname(1:end-9), '_drift.txt'];
    while exist(DriftFile)~=2
        disp('Cannot Find the Drift.txt File in This Folder.');
        judge=input('Do STORM Auto Drift Corr in Insight3, then Press Enter','s');
    end
    disp(DriftFile);Drift=importdata(DriftFile);
    % revised by Q.Peter Su on Nov 13 2013
    % display error for some pretty long movie
    FrameNo=size(Drift,1);
    if FrameNo<=1000
        DriftDisplay=Drift;
    else
        DriftDisplay=zeros([1000 4]);
        %DriftDisplayIndex=zeros([1000 1]);
        for i=1:1000
            DriftDisplayIndex=floor(i/1000*FrameNo);
            DriftDisplay(i,:)=Drift(DriftDisplayIndex,:);
        end
    end
 
    close all
    figure1=figure;
    axes1 = axes('Parent',figure1);
    scatter(DriftDisplay(:,2)*Xpixel/1000,DriftDisplay(:,3)*Ypixel/1000,3,DriftDisplay(:,1),'filled');
    title({[DriftFile];['Auto Drift Curve Generated by Insight3(sellected 1000 dots)']});
    xlabel('X Drift in um');ylabel('Y Drift in um');
    set(gca,'YDir','reverse');
    colorbar('peer',axes1);
    axis equal
    box on
    DriftFileOut1 = strcat(DriftFile(1:end-4),'.fig');
    saveas(gcf,DriftFileOut1,'fig');
    DriftFileOut2 = strcat(DriftFile(1:end-4),'.bmp');
    saveas(gcf,DriftFileOut2,'bmp');

    
    disp('-------------!!!!!!!!!!! Caution  !!!!!!!!!!!!-----------------')
    disp('-- Insight3 WILL NOT Save Corrected Coordinate Automatically --')
    disp('-- In Insight3, Sellect STORM -> Automatic STORM Drift Corr ---');
    disp('---------------------------------------------------------------')

    judge_drift=input('File -> Save Molecule List. then Back to Matlab, Press Enter - ','s');
    disp('---------------------------------------------------------------')
    disp('Drift Correcting ... ...');
    XYZt(:,1) = dataMList.xc;                 % X coordinate (pixel)
    XYZt(:,2) = dataMList.yc;                 % Y coordinate (pixel)  
    XYZt(:,1)=XYZt(:,1) * Xpixel/1000;       % X coordinate (Micrometer)
    XYZt(:,2)=XYZt(:,2) * Ypixel/1000;       % Y coordinate (Micrometer)
end



disp('Data Reading ... ...');
    %%%%%%%%%%%% Column 3 : Z coordinate in nm %%%%%%%%%%%%%%%%%%
    %%%%%%%%%%% will Calculate after colorjudge %%%%%%%%%%%%%%%%%
XYZt(:,3) = 0;                               % Z coordinate (Nanometer)

    %%%%%%%%%%%%  Column 4 : Intensity/Height  %%%%%%%%%%%%%%%%%%
XYZt(:,4) = dataMList.h;                     % Height of Gaussian Peak (Counts)
dataHeight= dataMList.h; 

    %%%%%%%%%%%% Column 5 : Ellipticity=Wx-Wy (nm)%%%%%%%%%%%%%%%
    dataWidth   = dataMList.w;               % SQRT(Wx * Wy) (Nanometer)
    dataAx      = dataMList.ax;              % Wx / Wy (Counts)
    dataWidthX  = dataWidth .* sqrt(dataAx); % Wx = sqrt(Wx * Wy) * sqrt(Wx / Wy)
    dataWidthY  = dataWidth ./ sqrt(dataAx); % Wy = sqrt(Wx * Wy) / sqrt(Wx / Wy)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% In Insight3 the Wx and Wy are reversed, so, change them %%%%%%%%
    temp      = dataWidthX;dataWidthX = dataWidthY;dataWidthY = temp; %%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XYZt(:,5) = dataWidthX - dataWidthY;         % Ellipticity = Wx - Wy (Nanometer)
dataEllip = dataWidthX - dataWidthY;

    %%%%%%%%%%%% Column 6 : Wx=2*sigma of Gaussian %%%%%%%%%%%%%%
    %%%%%%%%%%%% Column 7 : Wy=2*sigma of Gaussian %%%%%%%%%%%%%%
XYZt(:,6) = dataWidthX;                      % X width = 2*sigma of Gassian peak
XYZt(:,7) = dataWidthY;                      % Y width = 2*sigma of Gassian peak

     %%%%%%%%%%%% Column 8 : Background of Gaussian %%%%%%%%%%%%%
XYZt(:,8) = dataMList.bg;                    % Background of Gassian fit

     %%%%%%%%%%%% Column 9 : Nothing N/A %%%%%%%%%%%%%%%%%%%%%%%%
dataLength= dataMList.length;
dataLength = double(dataLength);
XYZt(:,9) = dataLength; 
disp('---------------------------------------------------------------');
judge_length=input('Do you want to delete dots repeatly localized? Press Enter - ','s');
if isempty(judge_length)
    frame_length=input('Type in a frame length for delete (usually 3 frames) - ');
    index_length=find(dataLength>frame_length);
    LengthStr=strcat('There are -',num2str(length(index_length)),'- out of -',num2str(length(dataLength)),'- dots have been repeatly localized for more than -',num2str(frame_length),'- times');
    disp('*********************************************************************************************')
    disp(LengthStr);
    disp('*********************************************************************************************')
end

     %%%%%%%%%%%% Column 10: Nothing N/A %%%%%%%%%%%%%%%%%%%%%%%%
XYZt(:,10)= 0; 

     %%%%%%%%%%%% Column 11: Ellipticity Ratio %%%%%%%%%%%%%%%%%%
XYZt(:,11)= dataAx;

     %%%%%%%%%%%% Column 12: Frame Number %%%%%%%%%%%%%%%%%%%%%%%
XYZt(:,12)= dataMList.frame; 
dataFrame = XYZt(:,12);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

XYZt(index_length,:)=[];
dataAx(index_length,:)=[];
dataEllip(index_length,:)=[];
dataFrame(index_length,:)=[];
dataHeight(index_length,:)=[];
dataLength(index_length,:)=[];
dataWidth(index_length,:)=[];
dataWidthX(index_length,:)=[];
dataWidthY(index_length,:)=[];
%%

if min(dataEllip)==0 && max(dataEllip)==0
    Dimension='2';
    disp('-------------------------------------------------------------')
    disp('This is 2D STORM imaging, right? If not, please stop and check it...');
    judge_color='z';
    color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
    while ~color_continue
        judge_color=input('Type -r- for Red Channel, -g- for Green, -b- for Blue - ','s');
        color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
        if isempty(judge_color)
            color_continue=0;
        end
    end
    XYZt(:,3)=0;    % for color infomation only
else
    Dimension='3';
    disp('-------------------------------------------------------------')
    disp('This is 3D STORM imaging, right? If not, please stop and check it...');
    judge_input23=input('Type in 2 for input 2D data only, 3 for input 3D data - ');
    if judge_input23==2
        Dimension='2';
        judge_color='z';
        color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
        while ~color_continue
            judge_color=input('Type -r- for Red Channel, -g- for Green, -b- for Blue - ','s');
            color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
            if isempty(judge_color)
                color_continue=0;
            end
        end
        XYZt(:,3)=0;    % for color infomation only
    else    % input 3D info from insight or calculate in matlab
        judge_color='z';
        color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
        while ~color_continue
            judge_color=input('Type -r- for Red Channel, -g- for Green, -b- for Blue - ','s');
            color_continue=judge_color=='r' | judge_color=='g' | judge_color=='b';
            if isempty(judge_color)
                color_continue=0;
            end
        end
        disp('Wanna Input Z coordinate from Insight3 or Calculate by Matlab???');
        judge_Insight_Z=input('Type in -i- for Insight3, or -m- for Matlab - ','s');
        if judge_Insight_Z=='i'
            if judge_file=='c'
                XYZt(:,3)=dataMList.zc;
            else
                XYZt(:,3)=dataMList.z;
            end
            Zrange=[min(XYZt(:,3)) max(XYZt(:,3))];
            % modified on Nov 10 2016 by Q.Peter SU
        elseif judge_Insight_Z=='m'
            if judge_color=='r'
                dataZ=p_R(1) * dataEllip.^3 + ...
                      p_R(2) * dataEllip.^2 + ...
                      p_R(3) * dataEllip.^1 + ...
                      p_R(4) ;
                XYZt(:,3)=dataZ;
            elseif judge_color=='g'
                dataZ=p_G(1) * dataEllip.^3 + ...
                      p_G(2) * dataEllip.^2 + ...
                      p_G(3) * dataEllip.^1 + ...
                      p_G(4) ;
                XYZt(:,3)=dataZ;
            else
                disp('*******************************');
                disp('Wrong Color Channle Information');
                disp('*******************************');
            end

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%% Delete Abnormal Dots out of Z Range (Zmin ~ Zmax) %%%%%%%%%%%%
            disp('-------------------------------------------------------------')
            ZrangeStr=strcat('Zmin=',num2str(min(XYZt(:,3))),'-','Zmax=',num2str(max(XYZt(:,3))));
            disp(ZrangeStr);
            disp('for the former data in BIOPIC No3 STORM, range should be -300~250')
            Zrange=input('Type in Z Range [Zmin Zax], or press Enter for all - ');
            if isempty(Zrange)==1
                Zrange=[min(XYZt(:,3)) max(XYZt(:,3))];
            else
            end
            index1 = find(XYZt(:,3)<= Zrange(1));
            index2 = find(XYZt(:,3)>= Zrange(2));
            index = [index1 ; index2];
            XYZt(index,:)=[];
            XYZt(:,3) = XYZt(:,3) * 0.72;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Zrange=[-700 700];
index_Z1 = find(XYZt(:,3)<= Zrange(1));
index_Z2 = find(XYZt(:,3)>= Zrange(2));
index_Z = [index_Z1 ; index_Z2];
XYZt(index_Z,:)=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

moleculeNo=size(XYZt,1);
strout=['Totally, there are -' num2str(moleculeNo) '- out of -'...
    num2str(size(dataMList.x,1)) '- dots/localizations'];
disp('-------------------------------------------------------------');
disp(strout);
disp('-------------------------------------------------------------');


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
judge_BF=input('Delete Dots in BF? Press Enter for Sure or Any Other Key to Skip - ','s');
if isempty(judge_BF)
    [FileNameIndex, PathNameIndex] = uigetfile('*BfIndex.txt',...
    'Index File Output by ST22_IndexBF_and_GenDrift','MultiSelect','off');
    AfnameIndex=char(FileNameIndex);
    disp('-------------------------------------------------------------');
    cd(PathNameIndex);disp(PathNameIndex);disp(FileNameIndex);
    BFIndex = importdata(AfnameIndex);
    BFIndexDiff=[1 diff(BFIndex)];
    BFIndexDiffDiff=[0 diff(BFIndexDiff)];
    BFIndexCycle=find(BFIndexDiffDiff>0);
    BFIndexNo=length(BFIndexCycle);
    rangeDel=2
    BFIndexDel=BFIndex;
    
    tic
    % for the first BF cycle only delete later frames
    indextemp=[];
    for i=1:rangeDel
        indextemp=[indextemp BFIndex(BFIndexCycle(1)-1)+i];
    end
    BFIndexDel=[BFIndexDel indextemp];
    % for the last BF cycle only delete former frames
    indextemp=[];
    for i=1:rangeDel
        indextemp=[indextemp BFIndex(BFIndexCycle(BFIndexNo))-i];
    end
    BFIndexDel=[BFIndexDel indextemp];
    % for the rest frames
    for j=2:(BFIndexNo-1)
        indextemp=[];
        for i=1:rangeDel
            indextemp=[indextemp BFIndex(BFIndexCycle(j-1))-i];
            indextemp=[indextemp BFIndex(BFIndexCycle(j)-1)+i];
        end
        BFIndexDel=[BFIndexDel indextemp];
    end
    
    BFIndexDel=sort(BFIndexDel);
    
    BFIndexDelIndex=[];
    for k=1:length(BFIndexDel)
        indextemp=find(XYZt(:,12)==BFIndexDel(k));
        BFIndexDelIndex=[BFIndexDelIndex;indextemp];
        if mod(k,2000)==0
            strout=['Frame No. = ',num2str(k),'/',num2str(length(BFIndexDel)),'-Frames'];
            disp(strout);
         end
    end

%     tic
%     rangeDel=2
%     for i=1:length(BFIndex)
%         range=BFIndex(i)-rangeDel : BFIndex(i)+rangeDel;        
%         for j=1:(rangeDel*2+1)
%             indextemp=find(XYZt(:,12)==range(j));
%             BFindexDel=[BFindexDel;indextemp];
%         end
%         if mod(i,2000)==0
%             strout=['Frame No. = ',num2str(i),'/',num2str(length(BFIndex)),'-Frames'];
%             disp(strout);
%         end
%     end
%     toc

    
    XYZt(BFIndexDelIndex,:)=[];
    dataLength(BFIndexDelIndex,:)=[];
    deletepoint=size(dataHeight,1)-size(XYZt,1);
    toc
    
    str=['Totally, there are -' num2str(deletepoint) '- dots deleted from the BF frames.'];
    disp(str);
    
%     dataAx(BFindexDel,:)=[];   
%     dataEllip(BFindexDel,:)=[];   
%     dataFrame(BFindexDel,:)=[];   
%     dataHeight(BFindexDel,:)=[];   
%     dataLength(BFindexDel,:)=[];   
%     dataWidth(BFindexDel,:)=[];   
%     dataWidthX(BFindexDel,:)=[];   
%     dataWidthY(BFindexDel,:)=[];   
%     dataZ(BFindexDel,:)=[];   
   
end

% save plot time
%%
 tic
 
if size(XYZt,1)>2000000
    disp('-------------------------------------------------------------')
    moleculeStr=strcat('There are -',num2str(size(XYZt,1)),'- localizations left');
    disp(moleculeStr);
    MolLim=input('Please Type in a Smaller Number for Display (Enter for All) - ');
   
     if isempty(MolLim)==1
         MolLim=size(XYZt,1);    
     else
     end
else
    MolLim=size(XYZt,1);
end

close all
scrsz = get(0,'ScreenSize');
figure('OuterPosition',[1 0.5*scrsz(4) 0.75*scrsz(4) 0.5*scrsz(4)]);
                        %[left bottom width height]
whitebg([0 0 0]);hdl1=gcf;clf;
figure('OuterPosition',[1 1  0.75*scrsz(4) 0.5*scrsz(4)]);
whitebg([0 0 0]);hdl2=gcf;clf;
figure(hdl1)
if Dimension=='2'
    plot3k([XYZt(1:MolLim,1)*1000,XYZt(1:MolLim,2)*1000,XYZt(1:MolLim,3)],...
            XYZt(1:MolLim,3),[-100 100],{'o',1});
        view([0 -90]);
elseif Dimension=='3'
    plot3k([XYZt(1:MolLim,1)*1000,XYZt(1:MolLim,2)*1000,XYZt(1:MolLim,3)],...
        XYZt(1:MolLim,3),Zrange,{'o',1});
end
%scatter3(XYZt(1:MolLim,1),XYZt(1:MolLim,2),XYZt(1:MolLim,3)/1000,1,XYZt(1:MolLim,3)/1000,'filled');
xlabel('X Coordinate (nm)');ylabel('Y Coordinate (nm)');
zlabel('Z Coordinate (nm)');
axis equal


if judge_file=='u'
    if isempty(judge_BF)
        title({['Final ' Dimension 'D HR Image (Without Auto-Insight3-Correction), Delete BF.'];...
            [num2str(MolLim),'-out of-',num2str(moleculeNo),'-Localizations after Deleting BF Dots.'];...
            ['Colormap (Blue to Red) Means Z Position (Low to High) in nm.']});
    else
        title({['Final ' Dimension 'D HR Image (Without Auto-Insight3-Correction), Without Del BF.'];...
            [num2str(MolLim),'-out of-',num2str(moleculeNo),'-Localizations.'];...
            ['Colormap (Blue to Red) Means Z Position (Low to High) in nm.']});
    end
elseif judge_file=='c'
    if isempty(judge_BF)
        title({['Final ' Dimension 'D HR Image (With Insight3 Auto-Insight3-Correction), Delete BF.'];...
            [num2str(MolLim),'-out of-',num2str(moleculeNo),'-Localizations after Deleting BF Dots.'];...
            ['Colormap (Blue to Red) Means Z Position (Low to High) in nm.']});
    else
        title({['Final ' Dimension 'D HR Image (With Insight3 Auto-Insight3-Correction), Without Del BF.'];...
            [num2str(MolLim),'-out of-',num2str(moleculeNo),'-Localizations.'];...
            ['Colormap (Blue to Red) Means Z Position (Low to High) in nm.']});
    end
end
%set(gca,'YDir','reverse');colorbar;


% 06th July 2017 by Q.Peter SU
if isempty(judge_length)
    outputfilename = strcat(Afname(1:end-4),'-',num2str(frame_length),'Rep');
else
    outputfilename = Afname(1:end-4);
end

if judge_file=='c'    
    if isempty(judge_BF)
        figurefilename1_1 = strcat(outputfilename,'-',Dimension,'D-AC-DelBF1', '.fig');
        figurefilename1_2 = strcat(outputfilename,'-',Dimension,'D-AC-DelBF1', '.bmp');
    else
        figurefilename1_1 = strcat(outputfilename,'-',Dimension,'D-AC1', '.fig');
        figurefilename1_2 = strcat(outputfilename,'-',Dimension,'D-AC1', '.bmp');
    end
elseif judge_file=='u'
    if isempty(judge_BF)
        figurefilename1_1 = strcat(outputfilename,'-',Dimension,'D-UC-DelBF1', '.fig');
        figurefilename1_2 = strcat(outputfilename,'-',Dimension,'D-UC-DelBF1', '.bmp');
    else
        figurefilename1_1 = strcat(outputfilename,'-',Dimension,'D-UC1', '.fig');
        figurefilename1_2 = strcat(outputfilename,'-',Dimension,'D-UC1', '.bmp');
    end
end
disp(figurefilename1_1);
saveas(gcf,figurefilename1_1,'fig');
disp(figurefilename1_2);
saveas(gcf,figurefilename1_2,'bmp');




figure(hdl2)
whitebg([0 0 0]);

%%%%%%%%%%%%%%%%%%%%%%%%% sub fig 1 : width to peak %%%%%%%%%%%%%%%%%%%%%%
subplot(2,3,1);
%scatter(dataWidth(1:MolLim),dataHeight(1:MolLim),1,'filled');
dataWidthPlot=XYZt((1:MolLim),6).*XYZt((1:MolLim),7);
dataWidthPlot=sqrt(dataWidthPlot);
plot(dataWidthPlot(1:MolLim),XYZt((1:MolLim),4),'MarkerSize',1,'Marker','.',...
    'LineStyle','none','Color',[1 1 0]);
hold on
plot([300 250 650 300],[16000 1500 1500 16000],'Marker','o','LineWidth',2,...
    'LineStyle','--','Color',[1 0 0]);
hold on
plot([400 200 850 400],[10000 1000 1000 10000],'Marker','o','LineWidth',2,...
    'LineStyle','--','Color',[0 1 0]);
hold off
xlim([100 1000]);ylim([0 20000]);
xlabel('Width of Gaussian Fit');ylabel('Height of the Gaussian Fit');
title('Width vs Intensity (red2D-green3D)');grid on;zoom on;

%%%%%%%%%%%%%%%%%%%%%%%%% sub fig 2 : hist of height %%%%%%%%%%%%%%%%%%%%%%
subplot(2,3,2);
hist(dataHeight(1:MolLim),100);xlim([0 12000]);xlabel('Intensity (A.U.)');
ylabel('Occurrence');view([90 -90]);title('Hist of Intenstiy of Every Dot'); 
grid on;zoom on;

%%%%%%%%%%%%%%%%%%%%%%%%% sub fig 3 : hist of ellipticity %%%%%%%%%%%%%%%%%%
subplot(2,3,3)
hist(XYZt(1:MolLim,3),100);xlim([-750 750]);xlabel('Z Position (nm)');
ylabel('Occurrence');title('Hist of Z Position');grid on;zoom on;

%%%%%%%%%%%%%%%%%%%%%%%%% sub fig 4 : hist of width %%%%%%%%%%%%%%%%%%
subplot(2,3,4);
hist(dataWidthPlot(1:MolLim),100);xlim([100 1000]);
xlabel('Width of Gaussian Fit(nm)');
ylabel('Occurrence');title('Hist of Width (=SQRT(Wx*Wy))');grid on;zoom on;
set(gca,'YDir','reverse');

%%%%%%%%%%%%%%%%%%%%%%%%% sub fig 5 : hist of dots in frame %%%%%%%%%%%%%%%
subplot(2,3,5);
hist(XYZt((1:MolLim),12),max(XYZt((1:MolLim),12)));xlabel('Frame No.(A.U.)');
ylabel('Dots Counts(A.U.)');title('Dots Counts in Every Frame');grid on;zoom on;

%%%%%%%%%%%%%%%%%%%%%%%%% sub fig 6 : hist of dots length %%%%%%%%%%%%%%%%
subplot(2,3,6);
hist(dataLength(1:MolLim),max(dataLength));xlabel('Dots Length(A.U.)');xlim([0 30]);
ylabel('Occurrence');title('Hist of Length of Every Dot');grid on;zoom on;

% outputfilename = Afname(1:end-4);
if judge_file=='c'    
    if isempty(judge_BF)
        figurefilename2_1 = strcat(outputfilename,'-',Dimension,'D-AC-DelBF2', '.fig');    
        figurefilename2_2 = strcat(outputfilename,'-',Dimension,'D-AC-DelBF2', '.bmp');    
    else
        figurefilename2_1 = strcat(outputfilename,'-',Dimension,'D-AC2', '.fig'); 
        figurefilename2_2 = strcat(outputfilename,'-',Dimension,'D-AC2', '.bmp'); 
    end
elseif judge_file=='u'    
    if isempty(judge_BF)
        figurefilename2_1 = strcat(outputfilename,'-',Dimension,'D-UC-DelBF2', '.fig');    
        figurefilename2_2 = strcat(outputfilename,'-',Dimension,'D-UC-DelBF2', '.bmp');    
    else
        figurefilename2_1 = strcat(outputfilename,'-',Dimension,'D-UC2', '.fig'); 
        figurefilename2_2 = strcat(outputfilename,'-',Dimension,'D-UC2', '.bmp'); 
    end
end
disp(figurefilename2_1);
saveas(gcf,figurefilename2_1,'fig');
disp(figurefilename2_2);
saveas(gcf,figurefilename2_2,'bmp');


%Output for DriftCorrect and precisionstat.m
% outputfilename = Afname(1:end-4);
if judge_file=='c'    
    if isempty(judge_BF)
        outDATfilename  = strcat(outputfilename,'-',Dimension,'D-AC-DelBF', '.dat');    
    else
        outDATfilename  = strcat(outputfilename,'-',Dimension,'D-AC', '.dat');    
    end
elseif judge_file=='u'
    if isempty(judge_BF)
        outDATfilename  = strcat(outputfilename,'-',Dimension,'D-UC-DelBF', '.dat');    
    else
        outDATfilename  = strcat(outputfilename,'-',Dimension,'D-UC', '.dat');    
    end
end

save(outDATfilename, 'XYZt', '-ascii', '-tabs');
disp(outDATfilename);
toc
disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')

%% Part III - Drift Correction for Single-color if necessary
if judge_file=='c'
    % nothing
elseif judge_file=='u'   
    disp('-------------------------------------------------------------');
    judge_Drift=input('Run DriftCorr? Press Enter for Sure or Any Other Key to Skip - ','s');
    if isempty(judge_Drift)
        ST05_DriftCorrect_3D(outDATfilename);
        outDATfilename=strcat(outDATfilename(1:end-4),'-C2D.dat');
        disp(outDATfilename);
    end
end
disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')

%% Part III - Drift Correction for Dual-color if necessary 
if judge_color=='g'
    disp('-------------------------------------------------------------');
    disp('You are Analyzing Green Channel, if This is a Dual-Color Imaging,');
    disp('Type in -d- for Dual-Color, Then Do the Two-Color DriftCorr - ');
    judge_color_dual=input('Or Just Press Enter Key to Skip This Step - ','s');
    if judge_color_dual=='d'   
        ST05_DriftCorrect_3D(outDATfilename);
        outDATfilename=strcat(outDATfilename(1:end-4),'-C2D.dat');
        disp(outDATfilename);
    end
end

if judge_color=='b'
    disp('-------------------------------------------------------------');
    disp('You are Analyzing Blue Channel, if This is a Dual-Color Imaging,');
    disp('Type in -d- for Dual-Color, Then Do the Two-Color DriftCorr - ');
    judge_color_dual=input('Or Just Press Enter Key to Skip This Step - ','s');
    if judge_color_dual=='d'   
        ST05_DriftCorrect_3D(outDATfilename);
        outDATfilename=strcat(outDATfilename(1:end-4),'-C2D.dat');
        disp(outDATfilename);
    end
end

disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')

%% Part IV - Delete BF frames
% disp('-------------------------------------------------------------');
% judge_BF=input('Delete Dots in BF? Press Enter for Sure or Any Other Key to Skip - ','s');
% if isempty(judge_BF)    
%     ST11_2_DeleteBFframe(outDATfilename);
%     outDATfilename=strcat(outDATfilename(1:end-4),'-DelBF.dat');
% end

% %% Part IV - Delete Points  Modified by Q Peter Su on Feb 28th 2017
% disp('-------------------------------------------------------------');
% judge_Point=input('Delete UnExpected Dots? Press Enter for Sure or Any Other Key to Skip - ','s');
% if isempty(judge_Point)    
%     ST11_1_DeletePoint(outDATfilename)
%     outDATfilename=strcat(outDATfilename(1:end-4),'-DelP1.dat');
%     disp(outDATfilename);
% end
% disp('------------------------------------------------------------------')
% disp('Date and Time Now ')
% disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
% disp('------------------------------------------------------------------')

%% Part IV - Delete Points  Modified by Q Peter Su on Feb 28th 2017
disp('-------------------------------------------------------------');
judge_Point=input('Delete UnExpected Dots? Press Enter for Sure or Any Other Key to Skip - ','s');
if isempty(judge_Point)    
    ST11_1_DeletePoint(outDATfilename)
    outDATfilename=strcat(outDATfilename(1:end-4),'-DelP1.dat');
    disp(outDATfilename);
end
disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')

%% Part V - precisionstat
disp('-------------------------------------------------------------');
judge_Precis=input('Run precisionstat? Press Enter for Sure or Any Other Key to Skip - ','s');
if isempty(judge_Precis)
    minclustersize=[];
    while isempty(minclustersize)
        minclustersize=input('Type in the MinClusterSize for Precision Calculation. - ');
    end
    ST07_precisionstat_3D(outDATfilename,minclustersize);
    outDATfilename=strcat(outDATfilename(1:end-4),'-Thrd',...
                          num2str(minclustersize),'-Filt-DAT.dat');
    disp(outDATfilename);
end
disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')


%% Part IV - Delete Points  Modified by Q Peter Su on Feb 28th 2017
disp('-------------------------------------------------------------');
judge_Point=input('Delete UnExpected Dots? Press Enter for Sure or Any Other Key to Skip - ','s');
if isempty(judge_Point)    
    ST11_1_DeletePoint(outDATfilename)
    outDATfilename=strcat(outDATfilename(1:end-4),'-DelP1.dat');
    disp(outDATfilename);
end
disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')

%% Part VI - fastmap
disp('-------------------------------------------------------------');
judge_fastmap=input('Run fastmap.m? Press Enter for Sure or Any Other Key to Skip - ','s');
if isempty(judge_fastmap)
    binsize=5;
    strout=['Binsize of the Rendered Image = ' num2str(binsize)];
    disp(strout)
    resolution=25;
    strout=['Resolution of the STORM Image = ' num2str(resolution)];
    disp(strout)
    ST08_1_fastmap(outDATfilename,Xpixel,Ypixel,binsize,resolution);
end
disp('------------------------------------------------------------------')
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')

%% the end of the script
disp('-------------------------------------------------------------')
disp('Now, Please enjoy your STORM image, bye...')
disp(outDATfilename);
disp('Date and Time Now ')
disp(datestr(now,'mmm_dd_yyyy-HH_MM_SS'));
disp('------------------------------------------------------------------')
diary off


end

