Integrating Images

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 $\times$ 128 image. (Of course, only about $\frac{1}{3}$ 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:

$$\int_{t_{1}}^{t_{n}} f(t) dt \approx \sum_{i=1}^{n} f(t_{i}) (\Delta t)_{i},$$ (1)

where f(t) is one time-activity curve and the $\Delta t$'s are the frame lengths. Since we have f(t) 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 $\Delta t$ 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 x-points, and a vector or matrix of y points corresponding to each x-point. If the y's are given as a matrix, then each row corresponds to one xpoint. 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 x-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².