% function to create an averaged image from STORM movies 
% as a representation of conventional image

% copyright @Rahul Mar 08, 2012
%modified from AvImage to calculate average images per N number of frames 

function AvImages=ST03_AvImage2()
clc
clear all;
close all;
[FileName,PathName] = uigetfile(... 
    {'*.tif',  'All Data Files (*.tif)';...
    '*RSI*.tif',' All Sheet stacks(*RSI*.tif)'; ...
    %'*RSI*a.tif','Sheet stack a (*a.tif)'; ...
    %'*RSI*b.tif','Sheet stack b (*b.tif)'; ...
    %'*Epi*a.tif','Epi stack a (RSI*a.tif)'; ...
    %'*Epi*b.tif','Epi stack b (Epi*b.tif)'; ...
    '*Epi*.tif',' All Epi stacks(*Epi*.tif)'}, ...
    'Select files to use for average image generation',...
    'MultiSelect', 'on');

disp('Which set of images to analyze?');
setofimages = input('all for all(default).  ', 's');
cd(PathName);
numberfiles=length(FileName);
disp('How many frames to pool for averaging? ');
pooledframes = input('Use integer multiple of movie length.  ');
movielength = input('How many frames per movie?  ');
inc=1;n=1;
count =0;
subim = ceil(movielength/pooledframes);
AvImages = cell(1,numberfiles*subim); 

while (n<=numberfiles)
     filename = char(FileName(n));   
if (filename(end-3:end)=='.tif')  
    disp('\n');
    disp(filename); %revised by QPeter Su on oct 18
    strgcount = ['*** File number  ' num2str(n) '  of  ' num2str(numberfiles) '  files ***'];
    disp(strgcount);
    
   switch lower(setofimages)
       case 'ab'
        TT = (filename(end-4)=='a') | (filename(end-4)=='b');
        if (TT)
            count = count +1;
            start = 1+(count-1)*subim;
            tend = count*subim;
            AvImages(1,start:tend) = ZaveImage2(filename, subim);
        end
     
    case 'ori'
      
      TT = (filename(end-4)~='a') & (filename(end-4)~='b');
      if (TT)
         count = count +1;
          start = 1+(count-1)*subim;
            tend = count*subim;
         AvImages(1,start:tend) = ZaveImage2(filename, subim);
      end
    
    case 'all'     
        count = count +1;
        start = 1+(count-1)*subim;
        tend = count*subim;
        AvImages(1,start:tend) = ZaveImage2(filename, subim);
     
    otherwise
        display('No files selected for analysis!');
     
   end     
end
n=n+inc;
end 
%% select ROIs with beads
Bd = input('Do you want to save ROIs with beads? press y or n', 's');

if Bd == 'y'
    disp('Merged RGB image shows start(Red), middle(Green) and end(Blue)');
    a(:,:,1) = AvImages{1,1}/max(AvImages{1,1}(:));
    a(:,:,2) = AvImages{1,ceil(numberfiles*subim/2)}/max(AvImages{1,ceil(numberfiles*subim/2)}(:));
    a(:,:,3) = AvImages{1,numberfiles*subim}/max(AvImages{1,ceil(numberfiles*subim)}(:));
    roisel = [];
    imshow(a);
    title('Select regions to crop and save bead movies');
    
    while Bd == 'y'
        h = imrect;
        position = wait(h);
        pos = getPosition(h);
        roisel = [roisel pos'];
        Bd = input('Do you want to select another bead? press y or n', 's');
    end
end
    [dim noofbeads] = size(roisel);
    roisel = uint16(roisel);
%% save a tiff movie of averaged file
outputfilename =input('Enter output file name for full frame movie?', 's');
outfilename  = strcat(outputfilename, '.tif'); 
bds = [1:1:noofbeads];
bdno = num2str(bds');
beadfiles = [repmat(outputfilename,noofbeads,1) bdno repmat('b.tif',noofbeads,1)];

if exist(outfilename, 'file')
       delete(outfilename);
end

for cc = 1:length(AvImages)
   for bno = 1:noofbeads
       ROI = uint16(AvImages{1,cc}(roisel(2,bno):(roisel(2,bno)+roisel(4,bno)),roisel(1,bno):(roisel(1,bno)+roisel(3,bno))));
       imwrite(ROI, beadfiles(bno,:),'tif', 'Compression', 'none', 'WriteMode', 'append'); 
   end 
end
end

