function [Abl, intValues] = EK_getFRAPvals(tStack)

% Auto detect FRAP area and intensity values as a function of time. 
% tStacks is a movie stack, size(Xpixels, Ypixels, nframes). 



% THIS PART OBSOLETE. 
% relBleach int is a correction factor to rescale the image taken at the
% end of bleach period. Since intensity is typically 20-100 times higher
% than read-out intensity, the mean pix values are also higher either by
% the same amount (if bleach depth is not 100%) or by some other amount (if
% bleach depth =100%). 
% If bleach depth is <100%, then corrected value = mean pix val/relBleachInt.  
% If bleach depth is 100% (i.e. the rescaled pix values < noise), then the
% corrected value is = noise. 
% example: noise is 104. Frame (usually 6th frame) acquired just before
% bleach ended gives 248, and the relBleachInt=20. If the bleach depth
% were<100%, then the corrected value would be 248/20=12.4, but since that
% is < noise, bleach depth must be 100%. So the corrected value is = noise
% = 104. 
%
% Do in future: extract relBleachInt from metadata. 
% For now, need to lookup in the original files (need to know laser int). 
%


%tStack replaces series1
%series1 = bfopenOutput{1, 1}; % series are different channels, or z-position
                        % for a time lapse greyscale image, series1 is a
                        % nf-by-2 cell array where nf is the number of
                        % frames. 
                        
%fr1 = series1{1,1}; % get frame 1. This is double, but scaled as if int16
fr1 = tStack(:,:,1); %frame 1
fr1 = uint16(fr1); % convert to int16 (what it should really be)
fr1d = im2double(fr1); % convert back to double, but now with proper scaling ([0 1])
fr1a = imadjust(fr1d, [min(fr1d(:)) max(fr1d(:))], []); %stretch limits to  [0,1]

[level, EM] = graythresh(fr1a); % global image threshold 
                               % using Otsu's algorithm 
                               % for this to work well, need adjusted
                               % image to full range. 
msg=sprintf('level and effectiveness metric used for threshold are %g %g', level, EM);
disp(msg);          
BW = im2bw(fr1a, level); 
%imtool(fr1)
%imtool(BW)

se = strel('disk',6); %radius 6 pix
BWc = imclose(BW,se);
%se2 = strel('disk', 6);
BWe = imerode(BWc,se);

%% visualize goodness of detection: 
% C = imfuse(BWe, fr1); 
% imtool(C);
% % BWoutline = bwperim(BWe);
% % Segout = fr1;
% % Segout(BWoutline) = 255;
% % figure, imshow(Segout), title('outlined original image');
% % imtool(BWe)

%% bleached area, in pixels: 
s = regionprops(BWe, fr1, 'Area'); % get bleached area, pix
Abl = cat(1, s.Area); % s is structure, extract into array

% check only 1 value in area
if size(Abl,1)>1
    error('more than 1 region detected for bleached area!!');
end

nf= size(tStack, 3); % number frames
intValues = zeros(nf,1);

for k=1:nf
    fr_k = uint16( tStack(:,:,k) ); % k'th frame
    FRAPstats = regionprops(BWe, fr_k, 'MeanIntensity'); % do measurements on 
                                             % original, unadjusted image 
    intValues(k)=FRAPstats.MeanIntensity;                                         
end
 
% THIS PART BELOW IS OBSOLETE. IT ONLY APPLIES TO THE NIKON TIRFM WHICH HAD
% A BUG
% %correct for the frame acquired just before the end of the bleach period: 
% % find that frame, if any: 
% 
% if max(intValues)>intValues(1)
%     bfrind = find( intValues == max(intValues) ); %find fr acq. durign bleach
%     if intValues(bfrind)/relBleachInt < Fnoise
%         intValues(bfrind) = Fnoise; 
%     else
%         intValues(bfrind)=intValues(bfrind)./relBleachInt;
%     end
% else
%     %nothing to correct
% end





