%% Filter Higher Ordered Domain of STORM clusters or dots in 3D
%% Wrote by Q. Peter Su on Oct 26 2012
%  based on Nearestneighbour.m and AlphaHull.m
%  higher ordered domain means cluster of clusters or fibers of clusters
% 
%% Input:
%       1. precision *.mat file for filtered clusters output by precision
%        or filtered *.dat file for filtered   dots   output by precision
%       2. distance between clusters (or dots) in a single domain
%       3. threshold numbers of clusters or dots in a single domain
%       4. clusterpicked*.mat file for localizations/dots in clusters
%
%% Output:
%       1. no. of input clusters or dots
%       2. total no. of domains
%          domains is a cell: Row 1 : clusters X Y Z coordinate in a domain
%                                      X in um. Y in um. Z in um.
%                             Row 2 : clusters number in a domain
%                             Row 3 : domain size of clusters in [X Y Z] in um
%                             Row 4 : domain area in square pixel
%                             Row 5 : reserved domain volume in cubic pixel
%                             Row 6 : Slope & Intercept of the clusters in domain
%                             Row 7 : minX=0 minY=0
%                                     Rotated dots/clusters in domain
%                             Row 8 : rotated domain size of clusters in [X Y Z] in um
%
%                             Row 11: Dots XYZ in domains
%                             Row 12: Dots number in a domain
%                             Row 13: domain size of dots in [X Y Z] in um
%                             Row 17: rotated dots XYZ coordinate
%       3. higher ordered domains number within the threshold
%          scatter3 all the domains within the limit.
%       4.1 Domain #D       *.fig   file of all domained clusters
%       4.2 Domain2D-color  *.fig   file of all domained clusters
%       4.3 Domain #D Hist  *.fig   file of all domains in X Y Z range
%                                           in domain slope
%                                           in cluster no and in area 
%       4.4 Domain #D       *.dat   file of all domained clusters in X Y Z
%                             Column 1 - X in um
%                             Column 2 - Y in um
%                             Column 3 - Z in um
%                             Column 4 - Domain count no.
%                             Column 5 - Domain count layer
%                             Column 6 - Cluster No. in precisionstat.m
%                             Column 7 - Number of Clusters in this domain 
%       4.5 DomainDots      *.dat   file of all dots in domained clusters 
%                                        for fastmap.m rendering
%                             Column 1 - X in um
%                             Column 2 - Y in um
%                             Column 3 - Z in um
%                                  
%% Revised by Q.Peter Su on Oct 27 2012
%  avoid repeat domain counting
%% Revised by Q.Peter Su on Nov 25 2012
%  1.             size(domains{1,k},1)>=domainsize
%     Don't use   length(domains{1,k})>=domainsize
%  2. Histogram title ' X mean - ... '
%                     ' X  std - ... '    
%  3. Add dimensionjudge 2D or 3D
%% Revised by Q.Peter Su on Dec 03 2012
%  New Matlab codes for calculate the area of the domain in squared pixel
%  domains{4,j}=ST40_DomainArea(domains{1,j},25); 
%% Revised by Q.Peter Su on Dec 26 2012
%  Use dendrogram (tree diagram) to calculate the domain
%  talked with Qi Ma
%% Revised by Q.Peter Su on Dec 28 2012
%  nXrange=10*ceil(max(HistXrange)-min(HistXrange))/std(HistXrange);
%  ref to Rahul's Merge_and_plotfiles ===> 10*ceil((max-min)/std)
%  subplot(3,2,1);hist(HistXrange,nXrange);
%% Revised by Q.Peter Su on Jan 02 2013
%  dotsjudge=input('Do you want to pick out the dots in these clusters?','s');
%  generate a *.dat file for those localizations/dots in domained clusters
%% Revised by Q.Peter Su on Jan 05-07 2013
%  Preview the domained clusters in red and all clusters in white
%  fastmap.m the dots or clusters
%% Revised by Q.Peter Su on Jan 09 2013
%  Input global variable 'distancenm' and 'dimensionjudge'
%% Revised by Q.Peter Su on Jan 13-15 2013
%  Rotate all the domains due to slope of clusters in the domain
%  Linear fit : p=polyfit(x,y,1);
%  Rotate Degree : D=atan(p(1));
%  Rotate Vector A=[cos(D) -sin(D) ; sin(D) cos(D)];
%  Display all the good domains by order of [number of clusters in domain]
%% Revised by Q.Peter Su on Jan 15 2013
%  colormap('hsv')
%% Revised by Q.Peter Su on Jan 18 2013
%  Analyze a Single Domain with Surface by calling AlphaHull.m
%  http://www.mathworks.cn/matlabcentral/fileexchange/32725-fast-alpha-hulls-parfor-enabled
%  with many triangles on the surface of a domain.
%% Revised by Q.Peter Su on Jan 20 2013
%  Row 11: Dots XYZ in domains 
%  Row 12: Dots number in a domain
%  Row 13: domain size of dots in [X Y Z] in um
%  Row 17: rotated dots XYZ coordinate
%% Revised by Q.Peter Su on Jan 20 2013
%  try
%      AlphaHull.m a domain
%  catch ME
%      disp('This Domain Cannot be Analysed by the Code, Soryy ! ');
%  end
%% revised by Q.Peter SU on Dec 08 2015  @ BIOPIC, PKU
%   add a black dots for each domain

%% Revised by Q.Peter SU on 07th July 2017  @ IBMD, UTS
%  fix the Z=0 problem
%  add gaussian fit for histogram

%% Revised by Q.Peter SU on 21st July 2017 @ IBMD, UTS
%  try domain size range, like [min max]-[3 10], rather than single threshold
%
%% Revised by Q.Peter SU on 27th Feb 2020 at CPC, USYD
%  Think about a output for all the data points after the batch analysis
%
%   data_output_temp=[];
%   data_output_temp=[HistXrange(2,:)', HistYrange(2,:)'];
%   data_output_temp(:,3)=(data_output_temp(:,1) + data_output_temp(:,2))/2;
%   data_output_temp(:,4)=Nearest_Dis;
%   data_output=[data_output; data_output_temp];
% 
%   end  % for file_i
%   matname2=strcat('0-diary-',datestr(now, 30),'.mat');
%   save([matname2],'data_output');
%   disp(matname2);
%  
%  Data_Output.mat
%   Column 1 - X range of domain in um
%   Column 2 - Y range of domain in um
%   Column 3 - Diameter of domain in um
%   Column 4 - NearestNeighbour distance of domains in um
%


%%
function ST10_1_Domain3D(distancenm,dimensionjudge)
%% clear
clc
%clear all;
close all;
%% 1. Data import
[FileName PathName] = uigetfile({'*prec*.mat','Precision Cluster *.mat file';...
                         '*.mat','All mat files';...
                         '*filter*.dat','Filtered Dots *.dat file';...
                         '*.dat','All dat files'},...
               'Choose the Filtered-precision Clusters *.mat file.',...
               'Multiselect', 'off');
disp('... ...');
cd(PathName);
disp('Cluster Center Coordinate File : ');disp(PathName);disp(FileName);
[FileName2 PathName2] = uigetfile({'*clus*.mat','Cluster *.mat file';...
                         '*.mat','All mat files'},...
               'Choose the Filtered-precision Clusters *.mat file.',...
               'Multiselect', 'off');
disp('... ...');
disp('Dots Coordinate File : ');disp(PathName2);disp(FileName2);
disp('Import Data ...ing... ')
tic
images = importdata([PathName FileName]);

if FileName2==0
    dots=cell([size(images,1),1]); % Revised by Q.Peter Su on Jan 20 2013
else
    dots = importdata([PathName2 FileName2]);
end

cd(PathName);
toc
%% change matrix due to different file types
% revised on Oct 27 2012
if strcmp(FileName(end-3:end),'.mat')
    centers=[]; 
    centers(:,1)=images(:,1);   % cluster center X    in um
    centers(:,2)=images(:,3);   % cluster center Y    in um
    centers(:,3)=images(:,5);   % cluster center Z    in um
    centers(:,6)=images(:,10);  % cluster no. output by precisionstat
                                % Revised by Q. Peter Su
                                % on Jan 02 2013
elseif strcmp(FileName(end-3:end), '.dat')
    centers=[]; 
    centers(:,1)=images(:,1);        % dots center X    in um
    centers(:,2)=images(:,2);        % dots center Y    in um
    centers(:,3)=images(:,3)/1000;   % dots center Z    in um
    centers(:,6)=(1:size(images,1))';
end
strout=strcat('Totally, There are -', num2str(size(centers,1)) , '- Dots/Clusters.');
disp(strout);
% sort centers by x position if needed
% mat file and dat file output by precision is already sorted by x position
[Y,I]=sort(centers(:,1));
centers=centers(I,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% centers Column-1 Column-2 Column-3  Column-4   Column-5   Column-6   %%%
%%%             X       Y        Z     domaincount   layer   cluster no. %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin<1
distancenm    =input('Distance between Dots in a Single Domain (in nm). - ');
end
distance=distancenm/1000;
domaincount=0;
if nargin<1
dimensionjudge=input('You Wanna Analyse the Domian in 2D or 3D?        - ');
end
centerstemp=centers;
%% Mainbody Two-Dimension 2D domain count
disp('------------------------------------------------------------------')
disp('Calculate Domain ...ing...');
tic
if dimensionjudge==2
while min(centers(:,4))==0
    %%%%%%%%%%%%%%%% New Domaincount %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    domaincount=domaincount+1;  % display the domaincount
    if mod(domaincount,100)==0
    strdomain=['Domain Count : '  num2str(domaincount)];
    disp(strdomain);
    end
    newdomain=find(centers(:,4)==0);
    centers(newdomain(1),4)=domaincount; % New domain 1st center
    layercount=1;
    centers(newdomain(1),5)=layercount;   
    %%%%%%%%%%%%%%% New Domaincount 1st layer finding %%%%%%%%%%%
    centerstemp(newdomain(1),1:3)=[100000 100000 100000];
    dis=0;
    while dis<=distance
        IDX=nearestneighbour(centers(newdomain(1),1:2)', centerstemp(:,1:2)');
        dis=sqrt((centerstemp(IDX,1)-centers(newdomain(1),1)).^2+...
                 (centerstemp(IDX,2)-centers(newdomain(1),2)).^2);
        if dis<=distance
           centerstemp(IDX,1:3)=[100000 100000 100000];
           centers(IDX,4)=domaincount;
           centers(IDX,5)=layercount+1;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    DomainFindIndex=find(centers(:,4)==domaincount & centers(:,5)==2);
    while ~isempty(DomainFindIndex)
        %%%%%%% find all the nearestneighbours < distance threshold
        %%%%%%% for specific Domaincount & layercount  
        for i=1:length(DomainFindIndex) % scan all the clusters in the centers
            centerstemp(DomainFindIndex(i),1:3)=[100000 100000 100000];
            dis=0;
            while dis<=distance
            IDX=nearestneighbour(centers(DomainFindIndex(i),1:2)', centerstemp(:,1:2)');
            dis=sqrt((centerstemp(IDX,1)-centers(DomainFindIndex(i),1)).^2+...
                     (centerstemp(IDX,2)-centers(DomainFindIndex(i),2)).^2);
            if dis<=distance
               centerstemp(IDX,1:3)=[100000 100000 100000];
               centers(IDX,4)=domaincount;
               centers(IDX,5)=layercount+1;
            end
            end
        end
        layercount=layercount+1;
        DomainFindIndex=find(centers(:,4)==domaincount ...
                     & centers(:,5)==layercount); 
     end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
strdomain=['Domain Count : '  num2str(domaincount)];
disp(strdomain);
%% Mainbody Three-Dimension 3D domain count
%peter=1
if dimensionjudge==3
while min(centers(:,4))==0
    %%%%%%%%%%%%%%%% New Domaincount %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    domaincount=domaincount+1;  % display the domaincount
    if mod(domaincount,100)==0
    strdomain=['Domain Count : '  num2str(domaincount)];
    disp(strdomain);
    end
    newdomain=find(centers(:,4)==0);
    centers(newdomain(1),4)=domaincount; % New domain 1st center
    layercount=1;
    centers(newdomain(1),5)=layercount;   
    %%%%%%%%%%%%%%% New Domaincount 1st layer finding %%%%%%%%%%%
    centerstemp(newdomain(1),1:3)=[100000 100000 100000];
    dis=0;
    while dis<=distance
        IDX=nearestneighbour(centers(newdomain(1),1:3)', centerstemp(:,1:3)');
        dis=sqrt((centerstemp(IDX,1)-centers(newdomain(1),1)).^2+...
                 (centerstemp(IDX,2)-centers(newdomain(1),2)).^2+...
                 (centerstemp(IDX,3)-centers(newdomain(1),3)).^2);
        if dis<=distance
           centerstemp(IDX,1:3)=[100000 100000 100000];
           %peter=peter+1
           centers(IDX,4)=domaincount;
           centers(IDX,5)=layercount+1;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    DomainFindIndex=find(centers(:,4)==domaincount & centers(:,5)==2);
    while ~isempty(DomainFindIndex)
        %%%%%%% find all the nearestneighbours < distance threshold
        %%%%%%% for specific Domaincount & layercount  
        for i=1:length(DomainFindIndex) % scan all the clusters in the centers
            centerstemp(DomainFindIndex(i),1:3)=[100000 100000 100000];
            dis=0;
            while dis<=distance
            IDX=nearestneighbour(centers(DomainFindIndex(i),1:3)', centerstemp(:,1:3)');
            dis=sqrt((centerstemp(IDX,1)-centers(DomainFindIndex(i),1)).^2+...
                     (centerstemp(IDX,2)-centers(DomainFindIndex(i),2)).^2+...
                     (centerstemp(IDX,3)-centers(DomainFindIndex(i),3)).^2);
            if dis<=distance
               centerstemp(IDX,1:3)=[100000 100000 100000];
               %peter=peter+1
               centers(IDX,4)=domaincount;
               centers(IDX,5)=layercount+1;
            end
            end
        end
        layercount=layercount+1;
        DomainFindIndex=find(centers(:,4)==domaincount ...
                     & centers(:,5)==layercount); 
     end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
strdomain=['Domain Count : '  num2str(domaincount)];
disp(strdomain);

%% Generate domains cells
domains=cell(20,domaincount);
for j=1:size(centers,1)
    n=centers(j,4); % n - Domain  No.
    m=centers(j,6); % m - Cluster No.
    domains{1,n}  = [domains{1,n} ; centers(j,:)];
    domains{11,n} = [domains{11,n}; dots{m}/1000];
end
%% Cluster number for every domain & Size of every domain in x y z
%clusterindomain=0;
clusterindomains=zeros([1,domaincount]);
Xrange=zeros([1,domaincount]);Yrange=zeros([1,domaincount]);Zrange=zeros([1,domaincount]);
XrangeR=zeros([1,domaincount]);YrangeR=zeros([1,domaincount]);ZrangeR=zeros([1,domaincount]);
DomainArea=zeros([1,domaincount]);%DomainVolume=zeros([1,domaincount]);
DomainSlope=zeros([1,domaincount]);DomainJieju=zeros([1,domaincount]);

for j=1:domaincount
    
    %%%%%%%%% Row 2 : Clusters number of every domain %%%%%%%%%%%%%%%
    clusterindomains(j)=size(domains{1,j},1);
    domains{2,j}=clusterindomains(j);
    
    %%%%%%%%% Row 3 : Xrange in um, Yrange in um, Zrange in um %%%%%%%
    Xrange(j)=max(domains{1,j}(:,1))-min(domains{1,j}(:,1));
    Yrange(j)=max(domains{1,j}(:,2))-min(domains{1,j}(:,2));
    Zrange(j)=max(domains{1,j}(:,3))-min(domains{1,j}(:,3));
    domains{3,j}=[Xrange(j) Yrange(j) Zrange(j)];
    
    %%%%%%%%% Row 4 : Domain Area in Square Pixel %%%%%%%%%%%%%%%%%%%% 
    domains{4,j}=ST40_DomainArea(domains{1,j},25); 
    DomainArea(j)=domains{4,j};
    
    %%%%%%%%% Row 5 : Domain Volume in Cubic Pixels %%%%%%%%%%%%%%%%%%
    %domains{5,j}=ST41_DomainVolume(domains{1,j},25); 
    %DomainVolume(j)=domains{5,j};
    
    %%%%%%%%% Row 6 : (Polyfit,1) fit of the clusters in domain %%%%%%
                
    if size(domains{1,j},1)>=3
        p=polyfit(domains{1,j}(:,1),domains{1,j}(:,2),1);
        domains{6,j}=p;
        DomainSlope(j)=atan(p(1)); % slope in radian
        DomainJieju(j)=p(2);
    else
        domains{6,j}=[0 0];
        DomainSlope(j)=0;
    end
    
    
    if DomainSlope(j)>=0        
        domains{7,j}(:,1)=domains{1,j}(:,1);
        domains{7,j}(:,2)=domains{1,j}(:,2)-DomainJieju(j);
        domains{7,j}(:,3)=domains{1,j}(:,3)-min(domains{1,j}(:,3));
    elseif DomainSlope(j)<0
        domains{7,j}(:,1)=domains{1,j}(:,1);
        domains{7,j}(:,2)=domains{1,j}(:,2)-DomainJieju(j);
        domains{7,j}(:,3)=domains{1,j}(:,3)-min(domains{1,j}(:,3));
    end
    
    %%%%%%%% Row 7 : Rotated dots/clusters in the domain %%%%%%%%%%%%%
    %%%%%%%%         Min(X)=0, Min(Y)=0                  %%%%%%%%%%%%%
    %%%%%%%%         Rotated by -atan(DomainSlope)       %%%%%%%%%%%%%
        
    
    if size(domains{1,j},1)>=3 && ~isnan(p(1))        
        if DomainSlope(j)>0
            roangle=(pi/2)+DomainSlope(j); % in radian
            Acluster=[cos(-roangle) -sin(-roangle);sin(-roangle) cos(-roangle)]*...
                     [(domains{7,j}(:,1))'; (domains{7,j}(:,2))' ]; 
            if ~FileName2==0
            Adots  = [cos(-roangle) -sin(-roangle);sin(-roangle) cos(-roangle)]*...
                     [(domains{11,j}(:,1))'; (domains{11,j}(:,2))' ];
            end            
        elseif DomainSlope(j)==0
            roangle=(pi/2); % in radian
            Acluster=[cos(roangle) -sin(roangle);sin(roangle) cos(roangle)]*...
                     [(domains{7,j}(:,1))'; (domains{7,j}(:,2))' ];
            % clockwise rotate the dots for 'roangle' degree 
            if ~FileName2==0
            Adots  = [cos(roangle) -sin(roangle);sin(roangle) cos(roangle)]*...
                     [(domains{11,j}(:,1))'; (domains{11,j}(:,2))' ];
            end            
        elseif DomainSlope(j)<0 
            roangle=(pi/2)+DomainSlope(j); % in radian
            Acluster=[cos(-roangle) -sin(-roangle);sin(-roangle) cos(-roangle)]*...
                     [(domains{7,j}(:,1))'; (domains{7,j}(:,2))' ];    
            if ~FileName2==0
            Adots  = [cos(-roangle) -sin(-roangle);sin(-roangle) cos(-roangle)]*...
                     [(domains{11,j}(:,1))'; (domains{11,j}(:,2))' ];
            end              
        end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% x=linspace(-2,2);y=x.^2; A=[cosd(-45) -sind(-45);sind(-45) cosd(-45)]*[x;y];    %%%%
%%%% plot(x,y,'r',A(1,:),A(2,:));axis equal;legend('原图像','顺时针旋转45°后的图像'); %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        domains{7,j}=[];
        domains{7,j}(:,1:2)=[(Acluster(1,:))' (Acluster(2,:))'];
        domains{7,j}(:,3:6)=domains{1,j}(:,3:6);
        
        Xmin=min(domains{7,j}(:,1));
        Ymin=min(domains{7,j}(:,2));
        Zmin=min(domains{7,j}(:,3));
        
        domains{7,j}(:,1)=domains{7,j}(:,1)-Xmin;
        domains{7,j}(:,2)=domains{7,j}(:,2)-Ymin;
        domains{7,j}(:,3)=domains{7,j}(:,3)-Zmin;
    
    else
        domains{7,j}=[0 0 0 domains{1,j}(1,4:6)];
        Xmin=min(domains{1,j}(:,1));
        Ymin=min(domains{1,j}(:,2));
        Zmin=min(domains{1,j}(:,3));
    end 
    
    %%%%%%%% Row 8 : Rotated Domain Size in X Y Z in um %%%%%%%%%%%%%
    XrangeR(j)=max(domains{7,j}(:,1))-min(domains{7,j}(:,1));
    YrangeR(j)=max(domains{7,j}(:,2))-min(domains{7,j}(:,2));
    ZrangeR(j)=max(domains{7,j}(:,3))-min(domains{7,j}(:,3));
    domains{8,j}=[XrangeR(j) YrangeR(j) ZrangeR(j)];
    
    
    
    
    if ~FileName2==0 
    %%%%%%%% Row 11 : Dots X Y Z in domains %%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%% Row 12 : Dots Number in domains %%%%%%%%%%%%%%%%%%%%%%%%
        domains{12,j}=size(domains{11,j},1);
    %%%%%%%% Row 13 : Dots range X Y Z in domains %%%%%%%%%%%%%%%%%%%
        if domains{12,j}~=0
            domains{13,j}=[max(domains{11,j}(:,1))-min(domains{11,j}(:,1)),...
                           max(domains{11,j}(:,2))-min(domains{11,j}(:,2)),...
                           max(domains{11,j}(:,3))-min(domains{11,j}(:,3))];
        end
    %%%%%%%% Row 14 : Dots Area   in domains %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%% Row 15 : Dots Volume in domains %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%% Row 16 :             in domains %%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %%%%%%%% Row 17 : Rotated Dots X Y Z %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if size(domains{1,j},1)>=3
            if domains{12,j}~=0
                domains{17,j}(:,1:2)= [(Adots(1,:))' (Adots(2,:))'];
                domains{17,j}(:,3)  = domains{11,j}(:,3);
        
                domains{17,j}(:,1)=domains{17,j}(:,1)-Xmin;
                domains{17,j}(:,2)=domains{17,j}(:,2)-Ymin;
                domains{17,j}(:,3)=domains{17,j}(:,3)-Zmin;
            end
        elseif size(domains{1,j},1)==1
            domains{17,j}=domains{11,j};
            if domains{12,j}~=0
                domains{17,j}(:,1)=domains{17,j}(:,1)-Xmin;
                domains{17,j}(:,2)=domains{17,j}(:,2)-Ymin;
                domains{17,j}(:,3)=domains{17,j}(:,3)-Zmin;
            end
        end  
   %%%%%%%% Row 18 : Domain Center X Y Z %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        domains{18,j}=[mean(domains{1,j}(:,1)) ...
                       mean(domains{1,j}(:,2)) ...
                       mean(domains{1,j}(:,3)) ];
    end    
end
toc
%% Preview the domained clusters with the same color
%   Revised by Q. Peter Su on Jan 04 2013
%   to speed up the process of domained clusters
disp('------------------------------------------------------------------')
judge='n';%input('Preview a Threshold of Domained Clusters? y or n - ','s');
          % de-activate by Q.Peter SU on Dec 08 2015
if judge=='y'
%     close all
%     judge='n';
%     while judge=='n'    
%     scrsz = get(0,'ScreenSize');
%     figure('OuterPosition',[1 0.5*scrsz(4) 0.75*scrsz(4) 0.5*scrsz(4)],'Color',...
%         [0.313725501298904 0.313725501298904 0.313725501298904]);
%     hdl1=gcf;    clf;
%     figure('OuterPosition',[1 1 0.75*scrsz(4) 0.5*scrsz(4)]);
%     hdl2=gcf;    clf;
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     domainsize=input('Type the Domain Size Threshold (Clusters/Domain) - ');
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%     figure(hdl1)
%     axes1 = axes('Parent',figure(hdl1),'ZColor',[1 1 1],...
%         'YDir','reverse','YColor',[1 1 1],...
%         'XColor',[1 1 1],...
%         'Color',[0 0 0]);
%     box(axes1,'on');    hold(axes1,'all');
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     DomainIndex=find(clusterindomains>=domainsize);    
%     higherdomain=length(DomainIndex);    
%     domainsDAT=[];
%     
%     HistXrange=zeros([2,higherdomain]);
%     HistYrange=zeros([2,higherdomain]);
%     HistZrange=zeros([2,higherdomain]);
%     HistClusterNo=zeros([2,higherdomain]);
%     HistDomainSlope=zeros([2,higherdomain]);
%     HistDomainArea=zeros([2,higherdomain]);
%     
%     xplotmin=min(centers(:,1)); xplotmax=max(centers(:,1));
%     yplotmin=min(centers(:,2)); yplotmax=max(centers(:,2));
%     
%     disp('... ...');
%     disp('Analyze Domain ...ing...');
%     disp('... ...');
%     
%     for j=1:higherdomain
%     %%%%%%% Data for Domains DAT %%%%%%%%%%%%%%%%%%%%%%
%         k=DomainIndex(j);
%         domainstemp=[domains{1,k}(:,1:6) repmat(domains{2,k},domains{2,k},1)];
%         domainsDAT=[domainsDAT; domainstemp];
%     %%%%%%% Data for Histogram %%%%%%%%%%%%%%%%%%%%%%%%%
%         HistXrange(1:2,j)     = [k;Xrange(k)];
%         HistYrange(1:2,j)     = [k;Yrange(k)];
%         HistZrange(1:2,j)     = [k;Zrange(k)];
%         HistClusterNo(1:2,j)  = [k;domains{2,k}];
%         HistDomainArea(1:2,j) = [k;domains{4,k}];
%         HistDomainSlope(1:2,j)= [k;DomainSlope(k)];
%     end
% 
%     scatter3(centers(:,1),centers(:,2),centers(:,3),2,[1 1 1]);
%     hold on
%     scatter3(domainsDAT(:,1),domainsDAT(:,2),domainsDAT(:,3),3,[1 0 0],'filled');
%     xlim([xplotmin xplotmax]);  ylim([yplotmin yplotmax]);  
%     axis equal
%     view(-60,60)
%     title({[FileName(1:end-4)];...
%            ['Total Clusters Numbers = ' num2str(length(images))];...
%            ['Distance between Clusters in one Domain = ' num2str(distancenm) 'nm'];...
%            ['Total Domains Numbers = ' num2str(domaincount)];...
%            ['Threshold of Clusters per Domain = ' num2str(domainsize)];...
%            ['Higher Ordered Domain No. = ' num2str(higherdomain)]},...
%            'FontWeight','bold','FontSize',10,'FontName','Arial',...
%            'Color',[1 1 1]);
%    
%     if dimensionjudge==2
%         fig1name=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
%                         '-Thrd',num2str(domainsize),'-Domain2D-Single.fig');
%         saveas(gcf, [fig1name], 'fig');
%     elseif dimensionjudge==3
%         fig1name=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
%                         '-Thrd',num2str(domainsize),'-Domain3D-Single.fig');
%         saveas(gcf, [fig1name], 'fig');
%     end
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     nXrange=10*ceil((max(HistXrange(2,:))-min(HistXrange(2,:)))/std(HistXrange(2,:)));
%     nYrange=10*ceil((max(HistYrange(2,:))-min(HistYrange(2,:)))/std(HistYrange(2,:)));
%     nZrange=10*ceil((max(HistZrange(2,:))-min(HistZrange(2,:)))/std(HistZrange(2,:)));
%     nclusterno=10*ceil((max(HistClusterNo(2,:))-min(HistClusterNo(2,:)))/std(HistClusterNo(2,:)));
%     ndomainslope=10*ceil((max(HistDomainSlope(2,:))-min(HistDomainSlope(2,:)))/std(HistDomainSlope(2,:)));
%     ndomainarea=10*ceil((max(HistDomainArea(2,:))-min(HistDomainArea(2,:)))/std(HistDomainArea(2,:)));
%     %%%%%%%%%%%%%%%%%%%%%  Revised by Q.Peter Su on Dec 28 2012    %%%%%%
%     figure(hdl2)
%           subplot(3,2,1);hist(HistXrange(2,:),nXrange);
%           xlabel('Xrange of Domains(um)','FontWeight','bold');
%           ylabel('Occurrence');
%           title({['X mean - ', num2str(mean(HistXrange(2,:)))];...
%                  ['X std - ',num2str(std(HistXrange(2,:)))]});
%           grid on;zoom on;
%       
%           subplot(3,2,2);hist(HistYrange(2,:),nYrange);
%           xlabel('Yrange of Domains(um)','FontWeight','bold');
%           ylabel('Occurrence');
%           title({['Y mean - ', num2str(mean(HistYrange(2,:)))];...
%                  ['Y std - ',num2str(std(HistYrange(2,:)))]});
%           grid on;zoom on;
%       
%           subplot(3,2,3);hist(HistZrange(2,:),nZrange);
%           xlabel('Zrange of Domains(um)','FontWeight','bold');
%           ylabel('Occurrence');
%           title({['Z mean - ', num2str(mean(HistZrange(2,:)))];...
%                  ['Z std - ',num2str(std(HistZrange(2,:)))]});
%           grid on;zoom on;
%           
%           subplot(3,2,4);hist(HistDomainSlope(2,:),ndomainslope);
%           xlabel('Slope of Domains(Radians)','FontWeight','bold');
%           ylabel('Occurrence');
%           title({['Slope mean - ', num2str(mean(HistDomainSlope(2,:)))];...
%                  ['Slope std - ',num2str(std(HistDomainSlope(2,:)))]});
%           grid on;zoom on;
%       
%           subplot(3,2,5);hist(HistClusterNo(2,:),nclusterno);
%           xlabel('Clusters or Dots no. in Domains','FontWeight','bold');
%           ylabel('Occurrence');
%           title({['Cluster mean - ', num2str(mean(HistClusterNo(2,:)))];...
%                  ['Cluster std - ',num2str(std(HistClusterNo(2,:)))]});
%           grid on;zoom on;
%           
%           subplot(3,2,6);hist(HistDomainArea(2,:),ndomainarea);
%           xlabel('Domain Area in Sq. Pixel (25nm * 25nm)','FontWeight','bold');
%           ylabel('Occurrence');
%           title({['DomainArea mean - ', num2str(mean(HistDomainArea(2,:)))];...
%                  ['DomainArea std - ',num2str(std(HistDomainArea(2,:)))]});
%           grid on;zoom on;
%      
%           whitebg([0 0 0]);  
%     if dimensionjudge==2
%         fig2name=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
%                         '-Thrd',num2str(domainsize),'-Domain2D-Hist.fig');
%         saveas(gcf, [fig2name], 'fig');
%     elseif dimensionjudge==3
%         fig2name=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
%                         '-Thrd',num2str(domainsize),'-Domain3D-Hist.fig');
%         saveas(gcf, [fig2name], 'fig');
%     end
%     
%     disp('... ...');
%     disp('Are You Happy with the Previewed Domains? - ')
%     judge=input('y for Yes and Save; n for Another Threshold.     - ','s');
%     if judge=='n'
%         close all
%         delete(fig1name);
%         delete(fig2name);
%     end
%     end
% 
%     if dimensionjudge==2
%         outfilename=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
%                            '-Thrd',num2str(domainsize),'-Domain2D-Dat.dat');
%         save(outfilename, 'domainsDAT', '-ascii', '-tabs');
%     elseif dimensionjudge==3
%         outfilename=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
%                            '-Thrd',num2str(domainsize),'-Domain3D-Dat.dat');
%         save(outfilename, 'domainsDAT', '-ascii', '-tabs');
%     end
end
%% Final Figures of colored clusters
disp('-------------------------------------------------------------------')
judge=input('Do You Wanna See the Colored Domains? y or n     - ','s');
if judge=='y'
    close all
    judge='n';
    while judge=='n'
       close all
       scrsz = get(0,'ScreenSize');
       figure('OuterPosition',[1 0.5*scrsz(4) 0.75*scrsz(4) 0.5*scrsz(4)],'Color',...
          [0.313725501298904 0.313725501298904 0.313725501298904]);
        hdl1=gcf;    clf;
        figure('OuterPosition',[1 1 0.75*scrsz(4) 0.5*scrsz(4)]);
        hdl2=gcf;    clf;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        domainsize=input('Type the Domain Size Threshold (Clusters/Domain) - ');
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        figure(hdl1)
        axes1 = axes('Parent',figure(hdl1),'ZColor',[1 1 1],...
            'YDir','reverse','YColor',[1 1 1],...
            'XColor',[1 1 1],...
            'Color',[0 0 0]);
        box(axes1,'on');    hold(axes1,'all');
        colormap('hsv')
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %%%%%%%%%% Revised by Q.Peter SU on 21st July 2017 %%%%%%%%%%%%%%%%%%
        %  try domain size range, like [3 10], rather than single threshold
        if size(domainsize,2)==1
            DomainIndex=find(clusterindomains>=domainsize);    
        elseif size(domainsize,2)==2
            DomainIndex=find(clusterindomains>=domainsize(1) & clusterindomains<=domainsize(2));
        end
        %%%%%%%%%% Revised by Q.Peter SU on 21st July 2017 %%%%%%%%%%%%%%%%%%
        
        higherdomain=length(DomainIndex);    
        domainsDAT=[];
        domainCenter=[];

        HistXrange=zeros([2,higherdomain]);
        HistYrange=zeros([2,higherdomain]);
        HistZrange=zeros([2,higherdomain]);
        HistClusterNo=zeros([2,higherdomain]);
        HistDomainSlope=zeros([2,higherdomain]);
        HistDomainArea=zeros([2,higherdomain]);   

        
        xplotmax=max(centers(:,1))*1.1;
        xplotmin=min(centers(:,1))-0.1*xplotmax; 
        yplotmax=max(centers(:,2))*1.1;
        yplotmin=min(centers(:,2))-0.1*yplotmax;
        % modified by Q.Peter SU on Dec 08 2015
        

        disp('... ...');
        disp('Analyze Domain ...ing...');
        disp('... ...');

        for j=1:higherdomain   
            k=DomainIndex(j);
            domainstemp=[domains{1,k}(:,1:6) repmat(domains{2,k},domains{2,k},1)];
            domainsDAT=[domainsDAT
                        domainstemp];   
            domainCenter=[domainCenter
                          domains{18,k} 0 0 0];
            HistXrange(1:2,j)     = [k;Xrange(k)];
            HistYrange(1:2,j)     = [k;Yrange(k)];
            HistZrange(1:2,j)     = [k;Zrange(k)];
            HistClusterNo(1:2,j)  = [k;domains{2,k}];
            HistDomainArea(1:2,j) = [k;domains{4,k}];
            HistDomainSlope(1:2,j)= [k;DomainSlope(k)];
        end
        %%%%%%%%%%% Revised by Q.Peter SU on 28th Feb 2020 %%%%%%%%%%
        Nearest_Index=nearestneighbour(domainCenter(:,1:3)');
        Nearest_Dis=sqrt((domainCenter(Nearest_Index,1)-domainCenter(:,1)).^2+...
                         (domainCenter(Nearest_Index,2)-domainCenter(:,2)).^2+...
                         (domainCenter(Nearest_Index,3)-domainCenter(:,3)).^2);
        
        
        scatter3(domainsDAT(:,1),domainsDAT(:,2),domainsDAT(:,3),...
                 5,domainsDAT(:,7),'filled');
             hold on
        scatter3(domainCenter(:,1),domainCenter(:,2),domainCenter(:,3),...
                 8,domainCenter(:,4:6),'filled');
             hold off
        xlim([xplotmin xplotmax]);  ylim([yplotmin yplotmax]); 
        axis equal
        grid off
        view(0,90)
        title({[FileName(1:end-4)];...
               ['Total Clusters Numbers = ' num2str(length(images))];...
               ['Distance between Clusters in one Domain = ' num2str(distancenm) 'nm'];...
               ['Total Domains Numbers = ' num2str(domaincount)];...
               ['Threshold of Clusters per Domain = ' num2str(domainsize)];...
               ['Higher Ordered Domain No. = ' num2str(higherdomain)];...
               ['Colormap from Red to Blue means cluster number from many to few']},...
               'FontWeight','bold','FontSize',10,'FontName','Arial',...
               'Color',[1 1 1]);
            colorbar('peer',axes1);
        if dimensionjudge==2
            fig1name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain2D-Color.fig');
            fig1name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain2D-Color.bmp');
            saveas(gcf, [fig1name1], 'fig');
            saveas(gcf, [fig1name2], 'bmp');
        elseif dimensionjudge==3
            fig1name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain3D-Color.fig');
            fig1name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain3D-Color.bmp');
            saveas(gcf, [fig1name1], 'fig');
            saveas(gcf, [fig1name2], 'bmp');
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        nXrange=10*ceil((max(HistXrange(2,:))-min(HistXrange(2,:)))/std(HistXrange(2,:)));
        nYrange=10*ceil((max(HistYrange(2,:))-min(HistYrange(2,:)))/std(HistYrange(2,:)));
        nZrange=10*ceil((max(HistZrange(2,:))-min(HistZrange(2,:)))/std(HistZrange(2,:)));
        nclusterno=10*ceil((max(HistClusterNo(2,:))-min(HistClusterNo(2,:)))/std(HistClusterNo(2,:)));
        ndomainslope=10*ceil((max(HistDomainSlope(2,:))-min(HistDomainSlope(2,:)))/std(HistDomainSlope(2,:)));
        ndomainarea=10*ceil((max(HistDomainArea(2,:))-min(HistDomainArea(2,:)))/std(HistDomainArea(2,:)));
        %%%%%%%%%%%%%%%%%%%%%  Revised by Q.Peter Su on Dec 28 2012    %%%%%%
        figure(hdl2)
          subplot(3,2,1);hist(HistXrange(2,:),nXrange);
          xlabel('Xrange of Domains(um)','FontWeight','bold');
          ylabel('Occurrence');
          title({['X mean - ', num2str(mean(HistXrange(2,:)))];...
                 ['X std - ',num2str(std(HistXrange(2,:)))]});
          grid on;zoom on;

          subplot(3,2,2);hist(HistYrange(2,:),nYrange);
          xlabel('Yrange of Domains(um)','FontWeight','bold');
          ylabel('Occurrence');
          title({['Y mean - ', num2str(mean(HistYrange(2,:)))];...
                 ['Y std - ',num2str(std(HistYrange(2,:)))]});
          grid on;zoom on;
          
          %%%%%%%%%%% Revised by Q.Peter SU on 07th July 2017 %%%%%%%%%%
          if max(HistZrange(2,:))==0
              %%%%%%%%%%% Revised by Q.Peter SU on 28th Feb 2020 %%%%%%%%%%
              subplot(3,2,3);hist(Nearest_Dis, nXrange);
              xlabel('NearestNeighbour Distance of Domains(um)','FontWeight','bold');
              ylabel('Occurrence');
              title({['NearestNeighbour Dis mean - ', num2str(mean(Nearest_Dis))];...
                     ['NearestNeighbour Dis std - ',  num2str(std(Nearest_Dis))]});
              grid on;zoom on;
              % subplot(3,2,3);hist(HistZrange(2,:));
              % xlabel('Zrange of Domains(um)','FontWeight','bold');
              % ylabel('Occurrence');
              % title({['Z mean - ', num2str(mean(HistZrange(2,:)))];...
              %        ['Z std - ',num2str(std(HistZrange(2,:)))]});
              % grid on;zoom on;
          else
              subplot(3,2,3);hist(HistZrange(2,:),nZrange);
              xlabel('Zrange of Domains(um)','FontWeight','bold');
              ylabel('Occurrence');
              title({['Z mean - ', num2str(mean(HistZrange(2,:)))];...
                  ['Z std - ',num2str(std(HistZrange(2,:)))]});
              grid on;zoom on;
          end

          subplot(3,2,4);hist(HistDomainSlope(2,:),ndomainslope);
          xlabel('Slope of Domains(Radians)','FontWeight','bold');
          ylabel('Occurrence');
          title({['Slope mean - ', num2str(mean(HistDomainSlope(2,:)))];...
                 ['Slope std - ',num2str(std(HistDomainSlope(2,:)))]});
          grid on;zoom on;

          subplot(3,2,5);hist(HistClusterNo(2,:),nclusterno);
          xlabel('Clusters or Dots no. in Domains','FontWeight','bold');
          ylabel('Occurrence');
          title({['Cluster mean - ', num2str(mean(HistClusterNo(2,:)))];...
                 ['Cluster std - ',num2str(std(HistClusterNo(2,:)))]});
          grid on;zoom on;

          subplot(3,2,6);hist(HistDomainArea(2,:),ndomainarea);
          xlabel('Domain Area in Sq. Pixel (25nm * 25nm)','FontWeight','bold');
          ylabel('Occurrence');
          title({['DomainArea mean - ', num2str(mean(HistDomainArea(2,:)))];...
                 ['DomainArea std - ',num2str(std(HistDomainArea(2,:)))]});
          grid on;zoom on;

          whitebg([0 0 0]);

        if dimensionjudge==2
            fig2name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain2D-Hist.fig');
            fig2name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain2D-Hist.bmp');
            saveas(gcf, [fig2name1], 'fig');
            saveas(gcf, [fig2name2], 'bmp');
        elseif dimensionjudge==3
            fig2name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain3D-Hist.fig');
            fig2name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                        '-Thrd',num2str(domainsize),'-Domain3D-Hist.bmp');
            saveas(gcf, [fig2name1], 'fig');
            saveas(gcf, [fig2name2], 'bmp');
        end

        disp('... ...');
        disp('Are you happy with the domains?  ')
        judge=input('y for Yes and save; n for Another Threshold.     - ','s');
        if judge=='n'
            close all
            delete(fig1name1);  delete(fig2name1);
            delete(fig1name2);  delete(fig2name2);
            
        end
    end

    if dimensionjudge==2
        outfilename=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain2D-Dat.dat');
        save(outfilename, 'domainsDAT', '-ascii', '-tabs');
        outfilenameMAT=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain2D-Mat.mat');
        save (outfilenameMAT,'-mat');
    elseif dimensionjudge==3
        outfilename=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain3D-Dat.dat');
        save(outfilename, 'domainsDAT', '-ascii', '-tabs');
        outfilenameMAT=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain3D-Mat.mat');
        save (outfilenameMAT,'-mat');
    end
    disp(fig1name1);    disp(fig1name2);
    disp(fig2name1);    disp(fig2name2);
    disp(outfilename);  disp(outfilenameMAT);
end
%% Pick out the dots in the domained clusters
disp('-------------------------------------------------------------------')
dotsjudge=input('Wanna to pick out the dots in these clusters?    - ','s');
if dotsjudge=='y'
    % dots = importdata([PathName2 FileName2]);
    %        Column 1 : X in nano-meter
    %        Column 2 : Y in nano-meter
    %        Column 3 : Z in nano-meter
    cd(PathName2);
    dotsDAT=[];
    for i=1:size(domainsDAT,1)
        dotstemp=dots{domainsDAT(i,6)};
        dotstemp(:,1:3)=dotstemp(:,1:3)/1000;  % change to um for fastmap.m
        dotsDAT=[dotsDAT;dotstemp];
    end

    dotsDAT(:,12)=0;

    if dimensionjudge==2
        outfilename=strcat(FileName2(1:end-4),'-',num2str(distancenm),'nm',...
                           '-Thrd',num2str(domainsize),'-Dom2Ddots.dat');
        save(outfilename, 'dotsDAT', '-ascii', '-tabs');
    elseif dimensionjudge==3
        outfilename=strcat(FileName2(1:end-4),'-',num2str(distancenm),'nm',...
                           '-Thrd',num2str(domainsize),'-Dom3Ddots.dat');
        save(outfilename, 'dotsDAT', '-ascii', '-tabs');
    end
    disp(outfilename);
end
%% fastmap the domained localizations/dots or clusters
disp('-------------------------------------------------------------------')
judge=input('Wanna Fastmap the Domained Localizations/Dots?   - ','s');
if judge=='y'
    ST08_1_fastmap()
end
%% display all the good domains by order of size
disp('-------------------------------------------------------------------')
judge=input('Wanna See Domains by the Order of Size? y or n   - ','s');
if judge=='y'
    close all
    
       scrsz = get(0,'ScreenSize');
       figure('OuterPosition',[1 0.5*scrsz(4) 0.75*scrsz(4) 0.5*scrsz(4)],'Color',...
          [0.313725501298904 0.313725501298904 0.313725501298904]);
        hdl1=gcf;    clf;
        figure('OuterPosition',[1 1 0.75*scrsz(4) 0.5*scrsz(4)]);
        hdl2=gcf;    clf;
    
    PlotCellNo   = ceil(sqrt(size(HistClusterNo,2)));
    PlotCellSize = 2 * max([max(HistXrange(2,:)) max(HistYrange(2,:))]);
    PlotZero=zeros([PlotCellNo^2,3]);
    
    k=1;
    PlotZero=zeros([PlotCellNo,3]);
    for i=1:PlotCellNo
        for j=1:PlotCellNo
            PlotZero(k,:)=[(j-1)*PlotCellSize,(PlotCellNo-i)*PlotCellSize,0];
            k=k+1;
        end
    end 
    [Y,SizeIndex]=sort(HistClusterNo(2,:));
    RoPlot=[];
    
    HistXrangeR=zeros([2,size(HistClusterNo,2)]);
    HistYrangeR=zeros([2,size(HistClusterNo,2)]);
    HistZrangeR=zeros([2,size(HistClusterNo,2)]);
    HistClusterNoR=zeros([2,size(HistClusterNo,2)]);
    HistDomainAreaR=zeros([2,size(HistClusterNo,2)]);   
    
    figure(hdl1)
    axes1 = axes('Parent',figure(hdl1),'ZColor',[1 1 1],...
        'YColor',[1 1 1],...
        'XColor',[1 1 1],...
        'Color',[0 0 0]);
    box(axes1,'on');    hold(axes1,'all');
    colormap('hsv')
    tubecount=0;
    for i=1:size(HistClusterNo,2)
        j=HistClusterNo(1,SizeIndex(i));
        RoPlotTemp=[domains{7,j}(:,1)+PlotZero(i,1),...
                    domains{7,j}(:,2)+PlotZero(i,2),...
                    domains{7,j}(:,3)+PlotZero(i,3)];
        Area=repmat(Y(i),size(domains{7,j},1),1);
                
        RoPlot=[RoPlot;[RoPlotTemp,Area]];
                
        HistXrangeR(1:2,i)     = [j;XrangeR(j)];
        HistYrangeR(1:2,i)     = [j;YrangeR(j)];
        HistZrangeR(1:2,i)     = [j;ZrangeR(j)];
        HistClusterNoR(1:2,i)  = [j;domains{2,j}];
        HistDomainAreaR(1:2,i) = [j;domains{4,j}];
        if domains{8,i}(1,2)/domains{8,i}(1,1)>1.5
            tubecount=tubecount+1;
        end
    end
     
    tubecount
    scatter3(RoPlot(:,1),RoPlot(:,2),RoPlot(:,3),...
             3,RoPlot(:,4),'filled');
    whitebg([0 0 0]);
    grid off
    axis equal
    view(0,90)
    zoom on
    title({[FileName(1:end-4)];...
           ['Total Clusters Numbers = ' num2str(length(images))];...
           ['Distance between Clusters in one Domain = ' num2str(distancenm) 'nm'];...
           ['Total Domains Numbers = ' num2str(domaincount)];...
           ['Threshold of Clusters per Domain = ' num2str(domainsize)];...
           ['Higher Ordered Domain No. = ' num2str(higherdomain)];...
           ['Colormap from Blue to Red means cluster number from few to many']},...
           'FontWeight','bold','FontSize',10,'FontName','Arial',...
           'Color',[1 1 1]); 
    colorbar('peer',axes1);
    
    if dimensionjudge==2
        fig3name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain2D-Sort.fig');
        fig3name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain2D-Sort.bmp');
        saveas(gcf, [fig3name1], 'fig');
        saveas(gcf, [fig3name2], 'bmp');
        
    elseif dimensionjudge==3
        fig3name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain3D-Sort.fig');
        fig3name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain3D-Sort.bmp');
        saveas(gcf, [fig3name1], 'fig');
        saveas(gcf, [fig3name2], 'bmp');
    end
    
    RoPlot(:,12)=0;

    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    nXrangeR=10*ceil((max(HistXrangeR(2,:))-min(HistXrangeR(2,:)))/std(HistXrangeR(2,:)));
    nYrangeR=10*ceil((max(HistYrangeR(2,:))-min(HistYrangeR(2,:)))/std(HistYrangeR(2,:)));
    nZrangeR=10*ceil((max(HistZrangeR(2,:))-min(HistZrangeR(2,:)))/std(HistZrangeR(2,:)));
    nclusternoR=10*ceil((max(HistClusterNoR(2,:))-min(HistClusterNoR(2,:)))/std(HistClusterNoR(2,:)));
    ndomainareaR=10*ceil((max(HistDomainAreaR(2,:))-min(HistDomainAreaR(2,:)))/std(HistDomainAreaR(2,:)));
    %%%%%%%%%%%%%%%%%%%%%  Revised by Q.Peter Su on Dec 28 2012    %%%%%%
    
    figure(hdl2)
      
      subplot(3,2,1);
      hist(HistXrangeR(2,:),nYrangeR);
      [HistXrangeR_n HistXrangeR_x]=hist(HistXrangeR(2,:),nXrangeR);
%       hold on
%       mu=[];sigma=[];x=[];
%       [mu,sigma]=normfit(HistXrangeR(2,:));
%       x=HistXrangeR(2,:);
%       ezplot(@(x)normpdf(x,mu,sigma))
%       hold off
      xlim([0 inf])
      ylim([0 inf])
      xlabel('XrangeR of Domains(um)','FontWeight','bold');
      ylabel('Occurrence');
      title({['X mean - ', num2str(mean(HistXrangeR(2,:)))];...
             ['X std - ',num2str(std(HistXrangeR(2,:)))]});
      grid on;zoom on;
      
%       subplot(3,2,1);hist(HistXrangeR(2,:),nXrangeR);
%       [HistXrangeR_n HistXrangeR_x]=hist(HistXrangeR(2,:),nXrangeR);
%       xlabel('XrangeR of Domains(um)','FontWeight','bold');
%       ylabel('Occurrence');
%       title({['X mean - ', num2str(mean(HistXrangeR(2,:)))];...
%              ['X std - ',num2str(std(HistXrangeR(2,:)))]});
%       grid on;zoom on;
      
      subplot(3,2,2);
      hist(HistYrangeR(2,:),nYrangeR);
      [HistYrangeR_n HistYrangeR_x]=hist(HistYrangeR(2,:),nYrangeR);
%       hold on
%       mu=[];sigma=[];x=[];
%       [mu,sigma]=normfit(HistYrangeR(2,:));
%       x=HistYrangeR(2,:);
%       ezplot(@(x)normpdf(x,mu,sigma))
%       hold off
      xlim([0 inf])
      ylim([0 inf])
      xlabel('YrangeR of Domains(um)','FontWeight','bold');
      ylabel('Occurrence');
      title({['Y mean - ', num2str(mean(HistYrangeR(2,:)))];...
             ['Y std - ',num2str(std(HistYrangeR(2,:)))]});
      grid on;zoom on;
      
          
          
      %%%%%%%%%%% Revised by Q.Peter SU on 07th July 2017 %%%%%%%%%%
      if max(HistZrange(2,:))==0
          subplot(3,2,3);
          hist(HistZrangeR(2,:));
          [HistZrangeR_n HistZrangeR_x]=hist(HistZrangeR(2,:));
%           hold on
%           mu=[];sigma=[];x=[];
%           [mu,sigma]=normfit(HistZrangeR(2,:));
%           x=HistZrangeR(2,:);
%           ezplot(@(x)normpdf(x,mu,sigma))
%           hold off
          xlim([-inf inf])
          ylim([-inf inf])
          xlabel('ZrangeR of Domains(um)','FontWeight','bold');
          ylabel('Occurrence');
          title({['Z mean - ', num2str(mean(HistZrangeR(2,:)))];...
              ['Z std - ',num2str(std(HistZrangeR(2,:)))]});
          grid on;zoom on;
      else
          subplot(3,2,3);
          hist(HistZrangeR(2,:),nZrangeR);
          [HistZrangeR_n HistZrangeR_x]=hist(HistZrangeR(2,:),nZrangeR);
%           hold on
%           mu=[];sigma=[];x=[];
%           [mu,sigma]=normfit(HistZrangeR(2,:));
%           x=HistZrangeR(2,:);
%           ezplot(@(x)normpdf(x,mu,sigma))
%           hold off
          xlim([0 inf])
          ylim([0 inf])
          xlabel('ZrangeR of Domains(um)','FontWeight','bold');
          ylabel('Occurrence');
          title({['Z mean - ', num2str(mean(HistZrangeR(2,:)))];...
              ['Z std - ',num2str(std(HistZrangeR(2,:)))]});
          grid on;zoom on;
      end
                  
%       subplot(3,2,3);hist(HistZrangeR(2,:),nZrangeR);
%       [HistZrangeR_n HistZrangeR_x]=hist(HistZrangeR(2,:),nZrangeR);
%       xlabel('ZrangeR of Domains(um)','FontWeight','bold');
%       ylabel('Occurrence');
%       title({['Z mean - ', num2str(mean(HistZrangeR(2,:)))];...
%              ['Z std - ',num2str(std(HistZrangeR(2,:)))]});
%       grid on;zoom on;
      


      subplot(3,2,5);
      hist(HistClusterNoR(2,:),nclusternoR);
      [HistClusterNoR_n HistClusterNoR_x]=hist(HistClusterNoR(2,:),nclusternoR);
%       hold on
%       mu=[];sigma=[];x=[];
%       [mu,sigma]=normfit(HistClusterNoR(2,:));
%       x=HistClusterNoR(2,:);
%       ezplot(@(x)normpdf(x,mu,sigma))
%       hold off
      xlim([0 inf])
      ylim([0 inf])
      xlabel('Clusters or Dots no. in Domains','FontWeight','bold');
      ylabel('Occurrence');
      title({['Cluster mean - ', num2str(mean(HistClusterNoR(2,:)))];...
          ['Cluster std - ',num2str(std(HistClusterNoR(2,:)))]});
      grid on;zoom on;
      
      
      subplot(3,2,6);
      hist(HistDomainAreaR(2,:),ndomainareaR);
      [HistDomainAreaR_n HistDomainAreaR_x]=hist(HistDomainAreaR(2,:),ndomainareaR);
%       hold on
%       mu=[];sigma=[];x=[];
%       [mu,sigma]=normfit(HistDomainAreaR(2,:));
%       x=HistDomainAreaR(2,:);
%       ezplot(@(x)normpdf(x,mu,sigma))
%       hold off
      xlim([0 inf])
      ylim([0 inf])
      xlabel('Domain Area in Sq. Pixel (25nm * 25nm)','FontWeight','bold');
      ylabel('Occurrence');
      title({['DomainArea mean - ', num2str(mean(HistDomainAreaR(2,:)))];...
             ['DomainArea std - ',num2str(std(HistDomainAreaR(2,:)))]});
      grid on;zoom on;
     
      whitebg([0 0 0]);
      
    if dimensionjudge==2
        outfilename=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain2D-Sort.dat');
        save(outfilename, 'RoPlot', '-ascii', '-tabs');
        outfilenameMAT=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain2D-Sort.mat');
        save (outfilenameMAT,'-mat');
    elseif dimensionjudge==3
        outfilename=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain3D-Sort.dat');
        save(outfilename, 'RoPlot', '-ascii', '-tabs');
        outfilenameMAT=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                       '-Thrd',num2str(domainsize),'-Domain3D-Sort.mat');
        save (outfilenameMAT,'-mat');
    end   
    
      
    if dimensionjudge==2
        fig4name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain2D-Sort-Hist.fig');
        fig4name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain2D-Sort-Hist.bmp');
        saveas(gcf, [fig4name1], 'fig');
        saveas(gcf, [fig4name2], 'bmp');
    elseif dimensionjudge==3
        fig4name1=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain3D-Sort-Hist.fig');
        fig4name2=strcat(FileName(1:end-4),'-',num2str(distancenm),'nm',...
                    '-Thrd',num2str(domainsize),'-Domain3D-Sort-Hist.bmp');
        saveas(gcf, [fig4name1], 'fig');
        saveas(gcf, [fig4name2], 'bmp');
    end   
    
    
end

%% See a sinlge domain with triangles surface 
disp('-------------------------------------------------------------------')
judge=input('Wanna See a Single Domain with Surface? y or n   - ','s');
while judge=='y'
    figure
    surfacenum=input('Type in the No. of the Ordered Domain in the Pic - ');
    j=HistClusterNo(1,SizeIndex(surfacenum));
    RoPlotTemp=[domains{7,j}(:,1),domains{7,j}(:,2),domains{7,j}(:,3)];
    try
    [triHull, vbOutside, vbInside] = AlphaHull(RoPlotTemp(:,1:3),distance);
    SurfArea=0;
    for m=1:size(triHull,1)
        plot=zeros([4,3]);
        k=1;
        for n=1:3
            plot(k,1:3)=RoPlotTemp(triHull(m,n),1:3);
            k=k+1;
        end
        plot(4,1:3)=RoPlotTemp(triHull(m,1),1:3);
        
        plot3(plot(:,1),plot(:,2),plot(:,3));
        hold on
        
        a=sqrt( (plot(2,1)-plot(1,1))^2 + ...
                (plot(2,2)-plot(1,2))^2 + ...
                (plot(2,3)-plot(1,3))^2 );
        b=sqrt( (plot(3,1)-plot(2,1))^2 + ...
                (plot(3,2)-plot(2,2))^2 + ...
                (plot(3,3)-plot(2,3))^2 );
        c=sqrt( (plot(4,1)-plot(3,1))^2 + ...
                (plot(4,2)-plot(3,2))^2 + ...
                (plot(4,3)-plot(3,3))^2 );
        p=(a+b+c)/2;
        s=sqrt(p*(p-a)*(p-b)*(p-c));
        SurfArea=SurfArea+s; 
    end
    hold on
    plot3(RoPlotTemp(:,1),RoPlotTemp(:,2),RoPlotTemp(:,3),'ro');
    axis equal
    whitebg([0 0 0]);
    title({[FileName(1:end-4)];...
           ['Total Clusters Numbers = ' num2str(length(images))];...
           ['Distance between Clusters in one Domain = ' num2str(distancenm) 'nm'];...
           ['Total Domains Numbers = ' num2str(domaincount)];...
           ['Threshold of Clusters per Domain = ' num2str(domainsize)];...
           ['Higher Ordered Domain No. = ' num2str(higherdomain)];...
           ['This is the No. -' num2str(surfacenum) '- Higher Ordered Domain'];...
           ['The Surface Area = ' num2str(SurfArea) 'square micron']},...
           'FontWeight','bold','FontSize',10,'FontName','Arial',...
           'Color',[1 1 1]); 
    catch ME
        disp('This Domain Cannot be Analysed by the Code, Soryy ! ');
    end
    judge=input('Wanna See Another Domain with a Surface?  y or n - ','s');       
end



























%function
end

