- Check the input arguments. We want to make sure that we have
at least two input arguments (the filename and the slice number).
Then, for various numbers of input arguments, we want to set options
such as whether to display progress or perform blood delay and
dispersion correction. Finally, we check to ensure that the
argument passed as slice is a scalar (has a length of one).
if (nargin < 2)
help rcbf2
error ('Not enough input arguments');
elseif (nargin < 3)
progress = 0;
correction = 1;
elseif (nargin < 4)
correction = 1;
elseif (nargin > 4)
help rcbf2
error('Incorrect number of arguments.');
end
if (length(slice)~=1)
help rcbf2
error('<Slice> must be a scalar.');
end
- Get the images and image information. This includes all of the
frames for the slice in question, the frame start times, frame
lengths, mid-frame times, and the blood data. The blood data is
returned by resampleblood in a new time domain, resampled to
half second intervals. In order to keep our units consistent, we
convert the units of the PET images from
to 
. The units of the blood data will be
converted after cross-calibration correction is performed.
img = openimage(filename);
FrameTimes = getimageinfo (img, 'FrameTimes');
FrameLengths = getimageinfo (img, 'FrameLengths');
MidFTimes = FrameTimes + (FrameLengths / 2);
[g_even, ts_even] = resampleblood (img, 'even');
PET = getimages (img, slice, 1:length(FrameTimes));
PET = PET * 37 / 1.05; % convert to decay / (g_tissue * sec)
- Create the weighting functions, plus some integrals that are
used several times later in the program. PET_int1,
PET_int2, and PET_int3 are weighted integrals of the
PET data across frames (integrated with respect to time). Here we
also create a simple mask that excludes information outside of the
brain, and apply it to our integrals. This will ensure that areas
outside of the brain will be set to zero.
w2 = MidFTimes;
w3 = sqrt (MidFTimes);
ImLen = size(PET,1);
PET_int1 = C_trapz (MidFTimes, PET')';
PET_int2 = C_trapz (MidFTimes, PET' .* (w2 * ones(1,ImLen)))';
PET_int3 = C_trapz (MidFTimes, PET' .* (w3 * ones(1,ImLen)))';
% Now use PET_int1 to create a simple mask, and mask all three
% PET integrals. This does a good job of removing the
% outside-of-head data for CBF studies.
mask = PET_int1 > mean(PET_int1);
PET_int1 = PET_int1 .* mask;
PET_int2 = PET_int2 .* mask;
PET_int3 = PET_int3 .* mask;
- Correct the blood data for delay and dispersion, and apply the
cross-calibration factor. The correction will be accomplished by
performing a curve fitting on average activity across grey matter.
Therefore, getmask allows the user to select grey matter by
choosing a threshold value. The average grey matter activity is then
passed to correctblood. The units of the blood data are
converted from
to 
by multiplying by the cross-calibration factor.
Therefore, to maintain consistent units, we convert to 
.
XCAL = 0.11;
g_even = g_even*XCAL*37/1.05; % units are decay / (g_tissue * sec)
if (correction)
mask = getmask (PET_int1);
A = (mean (PET (find(mask),:)))';
[ts_even, Ca_even, delta] = correctblood ...
(A, FrameTimes, FrameLengths, g_even, ts_even, progress);
else
Ca_even = g_even;
end
- Generate three very useful lookup tables. These are
evaluated for 
, 
, and 
as the weights
(
). These lookup tables are used in all subsequent evaluations of
these three weighted integrals. The values of 
used in the
lookup are also chosen here.
k2_lookup = (-10:0.05:10) / 60;
[conv_int1,conv_int2,conv_int3] = findintconvo (Ca_even,ts_even,...
k2_lookup, MidFTimes, FrameLengths, 1, w2, w3);
- Generate some additional useful integrals. These are used more
than once in the subsequent code, and so are calculated in advance to
try and increase the speed slightly. Ca_mft is the blood
data averaged over each frame. This is taken as the value of the
blood data at the mid-frame time.
Ca_mft = nframeint (ts_even, Ca_even, FrameTimes, FrameLengths);
Ca_int1 = C_trapz(MidFTimes, Ca_mft);
Ca_int2 = C_trapz(MidFTimes, (w2 .* Ca_mft));
Ca_int3 = C_trapz(MidFTimes, (w3 .* Ca_mft));
- Generate the lookup table relating to values of the right
hand side of equation (10). We also calculate the left
hand side of equation (10), which will be used in the
generation of a image. Since we will be looking up values of
from values of rR (the right hand side of equation
(10)), we sort the table on rR.
rL = ((Ca_int3 .* PET_int1) - (Ca_int1 .* PET_int3)) ./ ...
((Ca_int3 .* PET_int2) - (Ca_int2 .* PET_int3));
rR = ((Ca_int3 * conv_int1) - (Ca_int1 * conv_int3)) ./ ...
((Ca_int3 * conv_int2) - (Ca_int2 * conv_int3));
% Now, we must have the k2/rR lookup table in order by rR; however,
% we also want to keep k2_lookup in the original order. This
% is because the first lookup uses rL as a lookup into rR to
% find k2 (which requires that rR be monotonic, ie. sorted) whereas
% subsequent lookups all use k2 to find conv_int{1,2,3} -- which
% requires that k2_lookup be monotonic. So k2_lookup will be the
% list of k2's in order, and k2_sorted will be the same list, but
% in order according to the sorted rR.
[rR,sort_order] = sort (rR);
k2_sorted = k2_lookup (sort_order);
- Generate the image, through a simple table lookup.
Values of are chosen by finding the value of where
rL and rR are equal.
k2 = lookup(rR, k2_sorted, rL);
- Generate the
image by evaluating the numerator of
equation (10). All of the time consuming calculations
have already been performed, and we can evaluate the 
terms
through table lookup.
K1_numer = ((Ca_int3*PET_int1) - (Ca_int1 * PET_int3));
K1_denom = (Ca_int3 * lookup(k2_lookup,conv_int1,k2)) - ...
(Ca_int1 * lookup(k2_lookup,conv_int3,k2));
K1 = K1_numer ./ K1_denom;
- Generate the
image by evaluating equation
(7) directly. Once again, we may get values for
the complicated part of the equation through simple table lookup.
V0 = (PET_int1 - (K1 .* lookup(k2_lookup,conv_int1,k2))) / Ca_int1;
- Clean up the images by removing NaN's and infinities (setting
them to zero).
nuke = find (isnan (K1));
K1 (nuke) = zeros (size (nuke));
nuke = find (isinf (K1));
K1 (nuke) = zeros (size (nuke));
nuke = find (isnan (V0));
V0 (nuke) = zeros (size (nuke));
nuke = find (isinf (V0));
V0 (nuke) = zeros (size (nuke));
- Finally, close the image file so that everything gets cleaned
up nicely.
closeimage (img);
