Next: FINDINTCONVO Up: Annotated Program Listings Previous: RCBF2

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) = [];


wolforth@pet.mni.mcgill.ca