%% Wrote by Q.Peter Su and Yujie Sun on Jan 08 2013
%  Part  I : Generate the BF Frame Index for ImageJ
%  Part II : Generate Driftave Matrix 

%% Input
%  1. value*.txt output by ImageJ-Menu-Image-Stacks-Plot Z-axis profile
%  2. IntThreshold, pick up all the bright field out of fluo frame
%  3. Crop*.tif output by ImageJ using the BF index
%  4. IntThreshold 2, pick up all the bright BF out of all BF frame
%  5. polynomial factor to fit the drift curve.

%% Output
%  1. FileName(1:end-4),'-BfIndex.txt'
%  2. FileName(1:end-4),'-FluoIndex.txt'
%  3. FileName(1:end-4),-new.tif
%  4. FileName(1:end-4),Poly#-drift.fig
%  5. FileName(1:end-4),Poly#-drift.mat

%% Revised by Q.Peter Su on Jan 11 2013
%  IntThreshold-1 ~1000 is used to pick out all the bright field (BF) images
%              out of original *.tif file
%  IntThreshold-2 ~3500 is used to pick out all the Bright BF images
%              out of all the BF frames

%% Revised by Q.Peter Su on Jan 15 2013
%  xi = interp1(t,fx,ti, 'spline');
%  yi = interp1(t,fy,ti, 'spline');
%  Avoid abnormal polynomial fitting.

%% Revised by Q.Peter Su on Jan 18 2013
%  fit the curve with interp1 : YI = INTERP1(X,Y,XI,METHOD)
%  ti = 1:frames;
%  xi = interp1(t,x_off,ti, 'cubic');
%  yi = interp1(t,y_off,ti, 'cubic');

%% Revised by Q.Peter Su on Jan 25 2013
%  ti = min(corr_offset(:,7)):max(corr_offset(:,7));

%% Revised by Q.Peter Su on Mar 13 2013
%  Add intensity and frame range input to avoid some bright fluo frames.
%  IndexBF=find(meanintensity(:,2)>IntThreshold(1) & ...
%               meanintensity(:,2)<IntThreshold(2) & ...
%               meanintensity(:,1)>FrameThreshold(1) & ...
%               meanintensity(:,1)<FrameThreshold(2));

%% Revised by Q.Peter Su on Apr 16 2013
%   1. 
%   thresholdbf   =input('Type in the IntThreshold of Bright BF image.[Low   High] - ');
%   thresholdframe=input('Type in the Frame Range of Bright BF image. [Early Late] - ');
%                   disp('Check the Generated *BF.tif File in ImageJ, Find Out ...   ');
%                   disp('Other Frame No. You DONOT Wanna Analyze, as a Row Vector.  ')
%   indexbftemp   =input('Here is an Example : [a1, a2, a3:a10, a22, a34:a44, a56] - ');
%
%   2.
%   to avoid the repeat frames of Bright Field which is caused by the
%   abnormal condition of shutter
%   Delete continuous 20 frame of Bright field

%% Revised by Q.Peter Su on Apr 28 2013
%axes1 = axes('Parent',figure1,'YDir','reverse',...

%% Revised by Q.Peter Su on Oct 22 2013
%indexbfContinue2=[];%find(indexbfContinue1(:,2)>=50);

%% Revised by Q.Peter Su on Oct 22 2013
%  dlmwrite(new1,indexbfout,'precision',6);
%  dlmwrite(new2,indexfluoout,'precision',6);

%% Revised by Q.Peter Su on Dec 05 2013
%  use several windows to average the BF frames to avoid the noise
%
%  window=input('Type in a Window for Average BF Frames (like 50 or 100) - ');
%  windowNo=floor(size(meanint,1)/window);
%
%  StackNew=zeros([size(Stack,1) size(Stack,2) windowNo]);
%  for i=1:windowNo
%      Itemp=Stack(:,:,((i-1)*window+1):(i*window));
%      I=median(Itemp,3);

%% Revised by Q.Peter Su on Dec 05 2013
%   avoid mistake between frame to frame, every step
%         imwrite(uint16(round(I)),new,'Compression','none','WriteMode','append');
%         output0 = dftregistration(fft2(StackNew(:,:,i)),fft2(StackNew(:,:,1)),precision);
%         outputy = output0(3);
%         outputx = output0(4);
%         output  = [output;[output0 outputy outputx]];

%% Revised by Q.Peter Su on Dec 05 2013
%   avoid missing drift info of very first and very last frame 
%        ti = 1:size(meanintensity,1);
%        %ti = min(corr_offset(:,7)):max(corr_offset(:,7));
%        xi = interp1(t,x_off,ti, 'cubic');
%        yi = interp1(t,y_off,ti, 'cubic'); 
        
%% revised by Q.PEter Su on Mar 29 2015 
%   to avoid some frames without BF, which are at the begining or end %%%%%
%   xi(1:t(1))=0;
%   yi(1:t(1))=0;
%   xi(t(end):end)=x_off(end);
%   yi(t(end):end)=y_off(end);
        
        

%%
function ST22_IndexBF_and_GenDrift_V2_2(Xpixel,Ypixel)
close all
%% Part I : Generate the BF Frame Index for ImageJ
%%%%%%%%%%%%%%  Import molecule-list*.txt file output by ImageJ  %%%%%%%%%
[FileName,PathName] = uigetfile('*value*.txt',...
    'Value of Mean *.txt File Output by ImageJ',...
    'MultiSelect', 'off');
cd(PathName);disp(PathName);disp(FileName);
Afname=char(FileName);
disp('Import Data ... ing ... ');
meanintensity = importdata(Afname);

%%%%%%%%%%%%%%%%%%%%% Preview the mean intensity  %%%%%%%%%%%%%%%%%%%%%%
close all
    plot(meanintensity(:,1),meanintensity(:,2),'MarkerFaceColor',[1 0 0],...
     'MarkerEdgeColor',[1 0 0],'Marker','o');
    xlabel('Frame No.');
    ylabel('Mean Intensity');
    title('Preview the Intensity of all Frames');
    whitebg([0 0 0]);
    zoom on
    
%%%%%%%%%%%%%%%%%%%% Input IntThreshold %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('---')
disp('To Separate Frames You Donnot Wanna Show in the Final Image, ')
disp('Including the Bright Field Iamges and the High Background Images.')

IntThreshold=[];
while ~(length(IntThreshold)==2)
    IntThreshold  =input('Type in the IntThreshold of the Intensity. [Low High] - ');
    if length(IntThreshold)==2
        disp(num2str(IntThreshold));
    else
        disp('Input Format is Wrong !!! ');
    end
end

FrameThreshold=[];
while ~(length(FrameThreshold)==2)
    FrameThreshold=input('Type in the FrameThreshold of the Frame. [Early Late] - ');
    if length(FrameThreshold)==2
        disp(num2str(FrameThreshold));
    else
        disp('Input Format is Wrong !!! ');
    end
end

IndexBF=find(meanintensity(:,2)>IntThreshold(1) & meanintensity(:,2)<IntThreshold(2)...
            & meanintensity(:,1)>=FrameThreshold(1) & meanintensity(:,1)<=FrameThreshold(2));
            % Revised by Q.Peter Su on Mar 13 2013
            
indexbfout=IndexBF(:)';  % output should be a row vector
IndexFluo=meanintensity;
IndexFluo(indexbfout,:)=[];
IndexFluo(:,2)=[];
indexfluoout=IndexFluo(:)';  % output should be a row vector
% save index file
new1=strcat(FileName(1:end-4),'-Thrld[',num2str(IntThreshold(1)),'-',num2str(IntThreshold(2)),']-BfIndex.txt');
new2=strcat(FileName(1:end-4),'-Thrld[',num2str(IntThreshold(1)),'-',num2str(IntThreshold(2)),']-FluoIndex.txt');
% Revised by Q.Peter Su on Oct 22 2013
dlmwrite(new1,indexbfout,'precision',6); % '6'means the precision of number
dlmwrite(new2,indexfluoout,'precision',6); % avoid '1.111 E6'
% Copy the BF index to clipboard
indexcopy=importdata(new1,'s');
clipboard('copy', indexcopy{1});
%% Part II : Generate Driftave Matrix 
disp('-- ');
disp('In ImageJ, Click the -Image- in the Menu List, ');
disp('Then -Stacks-, -Tools-, -Make Substack...-,');
disp('Paste the String of BF Index Already in Clipboard,');
disp('Donot Select the -Delete slices from original stack-');
judge=input('Finished? Press Enter Key to Continue - ','s');

% Select the 'Bright Field Only' tif stack
[FileName,PathName] = uigetfile({'*bf*.tif','Multilayer BF TIFF-Files (*.tif)'},...
                                 'Select the BrightField-Only Stack'); 
                                 %open dialog for *.stk files 
cd(PathName);
% Import the BF *.tif file
Stack=uint16([]);
disp('Import Data ...ing...');
if FileName~=0
    f=sprintf('%s%s',PathName,FileName); 
    c=1;
    while c~=0
        try 
            Stack(:,:,c)=imread(f,c);
            if mod(c,1000)==0
                disp(['Read Frame of -' num2str(c) '-']);
            end
        catch
            break
        end
        c=c+1;
    end
end
disp(['Read Frame of -' num2str(c-1) '-']);

meanint=zeros([size(Stack,3),1]);
%globalmax=10; globalmin=16000;
for i=1:size(Stack,3)
    I=double(Stack(:,:,i));
    meanint(i,1)=mean(mean(I));
    %if min(min(I))< globalmin
    %    globalmin=min(min(I));
    %end
    %if max(max(I))> globalmax
    %    globalmax=max(max(I));
    %end
end
meanint=[IndexBF meanint];

close all
plot(meanint(:,1),meanint(:,2),'MarkerFaceColor',[1 0 0],...
     'MarkerEdgeColor',[1 0 0],'Marker','o');
zoom on
xlabel('Frame No. in the Raw Movie');
ylabel('Mean Intensity of the BF Images');
title('Preview the IntThreshold of Bright BF image');
% input bf IntThreshold for bright BF
disp('------------------------------------------------------------------');

thresholdbf=[];
while ~(length(thresholdbf)==2)
    thresholdbf   =input('Type in the IntThreshold of Bright BF image.[Low   High] - ');
    if length(thresholdbf)==2
        disp(num2str(thresholdbf));
    else
        disp('Input Format is Wrong !!! ');
    end
end

thresholdframe=[];
while ~(length(thresholdframe)==2)
    thresholdframe=input('Type in the Frame Range of Bright BF image. [Early Late] - ');
    if length(thresholdframe)==2
        disp(num2str(thresholdframe));
    else
        disp('Input Format is Wrong !!! ');
    end
end
                disp('Check the Generated *BF.tif File in ImageJ, Find Out ...   ');
                disp('Other Frame No. You DONOT Wanna Analyze, as a Row Vector.  ')
indexbftemp   =input('Here is an Example : [a1, a2, a3:a10, a22, a34:a44, a56] - ');

% Revised by Q.Peter Su on Apr 16 2013
%  to avoid the repeat frames of Bright Field which is caused by the
%  abnormal condition of shutter

indexbfContinue1=diff(meanint(:,1));
indexbfContinue1(:,2)=0;
n=0;
for i=1:size(indexbfContinue1,1)
    if indexbfContinue1(i)==1
        n=n+1;
        indexbfContinue1(i,2)=n;
    else
        n=0;
        indexbfContinue1(i,2)=n;
    end
end
indexbfContinue2=[];%find(indexbfContinue1(:,2)>=50);
indexbfContinue=[];
for j=1:size(indexbfContinue2,1)
    range=[indexbfContinue2(j)-indexbfContinue1(indexbfContinue2(j),2) :...
           indexbfContinue2(j)];
    indexbfContinue=[indexbfContinue;range'];
end
%
indexbf2=[indexbfContinue;...
          indexbftemp';...
          find(meanint(:,2)>thresholdbf(2));...
          find(meanint(:,2)<thresholdbf(1));...
          find(meanint(:,1)>thresholdframe(2));...
          find(meanint(:,1)<thresholdframe(1))];

Stack(:,:,indexbf2)=[];
meanint(indexbf2,:)=[];  
disp('-----------------------------------------------------------------');
strout=strcat('Left Bright BF images Frame Number = ',' -',num2str(size(Stack,3)),'- Frames');
disp(strout);
disp('-----------------------------------------------------------------');

%%
%globalmax=10; globalmin=16000;
%for i=1:size(Stack,3)
%    I=double(Stack(:,:,i));
%    %meanint(i,1)=mean(mean(I));
%    if min(min(I))< globalmin
%        globalmin=min(min(I));
%    end
%    if max(max(I))> globalmax
%        globalmax=max(max(I));
%    end
%end
%%
window=[];
while isempty(window)
    window=input('Type in a Window for Average BF Frames (like 50 or 100) - ');
end
windowNo=floor(size(meanint,1)/window);
indexwindow=1:window:size(meanint,1);
indexwindow=indexwindow(1:windowNo);
new=sprintf('%s%s-new-avg%s.tif',PathName,FileName(1:end-4),num2str(window));
% Precision Input
%if nargin<1
    precision=100;
%end
output=[];
StackNew=zeros([size(Stack,1) size(Stack,2) windowNo]);
for i=1:windowNo
    Itemp=Stack(:,:,((i-1)*window+1):(i*window));
    I=median(Itemp,3);
    StackNew(:,:,i)=I;
    if i==1  %for the first frame
        imwrite(uint16(round(I)),new,'Compression','none');  
        strout=strcat('BF Frame No. -',num2str(i),'/',num2str(windowNo),'-');
        disp(strout);
        output=[0 0 0 0 0 0]; 
        %[error,diffphase,net_row(y)_shift,net_col(x)_shift,yshiftfromorigin,xshiftfromorigin ]
    else
        imwrite(uint16(round(I)),new,'Compression','none','WriteMode','append');
%       output0 = dftregistration(fft2(StackNew(:,:,i)),fft2(StackNew(:,:,i-1)),precision);
        output0 = dftregistration(fft2(StackNew(:,:,i)),fft2(StackNew(:,:,1)),precision);
%       outputy = output(i-1,5)+output0(3);
        outputy = output0(3);
%       outputx = output(i-1,6)+output0(4);
        outputx = output0(4);
        output  = [output;[output0 outputy outputx]];
        if mod(i,10)==0
        strout=strcat('BF Frame No. -',num2str(i),'/',num2str(windowNo),'-');
        disp(strout);
        end
    end
end

corr_offset=[output meanint(indexwindow,1)];

if nargin<1
    Xpixel=input('Please input X pixel size (in nm)   - ');
    Ypixel=input('Please input Y pixel size (in nm)   - ');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%frames=size(indexfluoout,2);
%%% frame no. change to fluo-only frames 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
x_off = Xpixel/1000*medfilt1(corr_offset(:,6),5);        % X in um
y_off = Ypixel/1000*medfilt1(corr_offset(:,5),5);        % Y in um
%x_off = corr_offset(:,6) * Xpixel/1000;        % X in um
%y_off = corr_offset(:,5) * Ypixel/1000;        % Y in um
t     = corr_offset(:,7);

%% interp fitting drift curve
close all

  ti = 1:size(meanintensity,1);
  %ti = min(corr_offset(:,7)):max(corr_offset(:,7));
  xi = interp1(t,x_off,ti, 'cubic');
  yi = interp1(t,y_off,ti, 'cubic'); 
  
  %%%%%%% revised on Mar 29 2015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%% to avoid some frames without BF, which are at the begining or end %%%%%
  xi(1:t(1))=0;
  yi(1:t(1))=0;
  xi(t(end):end)=x_off(end);
  yi(t(end):end)=y_off(end);
  %%%%%%% revised on Mar 29 2015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  
figure1 = figure('Color',...
    [0.313725501298904 0.313725501298904 0.313725501298904]); 
  axes1 = axes('Parent',figure1,'YDir','reverse',...
         'ZColor',[1 1 1],'YColor',[1 1 1],...
         'XColor',[1 1 1],...
         'Color',[0 0 0]);
  hold(axes1,'all');  
  plot(xi,yi,...
      'Parent',axes1,'MarkerFaceColor',[1 0 0],...
      'MarkerEdgeColor',[1 0 0],...
      'Marker','.',...
      'LineWidth',3,...
      'Color',[1 0 0]);
  title({[FileName(1:end-4)];...
         ['Drift Curve as Shown Below -- '];...
         ['Blue to Red Means Early to Late']},...
         'FontWeight','bold','FontSize',10,'FontName','Arial',...
         'color',[1,1,1]);
  xlabel('X in um');  
  ylabel('Y in um');
  %zlabel('Z in um');  
  hold on
  scatter(x_off,y_off,20,t,'filled');
  colorbar('peer',axes1);
  axis equal
  box on
  grid on
  hold off
  
%% generate driftave
    Driftave = zeros(size(ti,2), 5);
    Driftave(:,1) = ti;
    Driftave(:,2) = xi;   %x drift in um
    Driftave(:,3) = yi;   %y drift in um
    %Driftave(:,4) = zi;  %
    Driftave(:,4) = 0;    %z drift in um
%% auto save file
saveas(gcf, [FileName(1:end-4) '-Drift.fig'], 'fig');
saveas(gcf, [FileName(1:end-4) '-Drift.bmp'], 'bmp');
save([FileName(1:end-4) '-Drift.mat'],'Driftave');

%saveas(gcf, [FileName(1:end-4) '-Poly' num2str(polynomial) '-drift.fig'], 'fig');
%save([FileName(1:end-4) '-Poly' num2str(polynomial) '-drift.mat'],'Driftave');


%%
disp(FileName);
judge_dualcolor=input('-g- for Green, -b- for Blue, then Generate 2-Color Drift Info. - ','s');
if judge_dualcolor=='g'
    ST04_2_gen2Cdrift2_Red2Green(Xpixel,Ypixel,judge_dualcolor)
elseif judge_dualcolor=='b'
    ST04_2_gen2Cdrift2_Red2Green(Xpixel,Ypixel,judge_dualcolor)
end


end