- Check the input arguments. If the number of arguments is
incorrect, then print out the help information for the function,
followed by an explanatory error message. Additionally, we want to
ensure that the arguments themselves make sense. Therefore, we check
that slice is a scalar, by checking that its largest dimension
has a size of one (ie. it is a 1x1 matrix). If it is not, then we
print the help, followed by an error message explaining the problem.
if (nargin == 2)
progress = 0;
elseif (nargin ~= 3)
help rcbf1
error('Incorrect number of arguments.');
end
if (length(slice)~=1)
help rcbf1
error('<Slice> must be a scalar.');
end
- Get the images and image information. This includes all of the
frames for the slice being analyzed, the frame start times, the frame
lengths, the mid-frame times, and the blood data. The frame time
information is returned by simple enquiries with the data handle
img returned by openimage. The blood data is returned
by resampleblood, which gives the blood data in an evenly
sampled time domain (every half second). The blood data is then
cross-calibration corrected by multiplying by the factor XCAL.
In order to maintain consistent units throughout the analysis, this
data must then have its units converted from
to
. The actual PET images are then
retrieved, and the units are again converted to
. Any negative values present in the images are set to
zero.
img = openimage(filename);
FrameTimes = getimageinfo (img, 'FrameTimes');
FrameLengths = getimageinfo (img, 'FrameLengths');
MidFTimes = FrameTimes + (FrameLengths / 2);
[g_even, ts_even] = resampleblood (img, 'even');
% Apply the cross-calibration factor.
XCAL = 0.11;
g_even = g_even*XCAL*37/1.05; % units are decay / (g_tissue * sec)
Ca_even = g_even; % no delay/dispersion correction!!!
PET = getimages (img, slice, 1:length(FrameTimes));
PET = PET * 37 / 1.05; % convert to decay / (g_tissue * sec)
PET = PET .* (PET > 0); % set all negative values to zero
ImLen = size (PET, 1); % num of rows = length of image
- Calculate some useful integrals that are used several times in
the rest of the program. PET_int1, and PET_int2 are
weighted integrals of the PET data across frames (integrated with
respect to time). To get clean images out of the analysis, we wish to
set all points outside of the head to zero. Therefore, we create a
simple mask, and apply it to the weighted integrals.
PET_int1 = trapz (MidFTimes, PET')';
PET_int2 = trapz (MidFTimes, PET' .* (MidFTimes * ones(1,ImLen)))';
mask = PET_int1 > mean (PET_int1);
PET_int1 = PET_int1 .* mask;
PET_int2 = PET_int2 .* mask;
- Calculate the left hand side of equation (6).
rL = PET_int1 ./ PET_int2;
- Assign the
values to be used in the lookup table, and then
generate some useful weighted integrals.
k2_lookup = (0:0.02:3) / 60;
[conv_int1, conv_int2] = findintconvo (Ca_even, ts_even, k2_lookup,...
MidFTimes, FrameLengths, 1, MidFTimes);
- Calculate the right side of equation (6).
rR = conv_int1 ./ conv_int2;
- Generate the
image through table lookup.
k2 = lookup(rR, k2_lookup, rL);
- Generate the
image through table lookup and division.
k2_conv_ints = lookup (k2_lookup, conv_int1, k2);
K1 = PET_int1 ./ k2_conv_ints;
- Clean up the
image by setting any NaN's and infinities to
zero, and close the image data set.
nuke = find (isnan (K1) | isinf (K1));
K1 (nuke) = zeros (size (nuke));
closeimage (img);