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

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.

wolforth@pet.mni.mcgill.ca