loadstack.m

function stack = loadstack(filename, varargin)
%
% Load one of Peter's 2P stacks into memory. Note that each 'stack'
% actually contains two stacks -- if there are N frames in the TIFF
% then 1:N/2 is the green channel data and N/2:end is the red
% channel, which for us contains the ttl signal from pype.
%
% at least the left edge of the image is messed up, so a 1px region
% around the image is replaced with the mean..
%
% stack structure contains:
%   .g - green channel data (normalized 0-1 base on a 16bit input)
%   .ttl - pype's TTL gating pulse signal over time
%   .stim_onsets - start times (frame number) for each stimulus period
%   .stim_offsets - stop times (frame number) for each stimulus period
%

opts.lscan = 1;
opts.justinfo = 0;
opts.interleaved = 0;
opts.force = 0;
opts = pargin(opts, varargin);

stack.info = imfinfo(filename);
nframes = size(stack.info, 1);
if opts.justinfo
  fprintf('%d frames (%d green/%d red)\n', ...
          nframes, nframes/2, nframes/2);
  stack = [];
  return
end

mfilename = strrep(filename, '.tif', '.mat');
if ~opts.force && exist(mfilename, 'file')
  stack = load(mfilename);
  stack = stack.s;
  fprintf('loaded preprocessed .mat version\n');
  return
end

% loadtiff is WAY faster than imread..
s = double(loadtiff(filename));


% find frames where ttl line is high for the majority of the frame

if opts.lscan
  % stack generated by labview - 1st half green, 2nd red:
  %
  % split stack into green and red channels; red channel contains
  % TTL signal, so we're going to convert it to time domain and
  % then toss the image data below. There's some garbage around
  % the perimeter of the images, so we'll replace with the mean
  % to avoid edge problems as well.
  stack.g = s(:,:,1:(nframes/2));
  x = stack.g ./ bitshift(1,16);
  x(1,:,:) = NaN; x(end,:,:) = NaN;
  x(:,1,:) = NaN; x(:,end,:) = NaN;
  x(isnan(x)) = nanmean(unravel(x));
  stack.g(isnan(x)) = NaN;

  % pull normalized red channel to look for TTL signals
  r = s(:,:,(1+(nframes/2)):end);
  r = r ./ (max(r(:)) - min(r(:)));
  stack.ttl = [];
  for n = 1:size(r,3)
    if sum(unravel(r(:,:,n) > 0.50)) / ...
          (size(r,1) .* size(r,2)) > 0.90
      stack.ttl(n) = 1;
    else
      stack.ttl(n) = 0;
    end
  end
  if stack.ttl(1) == 1
    % this is not really right -- makes sync correct, but
    % first stimulus should really be discarded because onset
    % time is not really known.
    warning('ttl starting high, forcing low');
    stack.ttl(1:10) = 0
  end
  stack.stim_onsets = find(diff(stack.ttl) == 1);
  stack.stim_offsets = find(diff(stack.ttl) == -1);
  if length(stack.stim_onsets) > length(stack.stim_offsets)
    stack.stim_onsets = stack.stim_onsets(1:end-1);
  end
else
  if opts.interleaved
    stack.g = s(:,:,1:2:end);
  else
    stack.g = s;
  end
  stack.stim_onsets = [];
  stack.stim_offsets = [];
end
stack.filename = filename;


%% these should not be hardcoded!
stack.fps = 3.91;                       % frames per sec (Hz)
stack.xres = 1.0;                       % horizontal um/pixel
stack.yres = 1.0;                       % vertical um/pixel