%% Yujie Sun on Sep 11 2012
% Used to quantify the localization precision of STORM data by computing the
% spreading of each manually picked clusters. Yujie Sun 2012-09-11

% clusterlocation is a matrix that contains the coordinates of clusters from
% an analyzed image. When it is called, this function is used to compare
% precision between a pair of drift-corrected and uncorrected images. 
%
%% Revised by Q.Peter Su since Sep 2012 @ Harvard
%  petersq001@gmail.com
%
%% Revised by Q.Peter Su since Dec 2012 @ BIOPIC, PKU
%  petersq001@gmail.com
%
%% Input:
%   1. 
%   1. raw ST*.dat file of localizations/dots output by Merge_and_plotfile.m
%            *.dat file of localizations/dots output by Insight_to_Matlab.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.
%
%   2. 'Already had cluster location and want to do a comparison?'
%        default n (no)
%   3. 'Finding clusters: Manual pick or Automation?'
%        default a (auto)
%   4. Threshold of Min Size for Cluster 
%        Switching cycle of the dye molecule.
%
%% Output:
%   1. -Thrd#-Filtered-Auto-Quant.fig
%       fig of the unfiltered localizations/dots in yellow on the left
%       fig of the filtered   localizations/dots in green  on the right
%       fig of the filtered   clusters           in red    on the right
%           title({['stdx = ' num2str(xres,3) '   stdy = ' num2str(yres,3)]; ...
%           ['stdxy = ' num2str(xyres,3) '   stdz = ' num2str(zres,3)]; ...
%           ['cluster size threshold = ' num2str(minclustersize,2)];...
%           ['total dots = ' num2str(length(dots2plot), '%10.1d\n')];...
%           ['dots per cluster = ' num2str(dotspercluster,2)];...
%           ['cluster number = ' num2str(length(clusters), '%10.1d\n')];...
%           ['mean Intensity = ' num2str(meanInt, '%10.1d\n')];},...
%            'FontWeight','bold','FontSize',14,'FontName','Arial',...
%            'Color',[1 1 1]);
%
%   2. -Thrd#-Filtered-precisionstat-Auto.mat
%       'precision' Matrix of the information of all the filtered clusters
%                   total clusters no. in row and 10 columns
%       Column 1 : mean value of x axis - in um 
%       Column 2 : std value of x axis  - in nm
%       Column 3 : mean value of y axis - in um
%       Column 4 : std value of y axis  - in nm
%       Column 5 : mean value of z axis - in um
%       Column 6 : std value of z axis  - in nm
%       Column 7 : std value of x&y     - in nm
%       Column 8 : number of total dots 
%       Column 9 : mean intensity       
%       Column 10: cluster no.     - same as the clusterpicked.mat in row
%
%   3. -Thrd#-Filtered-clusterspicked-Auto.mat
%       'clusters' cell - total clusters no. in row and 1 Column
%        each cell is a cluster containing all the localizations/dots in it
%        in each cell
%        Column 1 : X in nano-meter
%        Column 2 : Y in nano-meter
%        Column 3 : Z in nano-meter
%
%   4. -Thrd#-Filtered-DAT.dat
%        Matrix of all the filtered localizations/dots
%        information is same as the output of Merge_and_plotfiles.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.
%
%% Revised by Q. Peter Su on Nov 16 2012 @ Harvard
%  abs(mol(:,3)-mol(1,3))<dz); 
%
%% Revised by Q. Peter Su on Jan 02 2013  @ BIOPIC, PKU
% Column 10: cluster no.
% for Domain3D.m identify dots in clusters in domains
%
%% Revised by Q.Peter Su on Jan 09 2013 @ BIOPIC, PKU
%  Input global variable 'minclustersize'
%
%% Revised by Yujie Sun on Jan 30 2013 @ BIOPIC, PKU
%  jj=1;
%  while abs(mol(jj,1)-mol(1,1))<3*dx && jj<size(mol,1)
%      jj=jj+1;
%  end
%  indx1 = find(abs(mol(1:jj,1)-mol(1,1))<dx &...
%               abs(mol(1:jj,2)-mol(1,2))<dy &...
%               abs(mol(1:jj,3)-mol(1,3))<dz);
%
%
%% Revised by Q.Peter Su on May 20 2013 @ BIOPIC, PKU
%  ST07_precisionstat_3D(Afname,minclustersize)
%
%% Marked by Q.Peter Su on Nov 21 2013 @ BIOPIC, PKU
%  XYZt_out=ST07_precisionstat_3D(FileName_or_XYZt_in, minclustersize)
%
%% Revised by Q.Peter Su on Dec 13 2013 @ BIOPIC, PKU
%  Fast 'precisionstat.m' by set a range of search for 'find' function
%  and, remember, DONOT use 'a(x,:)=[];', because 'a(x,:)=0;' is much much faster
%
%         pp=find(mol(:,1)>0);
%         while size(pp,1)>minclustersize
%              k=1;j=pp(k);jj=pp(k);
%              while abs(mol(jj,1)-mol(j,1))<2*dx
%                  k=k+1;jj=pp(k);% find all the nonzero dots within 2*dx
%              end                   
%              indx1 = find(abs(mol(j:jj,1)-mol(j,1))<2*dx &...
%                           abs(mol(j:jj,2)-mol(j,2))<2*dy &...
%                           abs(mol(j:jj,3)-mol(j,3))<2*dz);
%              indx1=indx1+j-1;
%
%% Revised by Q.Peter Su on Dec 16 2013 @ BIOPIC, PKU
%            while abs(mol(jj,1)-mol(j,1))<2*dx && pp(k+1)<totaldots
%                k=k+1;jj=pp(k); % find all the nonzero dots within 2*dx
%            end                   
%
%% Revised by Q.Peter SU on 19 Sep. 2017 @ HRI, USYD
%  Try to fix the problem of "Blade shape" like shape at the edge of clusters,
%     especially when the total molecule density is high.
%  IF a deleted dot is closed to (dis<dx=60nm) a dot already in a cluster,
%     combine them in same cluster.
%
%         Double_check_dis=(mol_del(i,1)-(mol_temp(clusterindx,1))).^2+...
%                          (mol_del(i,2)-(mol_temp(clusterindx,2))).^2+...
%                          (mol_del(i,3)-(mol_temp(clusterindx,3))).^2/1000000;
%         Double_check_indx= find(Double_check_dis==0);
%         Double_check_dis(Double_check_indx)=[];
%         clusterindx(Double_check_indx)=[];
%         Double_check_indx= Double_check_dis==min(Double_check_dis);
%         Double_check_cluster_no=mol_temp(clusterindx(Double_check_indx),13);
%         if Double_check_cluster_no~=0 && min(Double_check_dis)<dx
%             dots2plot_temp=[dots2plot_temp; mol_del(i,1:12)];
%             cluster_temp=[mol_del(i,1)*1000, mol_del(i,2)*1000, mol_del(i,3), mol_del(i,12), mol_del(i,4)];
%                               % X in nm          Y in nm            Z in nm      frame no.       intensity
%             clusters{Double_check_cluster_no,1}=[clusters{Double_check_cluster_no,1};cluster_temp];
%             Double_check_dot_no=Double_check_dot_no+1;
%
%% Revised by Q.Peter SU on 19 Sep. 2017 @ HRI, USYD
%  shorten the file name (sometimes followed by Domain3D)
%     saveas(gcf, [Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
%         '-Filt-Quant.fig'], 'fig');
%     saveas(gcf, [Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
%         '-Filt-Quant.bmp'], 'bmp');
%     save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
%         '-Filt-Prec.mat'],'precision');
%     save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
%         '-Filt-cluster.mat'],'clusters');
%     save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
%         '-Filt-DAT.dat'], 'dots2plot', '-ascii', '-tabs');
%
%% Revised by Q.Peter SU on 28 Feb 2020 at CPC, USYD
%         Double_check_indx= Double_check_dis==min(Double_check_dis); 
%         % find the dot nearest to the No.i dot, but there might be multiple
%         % dots with the same nearest distance.
%         Double_check_cluster_no=mol_temp(clusterindx(Double_check_indx),13);
%         Double_check_cluster_no=Double_check_cluster_no(1);
%         % sometimes, code will find multiple dots with same nearest distance.
%

function [xres, yres, zres, xyres, dotspercluster]=ST07_precisionstat_3D(Afname,minclustersize)
%clc
close all
%clear all;
%% impport data
if nargin<2
[FileName, Pathname]=uigetfile('*.dat',...
    'pick the STORM data to quantify localization precision',...
    'MultiSelect','off');
cd(Pathname);
Afname=char(FileName);

disp(Afname)
disp('Input Threshold of Min Size for Cluster:')
minclustersize=input('-2- times for Cy3; -5- times for Atto550 or Cy5 - ');
         %minimal number of dots per clusters
end

Data=[];
disp('Import Data ... ing ...');
tic
Data=importdata(Afname);
toc
disp('------------------------------------------------------------------')
strout=strcat('Totally, There are -', num2str(size(Data,1)), '- Localizations.');
disp(strout)
%% Calculate cluster
tic
disp('Calculate Clusters ...ing...');
[YY,I]=sort(Data(:,1)); %sort dots by x positions
Data=Data(I,:);
x=Data(:,1);       %X in um
y=Data(:,2);       %Y in um
z=Data(:,3)/1000;  %Z in nm to um
Intensity=Data(:,4);
meanIall=mean(Intensity);
totaldots=size(Data,1);

comparison='n';
%input('Already had cluster location and want to do a comparison? press y or n:  ','s');

if comparison=='y'
    [Fname, Pname]=uigetfile('*.mat','pick the cluster location file','MultiSelect','off');
    pos_cluster=load(Fname);
    pos=pos_cluster.pos;
    [precision, clusters, dots2plot]=clusterstat(pos,x,y,z,Intensity);
else
    pos=[];
    decision='a';%input('Finding clusters: Manual pick or Automation? press m or a:  ','s');
    if decision=='m'  %munually pick clusters for quantification
       fig1=figure;
       % Create axes
       axes1 = axes('Visible','on','Parent',fig1,'YDir','reverse','View',[0 90]);
       % Create plot
       p1=plot(x,y,'.r');
       axis equal;
       grid off;
       display('Zoom in a region for selection');
       pause;
       buttonpress=0;
       display('Draw rectangles to pick multiple clusters for analysis');
       display('Adjust box size and position, double-click in the box to confirm,');
       display(' and click anywhere in the figure to draw next box');
       display('After confirm the last box, press any letter button to finish');
       while buttonpress~=1
             h=imrect;
             position=wait(h); %Interactively place a rectangle by clicking and dragging. 
                               %Use wait to block the MATLAB command lines. 
                               %Double-click on the rectangle to resume execution of the MATLAB command line.
             if isvalid(h)
                 pos = [pos; getPosition(h)];
             else
                 return
             end
             buttonpress=waitforbuttonpress;
       end
       save([Afname '_clusterlocation.mat'],'pos');
       saveas(gcf, [Afname(1:end-4) '_manualpick'], 'fig');
       [precision, clusters, dots2plot]=clusterstat(pos,x,y,z,Intensity);

       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Finding clusters: Automation = a (default) %%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    elseif decision=='a' %automation
        %definition of a good cluster 
        dx = 0.06; %um;
        dy = 0.06; %um;
        dz = 0.12*1000; %nm;
        
        mol=Data;
        mol_del=[]; % revised by Q.Peter SU on 19 Sep 2017 for 'Blade shape' problem
        precision=[];     %Yujie Sun
        dots2plot=[];     %algorithm is modified from PulldatafromFiesta_5.m
        i=1;
        totaldotsStr=num2str(totaldots);
        pp=find(mol(:,1)>0);
        while size(pp,1)>minclustersize
             k=1;j=pp(k);jj=pp(k);
             while abs(mol(jj,1)-mol(j,1))<2*dx &&...
                   abs(mol(jj,1)-mol(j,1))<2*dy &&...
                   pp(k+1)<totaldots
                 k=k+1; % K is counting
                 jj=pp(k); % find all the nonzero dots within 2*dx and 2*dy
             end
             indx1 = find(abs(mol(j:jj,1)-mol(j,1))<2*dx &...
                          abs(mol(j:jj,2)-mol(j,2))<2*dy &...
                          abs(mol(j:jj,3)-mol(j,3))<2*dz);
             indx1=indx1+j-1;
       %%%%%%%%%%%%%%%%%%%%%% Cluster Picking %%%%%%%%%%%%%%%%%%%%%%%%%%%%                       
             if size(indx1,1)>minclustersize
                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 %%%% if Size(indx1,1)>minclustersize, do the folowing %%%%
                 %%%% use the median value dots as the center, find dots %%
                 %%%% around it within 60nm in X and Y, and 120nm in Z %%%%
                 %%%%%%%%%%%%%% This is a ture cluster of dots. %%%%%%%%%%%
                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                medscatter=(mol(indx1,1)-median(mol(indx1,1))).^2+...
                           (mol(indx1,2)-median(mol(indx1,2))).^2+...
                           (mol(indx1,3)-median(mol(indx1,3))).^2/1000000;
                medindx1=find(medscatter==min(medscatter));
                medindx=indx1(medindx1(1));
                %medianxy=mol(medindx,:);
                clusterindx = find(abs(mol(indx1,1)-mol(medindx,1))<dx &...
                                   abs(mol(indx1,2)-mol(medindx,2))<dy &...
                                   abs(mol(indx1,3)-mol(medindx,3))<dz);
                if length(clusterindx)>minclustersize
                   indx4=indx1(clusterindx);
                   indx1=indx4;   
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   numdots1=length(indx1);
                   meanx1=mean(mol(indx1,1));       %um
                   meany1=mean(mol(indx1,2));       %um
                   meanz1=mean(mol(indx1,3))/1000;  %um
                   stdx1=std(mol(indx1,1))*1000;    %nm
                   stdy1=std(mol(indx1,2))*1000;    %nm
                   stdz1=std(mol(indx1,3));         %nm
                   meanI=mean(mol(indx1,4));
                   precision=[precision; [meanx1 stdx1 meany1 stdy1 meanz1...
                              stdz1 sqrt(stdx1^2+stdy1^2) numdots1 meanI i]];
                     %Matrix 'Precision'
                     % Column 1 : mean value of x axis in um
                     % Column 2 : std value of x axis  -this is x resolution for every cluster in nm
                     % Column 3 : mean value of y axis in um
                     % Column 4 : std value of y axis  -this is y resolution for every cluster in nm
                     % Column 5 : mean value of z axis  in um
                     % Column 6 : std value of z axis  -this is z resolution for every cluster in nm
                     % Column 7 : std value of x&y    -this is xy resolution for every cluster in nm
                     % Column 8 : number of total dots
                     % Column 9 : mean intensity
                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                     % Column 10: cluster no.
                     % Revised by Q. Peter Su on Jan 02 2013 for Domain3D.m identify dots in clusters in domains
                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                           
                  clusters{i,1}=[mol(indx1,1)*1000, mol(indx1,2)*1000, mol(indx1,3), mol(indx1,12), mol(indx1,4)];
                                 % X in nm         Y in nm          Z in nm         frame no.        intensity
                                 % revised by Q.Peter Su on Dec 19 2013
                                 % revised by Q.Peter Su on Mar 02 2014
                  dots2plot=[dots2plot; mol(indx1,:)];
                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  Data(indx1,13)=i;
                  % add molecule cluster index
                  % revised by Q.Peter SU on 19 Sep 2017 for 'Blade shape' problem  
                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  i=i+1;
                  if mod(i,100)==0
                      strgg=['Number of Clusters Picking Out : ' num2str(i)];
                      disp(strgg);
                      remaindots=num2str(size(pp,1));
                      strg=['Remained Dots = ' remaindots '/' totaldotsStr];
                      disp(strg);
                  end                  
                  mol(indx1,1:3) = 0;
                  % Revised by Q.Peter Su on Dec 13 2013
                  %  Fast 'precisionstat.m' by set a range of search for 'find' function
                  %  and, remember, DONOT use 'a(x,:)=[];', because 'a(x,:)=0;' is much much faster
                  
                else
                  mol_del=[mol_del; mol(indx1,:)]; 
                  % revised by Q.Peter SU on 19 Sep 2017 for 'Blade shape' problem                  
                  mol(indx1,1:3) = 0;
                  % Revised by Q.Peter Su on Dec 13 2013
                  %  Fast 'precisionstat.m' by set a range of search for 'find' function
                  %  and, remember, DONOT use 'a(x,:)=[];', because 'a(x,:)=0;' is much much faster
                end
                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             else % if the remain dots number is smaller than threshold %%%%
                  %%%% the following comments was made on 19 Sep. 2017 %%%%%
                  %%%% when the molecule density is relatively high, %%%%%%%
                  %%%% some of "blade shapes" appear, this is artifact.%%%%%
                  %%%% So every dots should be calculated local density %%%%
                  %%%% after the deletion, like a re double check ??? %%%%%%
                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  mol(j,1:3)= 0;
                  mol_del=[mol_del; mol(indx1,:)];
                  % revised by Q.Peter SU on 19 Sep 2017 for 'Blade shape' problem  
             end
             pp=find(mol(:,1)>0);
        end
    end
end
toc

%% Re-check the deleted dots
% revised by Q.Peter SU on 19 Sep 2017 for 'Blade shape' problem  

stroo=['There are -' num2str(size(mol_del,1)) ' -Dots been deleted, double check them now ...'];
disp('-----------------------------------------------------------')
disp(stroo);

tic
mol_temp=Data;
Double_check_dot_no=0;
totaldotsStr2=num2str(size(mol_del,1));
dots2plot_temp=[];
for i=1:size(mol_del,1)
    clusterindx = find(abs(mol_del(i,1)-mol_temp(:,1))<dx &...
                       abs(mol_del(i,2)-mol_temp(:,2))<dy &...
                       abs(mol_del(i,3)-mol_temp(:,3))<dz);
    % find the distances between the No.i deleted dots and all the dots
    if mod(i,2000)==0
        remaindots=num2str(size(mol_del,1)-i);
        strg=['Remained Dots for Double Check = ' remaindots '/' totaldotsStr2];
        disp(strg);
    end
    % if the number of dots found within dx/dy/dz are larger than minclustersize,
    % then go ahead for check
    if length(clusterindx)>minclustersize
        Double_check_dis=(mol_del(i,1)-(mol_temp(clusterindx,1))).^2+...
                         (mol_del(i,2)-(mol_temp(clusterindx,2))).^2+...
                         (mol_del(i,3)-(mol_temp(clusterindx,3))).^2/1000000;
        Double_check_indx= find(Double_check_dis==0);
        Double_check_dis(Double_check_indx)=[];
        clusterindx(Double_check_indx)=[];
        Double_check_indx= Double_check_dis==min(Double_check_dis); 
        % find the dot nearest to the No.i dot, but there might be multiple
        % dots with the same nearest distance.
        Double_check_cluster_no=mol_temp(clusterindx(Double_check_indx),13);
        Double_check_cluster_no=Double_check_cluster_no(1);
        % sometimes, code will find multiple dots with same nearest distance.
        if Double_check_cluster_no~=0 && min(Double_check_dis)<dx
           % if the cluster index found is not 0, and the distance < dx,
           % then go ahead, add the deleted one into the dots2plot
            dots2plot_temp=[dots2plot_temp; mol_del(i,1:12)];
            cluster_temp=[mol_del(i,1)*1000, mol_del(i,2)*1000, mol_del(i,3), mol_del(i,12), mol_del(i,4)];
                              % X in nm          Y in nm            Z in nm      frame no.       intensity
            clusters{Double_check_cluster_no,1}=[clusters{Double_check_cluster_no,1};cluster_temp];
            Double_check_dot_no=Double_check_dot_no+1;
        end
    end
end
toc
dots2plot=[dots2plot; dots2plot_temp];
%% Create figure

% X Y axis limit
Xlimmin=min([min(x) min(y)]);
Xlimmax=max([max(x) max(y)]);
Ylimmin=min([min(x) min(y)]);
Ylimmax=max([max(x) max(y)]);

    disp('---------------------------------------------------------------')
    disp('Plot Filtered Localizations and Clusters ...ing...')
close all
    scrsz = get(0,'ScreenSize');
    figure1= figure('OuterPosition',[0.05*scrsz(3) 0.05*scrsz(4) 0.7*scrsz(3) 0.8*scrsz(4)],'Color',[0 0 0]);

% Create axes
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,[Xlimmin Xlimmax]);
ylim(axes1,[Ylimmin Ylimmax]);
hold(axes1,'all');

% Create plot3
plot3(x,y,z,'Marker','.','LineStyle','none','Color',[1 1 0]);
% Create title
surviveratio=100*length(dots2plot)/totaldots;
 title({[Afname(1:end-4)];...
       '   ';...
       ['total dots = ' num2str(totaldots,'%10.1d\n')];...
       ['mean Intensity = ' num2str(meanIall,'%10.1d\n')];...
       ['filter survive ratio = ' num2str(surviveratio,3) '%']},...
       'FontWeight','bold','FontSize',10,'FontName','Arial',...
       'Color',[1 1 1]);
   
% Create axes
axes2 = axes('Parent',figure1,'ZColor',[1 1 1],...
    'YDir','reverse',...
    'YColor',[1 1 1],... 
    'XColor',[1 1 1],...
    'Position',[0.26 0.09 0.9 0.66],...
    'DataAspectRatio',[1 1 1],...
    'Color',[0 0 0]);
xlim(axes2,[Xlimmin Xlimmax]);
ylim(axes2,[Ylimmin Ylimmax]);
hold(axes2,'all');

plot3(dots2plot(:,1),dots2plot(:,2),dots2plot(:,3)/1000,...
      'Marker','.','LineStyle','none','Color',[0 1 0]);
plot3(precision(:,1),precision(:,3),precision(:,5),...
    'MarkerFaceColor',[1 0 0],'MarkerSize',4,'Marker','o',...
    'LineStyle','none',...
    'Color',[1 0 0]);

xres=mean(precision(:,2));
yres=mean(precision(:,4));
zres=mean(precision(:,6));
xyres=mean(precision(:,7));
dotspercluster=mean(precision(:,8));
meanInt=mean(precision(:,9));
% Create title
title({['stdx = ' num2str(xres,3) '   stdy = ' num2str(yres,3)]; ...
     ['stdxy = ' num2str(xyres,3) '   stdz = ' num2str(zres,3)]; ...
     ['       cluster size threshold = ' num2str(minclustersize,2)];...
     ['       total dots = ' num2str(length(dots2plot), '%10.1d\n')];...
     ['       dots per cluster = ' num2str(dotspercluster,2)];...
     ['       cluster number = ' num2str((i-1), '%10.1d\n')];...
     ['       mean Intensity = ' num2str(meanInt, '%10.1d\n')];},...
     'FontWeight','bold','FontSize',10,'FontName','Arial',...
      'Color',[1 1 1]);
 
%% Save files
% revised by qpsu on sep 25 2012
% revised by qpsu on Oct 10 2012
if decision=='a' 
    saveas(gcf, [Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
        '-Filt-Quant.fig'], 'fig');
    saveas(gcf, [Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
        '-Filt-Quant.bmp'], 'bmp');
    save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
        '-Filt-Prec.mat'],'precision');
    save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
        '-Filt-cluster.mat'],'clusters');
    save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
        '-Filt-DAT.dat'], 'dots2plot', '-ascii', '-tabs');
else 
   saveas(gcf, [Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
          '-Filt-manualpick-quant.fig'], 'fig');
   saveas(gcf, [Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
          '-Filt-manualpick-quant.bmp'], 'bmp');
   save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
        '-Filt-precisionstat-manual.mat'],'precision');
   save([Afname(1:end-4) '-Thrd' num2str(minclustersize,2)...
        '-Filt-clusterspicked-manual.mat'],'clusters');

end
toc


% judge=input('wanna fastmap ? y or n - ','s');
% if judge=='y'
% ST08_1_fastmap()
% end 
end

function [precision, clusters, dots2plot]=clusterstat(pos,x,y,z,Intensity)
      clusters=cell(size(pos,1),1);
      precision=[];
      dots2plot=[];
      for j=1:size(pos,1)
        xleft=pos(j,1);
        xright=pos(j,1)+pos(j,3);
        ylow=pos(j,2);
        yhigh=pos(j,2)+pos(j,3);
        clusteridx=(x>xleft&x<xright&y>ylow&y<yhigh);
        cluster_x=x(clusteridx);
        cluster_y=y(clusteridx);
        cluster_z=z(clusteridx);
        cluster_I=Intensity(clusteridx);
        numdots=length(cluster_x);
        meanx=mean(cluster_x); %um
        meany=mean(cluster_y); %um
        meanz=mean(cluster_z); %um
        stdx=std(cluster_x)*1000; %nm
        stdy=std(cluster_y)*1000; %nm
        stdz=std(cluster_z); %nm
        meanI=mean(cluster_I);
        precision=[precision; [meanx stdx meany stdy meanz stdz sqrt(stdx^2+stdy^2) numdots meanI]];
        clusters{j,1}=[cluster_x, cluster_y, cluster_z];
        dots2plot=[dots2plot; [cluster_x, cluster_y, cluster_z]];        
      end      
end