We can view the matrix PET as a collection of 21 PET images
for the same brain slice but progressing in time, or as a collection
of 16,384 time-activity curves (TAC)-one for every pixel in the 128
128 image. (Of course, only about
of the
image data is actually data from inside the head. This matter will
be addressed when we consider image masking below.)
For purposes of integrating images, we will consider PET as a collection of TACs. Thus, we are simply performing 16384 numerical integrations. This may seem a daunting task, but MATLAB's vectorized approach to mathematics makes it surprisingly easy. In fact, let us calculate a rectangular integration using the formula:
where is one time-activity curve and the
's are the
frame lengths. Since we have
stored as a row
vector-remember, each column of PET is the entire image for
one frame, so each row of PET corresponds to a single pixel
across time-and
stored as a column vector, the
multiply/sum operation of equation (1) is simply a
dot product. And since a matrix multiplication is nothing more than
a sequence of dot products, all 16,384 TACs can be quite easily
integrated with a single MATLAB command:
PETsummed1 = PET * FrameLengths;
However, rectangular integration is a very crude approach to numerical integration. Trapezoidal integration provides a slightly more refined method, and is almost as easy to implement. Furthermore, it is provided by the MATLAB function trapz-which, if you care to examine it (try type trapz in MATLAB), is essentially a one-line operation as well. (Incidentally, one of the functions provided by EMMA is ntrapz. ntrapz currently has exactly the same restrictions and capabilities as trapz, but is considerably faster due to being written in C using MATLAB's C interface. For the purposes of this discussion, the two functions are equivalent.)
trapz takes two arguments: a vector of -points, and a vector
or matrix of
points corresponding to each
-point. If the
's
are given as a matrix, then each row corresponds to one
point. Note that this is different from the EMMA standard, where each
column of PET corresponds to one time point. Thus, we
will have to transpose PET using the MATLAB ' operator
before passing it to trapz. Also, trapz will return the
16,384 integrals as a row vector, so we must transpose the results to
get them back into the form we expect. Finally, trapz expects
points on the
-axis rather than widths of intervals. Thus, we will
pass MidFTimes to it, and not FrameLengths. The full
integration is performed by the single statement
PETsummed2 = trapz (MidFTimes, PET')';
Note: This method of calculating the integral is
extremely wasteful of memory due to transposing the PET images. For a
full discussion of MATLAB memory issues, please consult the online
document Controlling MATLAB Memory
Use
.