CORRECTBLOOD

1. Check the input arguments.

   
   if ((nargin < 5) | (nargin > 6))
      help correctblood
      error ('Incorrect number of input arguments.')
   end
   
   progress = 0;                   % defaults in case of no options vector -
   do_delay = 1;                   % may be overridden below
      
   if (nargin == 6)                % options vector
   
      if (length(options)>=1)      % options vector given; if it has an
         progress = options(1);    % element 1, that will be progress
      end
   
      if (length(options)>=2)      % delay correction toggle supplied
         do_delay = options(2);
      end
   
      if (length(options)>=3)      % value to use for delta supplied
         delta = options(3);
      else                         % no delta value given, so use 0
         delta = 0;                
      end
   
   end
   
   if (progress) 
      disp ('Showing progress');
   end
   
   if (~do_delay) 
      disp (['No delay-fitting will be performed; will use delta = ' ...
              int2str(delta)]);
   end

2. Do the initial setup.

   
   MidFTimes = FrameTimes + FrameLengths/2;
   first60 = find (FrameTimes < 60);       % all frames in first
   numframes = length(FrameTimes);         %  minute only
   
   tau = 4;                                % dispersion constant
   
   if (progress >= 2)
      figure;
      plot (ts_even, g_even, 'y:');
      title ('Blood activity: dotted=g(t), solid=g(t) + tau*dg/dt');
      drawnow
      hold on
   end

3. Do the dispersion correction.

   
   [smooth_g_even, deriv_g] = ...
        deriv (3, length(ts_even), g_even, (ts_even(2)-ts_even(1)));
   smooth_g_even(length(smooth_g_even)) = [];
   deriv_g(length(deriv_g)) = [];
   ts_even(length(smooth_g_even)) = [];
    
   g_even = smooth_g_even + tau*deriv_g;
   
   if (progress >= 2)
      plot (ts_even, g_even, 'r');
      drawnow
   end

4. Get ready to do delay correction.

   
   A = A (first60);                        % chop off stuff after 60 sec
   MidFTimes = MidFTimes (first60);        % first minute again
   
   if (progress >= 2)
      figure;
      plot (MidFTimes, A, 'or');
      hold on
      title ('Average activity across gray matter');
      old_fig = gcf;
      drawnow;
   end
   
   % Here are the initial values of alpha, beta, and gamma, in units of:
   %  alpha = (mL blood) / ((g tissue) * sec)
   %  beta = 1/sec
   %  gamma = (mL blood) / (g tissue)
   % Note that these differ numerically from Hiroto's suggested initial
   % values of [0.6, alpha/0.8, 0.03] only because of the different
   % units on alpha of (mL blood) / ((100 g tissue) * min).
   
   init = [.0001 .000125 .03];

5. Do the delay correction by performing several curve fittings.

   
   if (do_delay)
   
      if (progress), fprintf ('Performing fits...\n'), end
      deltas = -5:1:10;
      rss = zeros (length(deltas), 1);     % residual sum-of-squares
      params = zeros (length(deltas), 3);  % 3 parameters per fit
      options = [0 0.1];
   
      for i = 1:length(deltas)
         delta = deltas (i);
         if (progress), fprintf ('delta = %.1f', delta), end
   
         % Get the shifted activity function, g(t - delta), by shifting g(t)
         % to the right (ie. subtract delta from its actual times, ts_even)
         % and resample at the "correct" times ts_even).  Then do the 
         % three-parameter fit to optimise the function wrt. alpha, beta,
         % and gamma.  Plot this fit.  (Could get messy, but what the hell)
   
         shifted_g_even = lookup ((ts_even-delta), g_even, ts_even);
         g_select = find (~isnan (shifted_g_even));
   
         % Be really careful with the fitting.  If the algorithm you
         % choose makes args(2) a negative value, there will be
         % infinities in the result of b_curve, which will cause
         % the entire thing to bomb.
   
         options(2) = 1;
         options(3) = 1;
   
         final = fmins ('fit_b_curve', init, options, [], ...
                   shifted_g_even (g_select), ts_even (g_select), ...
                   A, FrameTimes, FrameLengths);
   
         params (i,:) = final;
         rss(i) = fit_b_curve (final, ...
                    shifted_g_even(g_select), ts_even(g_select), ...
                    A, FrameTimes, FrameLengths);
         init = final;
         if (progress)
            fprintf ('; final = [%g %g %g]; residual = %g\n', final, rss (i));
   
            if (progress >= 2)
               plot (MidFTimes, ...
                b_curve(final, shifted_g_even(g_select), ts_even(g_select), ...
                A, FrameTimes, FrameLengths));
               drawnow;
            end      % if graphical progress
         end      % if any progress
      end      % for delta
   
      [err, where] = min (rss);            % find smallest residual
      delta = deltas (where);              % delta for best fit
   
   end      % if do_delay

6. Now that we have a value for delta, perform the actual delay correction (shift the blood data). We must also remove NaN's from the new data. These appear if lookup cannot perform a lookup at a point (ie. there is no corresponding point in the table). They will therefore occur at the end points, where ts_even does not span (ts_even-delta).

   
   % At this point either we have performed the delay-correction fitting to 
   % get delta, or the caller set options(2) to zero so that delay-correction
   % was not explicitly done.  In this case, delta will have been set either
   % to zero or to options(3).
   
   Ca_even = lookup ((ts_even-delta), g_even, ts_even);
   
   nuke = find(isnan(Ca_even));
   Ca_even(nuke) = [];
   
   % Let's assume that the NaN's occur at the beginning or end of the data
   % (not in the middle), and that we can therefore modify ts_even without
   % screwing up the even time spacing.
   
   new_ts_even = ts_even;
   new_ts_even(nuke) = [];