RCBF1

1. 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

2. 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

3. 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;

4. Calculate the left hand side of equation (6).

   
   rL = PET_int1 ./ PET_int2;

5. 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);

6. Calculate the right side of equation (6).

   
   rR = conv_int1 ./ conv_int2;

7. Generate the image through table lookup.

   
   k2 = lookup(rR, k2_lookup, rL);

8. Generate the image through table lookup and division.

   
   k2_conv_ints = lookup (k2_lookup, conv_int1, k2);
   K1 = PET_int1 ./ k2_conv_ints;

9. 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);