iceTROPY_process.m

function [ flagImAnis, avAnisotropies, avIntensities ] = ...
    iceTROPY_process( fullFileDat, fullFileGf )
% iceTROPY_process evaluates anisotropy image of input data
%   Uses Registration, G-factor-calibration, and signal data
%   And returns some visualised outputs. 

% 1. Determine inputs
%    Image files and filtering

% Read in the raw image data and g-factor calibration
% Use fitsread() if the file is .fits, otherwise use imread
if( strcmp(fullFileGf(end-3:end),'fits') )
  imGcal = fitsread(fullFileGf);
else
  imGcal = imread(fullFileGf);
end
[~, nameDat, ~] = fileparts(fullFileDat); % Use for (console) output

% Read whether vertical flips and background subtraction are needed
flagFlipudGf    = evalin('base','flagFlipudGf');
flagFlipudDat   = evalin('base','flagFlipudDat');
flagSubBuffer   = evalin('base','flagSubBuffer');
flagVideo       = evalin('base','flagVideo');
flagShowVideo   = evalin('base','flagShowVideo');
flagVidFramewise= evalin('base','flagVidFramewise');
filterNumber    = evalin('base','filterNumber');
cameraAreaSetup = evalin('base','cameraAreaSetup');

% Read in regions of interest from iceTROPY_defineROI
% Here, the three initial ~,~,~ are for areaBG,areaPAR,areaPPD
% (which are presently not needed in this function)
% areaANISsmall is the region in which average anisotropy will be
% calculated. 
[ ~,~,~,areaANIS,areaANISsmall ] = iceTROPY_defineROI( cameraAreaSetup );

% Read in the image registration transform 
% (generated by "iceTROPY_registration.m" or loaded into base workspace)
mytform = evalin('base','mytform');

% Define smoothing filter (TO DO: MODIFY TO ALTERNATIVE CHOICES IF NEEDED
% Try a Gaussian filter (warning: this is digitised and nominally clipped)
gridXX = meshgrid(-3:3);
gridYY = gridXX';
myFsigSq = 0.6^2;
myFilter = exp(-(gridXX.^2 + gridYY.^2)/(2*myFsigSq) );


% 2. Process anisotropy data
%    (a) Determine G-factor image
%    (b) Determine r image
% (a) G-factor image 
[imGcalBG,imGcalPAR,imGcalPPD]=iceTROPY_segmentation(imGcal,flagFlipudGf);

BGgcal    = mean( imGcalBG(:) ); % Evaluate background
imGcalPAR = imGcalPAR - BGgcal; % Subtract background
imGcalPPD = imGcalPPD - BGgcal;
imGcalPAR = double(imGcalPAR);  % Make double, from (say) 16 or 32 bit
imGcalPPD = double(imGcalPPD);
% Align PAR onto PPD signal
imGcalPARreg = imtransform(imGcalPAR, mytform, 'Xdata', ...
                [1 size(imGcalPPD,2)], 'Ydata', [1 size(imGcalPPD,1)] );
% Filter the data (if required)
if(filterNumber == 1)
  imGcalPARregFil = conv2(imGcalPARreg, myFilter, 'same');
  imGcalPPDfil    = conv2(imGcalPPD, myFilter, 'same');
else
  imGcalPARregFil = imGcalPARreg; % Option 2 is no filter
  imGcalPPDfil    = imGcalPPD;
end
% G-factor Map
imGcalPARregFilCrop = imcrop(imGcalPARregFil, areaANIS);
imGcalPPDfilCrop    = imcrop(imGcalPPDfil, areaANIS);
G = imGcalPARregFilCrop ./ imGcalPPDfilCrop;


% 2 (b) Anisotropy (r) image 
%     (i) Determine whether video analysis is needed
%    (ii) Do anisotropy analysis once, or in a for loop for video
if(flagVideo ==0)
  numberOfFrames = 1;   % Process one frame 
else                    % Process movie data input
  myImInfo = imfinfo(fullFileDat,'tif');  % Assume Video is a TIF stack
  numberOfFrames = numel(myImInfo);
  if(flagVidFramewise ==0) % If we want to store the whole video data:
    imageAnisotropy = zeros( areaANIS(4)+1, areaANIS(3)+1, numberOfFrames);
    % This may run out of memory in some cases
    % Although it may be OK up to 100+ MB of input 
  else                     % Write output framewise into a folder
   
  end
end
avAnisotropies = zeros(numberOfFrames,1); % Pre-allocate memory for results
avIntensities  = zeros(numberOfFrames,1);

% Set the waitbar to update not more than 100 times...
dividerWaitbar = 10^(floor(log10(numberOfFrames))-1);
progress = waitbar(0,'Working on image stack');
% Now evaluate anisotropy:
for lpF = 1:numberOfFrames
     
  % Read in a frame from the image data file:
  if(flagVideo ==0)                             % If data is one frame,
   if( strcmp(fullFileDat(end-3:end),'fits') ) % and the file is .fits 
     imDat = fitsread(fullFileDat);
   else
     imDat = imread(fullFileDat);  % Use imread if not .fits (probably tif)
   end 
  else  % For video data assume the data is a tif stack
   imDat = imread(fullFileDat,'tif', lpF, 'Info',myImInfo);
  end

% TO DO: make BG subtraction work if wanted (needs another input)!
  if flagSubBuffer         % FROM DATA ONLY? G-factor should be brighter
   imBuffer = imread([myDirData,myFileBuffer]); % BG subtraction
   imDat = imDat - imBuffer;
  end 

[imDatBG,imDatPAR,imDatPPD]=iceTROPY_segmentation(imDat,flagFlipudDat);

BGdat    = mean( imDatBG(:) );     % Again, assumes flat background. 

imDatPAR = imDatPAR - BGdat;
imDatPPD = imDatPPD - BGdat;
imDatPAR = double(imDatPAR);
imDatPPD = double(imDatPPD);
% Data - Registration
imDatPARreg = imtransform(imDatPAR, mytform, 'Xdata', ...
                [1 size(imDatPPD,2)], 'Ydata', [1 size(imDatPPD,1)] );
 if(filterNumber == 1) % Filter the data if required
  imDatPARregFil = conv2(imDatPARreg, myFilter, 'same');
  imDatPPDfil    = conv2(imDatPPD, myFilter, 'same');
 else 
  imDatPARregFil = imDatPARreg; % Option 2 is no filter
  imDatPPDfil    = imDatPPD;
 end
  
imDatPARregFilCrop = imcrop(imDatPARregFil, areaANIS);
imDatPPDfilCrop    = imcrop(imDatPPDfil, areaANIS);
% Now evaluate anisotropy
% USING "r" for this frame of imageAnisotropy
r = (imDatPARregFilCrop - G.*imDatPPDfilCrop) ./ ...
    (imDatPARregFilCrop + 2*G.*imDatPPDfilCrop);

imageIntensity = (imDatPARregFilCrop + 2*G.*imDatPPDfilCrop);

imageAnisSmall = imcrop(r, areaANISsmall);
imageIntsSmall =  imcrop(imDatPARreg, areaANISsmall)+ ...
                  2 * imcrop(imDatPPD, areaANISsmall);
avAnisotropies(lpF) = mean(imageAnisSmall(:)); % Should be OK
avIntensities(lpF)  = mean(imageIntsSmall(:));
  
  if(flagVidFramewise)  % If we are not storing the image (stack) output
    imageAnisotropy(:,:,lpF) = avAnisotropies(lpF); % Not an image
    rT = r; % Truncate data to anisotropy range 0-1, for saving as double
    rT(rT<0) = 0;
    rT(rT>1) = 1;
    imwrite(rT,['vid/frame', int2str(lpF),'.tif'],'tif'); % Save tif
    
    iI = imDatPARregFilCrop + G.*imDatPPDfilCrop; % image Intensity
    imwrite(uint16(iI),['int/frame', int2str(lpF),'.png'],'png'); % Save tif
  else                  % Else store the output in MATLAB memory
    imageAnisotropy(:,:,lpF) = r;
  end
    
  if (round(lpF/dividerWaitbar)==lpF/dividerWaitbar) 
   waitbar(lpF/numberOfFrames,progress); 
  end
end % End of single frame analysis (or one instance of the analysis loop)
close(progress);

% 3. Visualise outputs
%  3.1 If data is a single frame: 
%    (a) Show G-factor image
%    (b) Show r image
%    (c) Plot anisotropy histogram
%  3.2 If data is a video:
%    (a) Show g-factor image
%    (b) Show r-image evolution
%    (c) Plot mean anisotropy and intensity over time

if(flagVideo==0)
figure(1)
 imagesc(G);
 title('G-factor','FontSize',14)
 colorbar
 axis image % Equal axis scale both direction, and tight fit around data

figure(2)
 imagesc(imageAnisotropy);
 colorbar
 title('Anisotropy image','FontSize',14)
 axis image
 caxis([0 0.4])

% figure(3)
% hist(imageAnisotropy(:), [0:0.01:0.4]);
% title('Anisotropy histogram','FontSize',14)
% xlim([0 0.4])

imageAnisSmall = imcrop(imageAnisotropy, areaANISsmall);
figure(3)
 hist(imageAnisSmall(:), 0:0.01:0.4 );
 title('Anisotropy histogram (crop)','FontSize',12)
 xlim([0 0.4]) % May need to rescale for 2-photon measurements
 ylabel('Pixels')
 xlabel('Anisotropy')
 set(gca,'FontSize',12);

figure(4)
 imagesc(imageIntensity);
 title('Total intensity','FontSize',14)
 colormap(gray)
 colorbar
 axis image % Equal axis scale both direction, and tight fit around data

anisMean = mean(imageAnisSmall(:));
disp(nameDat);
disp(anisMean);
anisStd = std(imageAnisSmall(:));
end

% 3.2 Video output:
if(flagVideo==1)
figure(1)
 imagesc(G);
 title('G-factor','FontSize',14)
 colorbar
 axis image

  if(flagShowVideo)
  figure(2)
   for lpF = 1:numberOfFrames
    imagesc( imageAnisotropy(:,:,lpF) );
    colorbar
    title('Anisotropy image','FontSize',14)
    caxis([0 0.4])  
    axis image
    drawnow; % Note: this can be very slow. 
   end
  end
figure(3)
  %   plot(1:numberOfFrames,avAnisotropies)
  %   title('Anisotropy versus time','FontSize',14) 
  %   xlabel('Frame number','FontSize',14) 
  %   ylabel('Average anisotropy','FontSize',14) 
  [AX,~,~]=plotyy(1:numberOfFrames,avIntensities, ...
                    1:numberOfFrames,avAnisotropies);
   % '~,~' replaces 'H1,H2' handles which are not used here
  title('Anisotropy versus time','FontSize',14) 
  legend('Intensity','Anisotropy')
  xlabel('Frame number','FontSize',14)
  set(get(AX(1),'Ylabel'),'String','Intensity','FontSize',14) 
  set(get(AX(2),'Ylabel'),'String','Anisotropy','FontSize',14) 
  
  
end

assignin('base','imageAnisotropy',imageAnisotropy);
assignin('base','imageIntensity',imageIntensity);
assignin('base','G',G);
assignin('base','avAnisotropies',avAnisotropies);
assignin('base','avIntensities',avIntensities);
assignin('base','numberOfFrames',numberOfFrames);
flagImAnis = 1;

iceTROPY_analyse_region;
end