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