Approximate Volumetric System Models for MicroSPECT

Approximate Volumetric System Models for MicroSPECT J. Gregora , S. Gleasonb , S. Kennelb , M. Paulusb , T. Bensona , J. Wallc University of Tennessee, Knoxville, TN 37996–3450 USA Oak Ridge National Laboratory, Oak Ridge, TN 37831-6010 USA c University of Tennessee Graduate School of Medicine, Knoxville, TN 37920 USA a

b

Abstract— A microSPECT system is being developed for the purpose of studying murine models of amyloidosis in vivo. The system is equipped with two detector heads, each consisting of a multi-anode photomultiplier tube coupled to a pixellated NaI(Tl) crystal array. Images are reconstructed using an OS-EM algorithm. In this paper, we describe the associated volumetric system models that we have developed in support of both pinhole and parallel hole collimation. These models, which are precomputed and stored to disk prior to reconstruction, are based on simple inner-product computations and a region-growing like search defined by the conic view that each detector pixel has of the image voxels. We provide illustrative experimental results based on phantom and mouse data.

I. I NTRODUCTION A microSPECT system is being developed for the purpose of studying the progression and regression of murine models of AA and AL-amyloidosis in vivo [1], [2]. The system is equipped with two detector heads, each consisting of a crossed-wire, multi-anode photomultiplier tube coupled to a pixellated NaI(Tl) crystal array [3]. Depending on the particular application, the detector heads can be individually configured to acquire either single pinhole or parallel hole collimated data. Images are reconstructed using the EM-ML algorithm [4] which we have implemented to support ordered subsets [5]. Let vectors x and y denote respectively the image and the projection data. Furthermore, let A be the system matrix that models the data acquisition process, and let vector c = AT 1 denote the column sums thereof. The iterative update equation for the EM-ML algorithm can then be expressed as: x(k+1) = (x(k) ⊗ AT (y ® Ax(k) )) ® c

(1)

where ⊗ and ® represent elementwise multiplication and division. An ordered subsets scheme results when the image update is based on using a sequence of subsets of the projection data as opposed to using all of it at once. Only the relevant portions of the system matrix are used at any given subiteration. Our implementation uses disjoint subsets that are formed by considering every ith element of the projection data and every ith row of the system matrix during the ith subiteration. The column sum vector is computed accordingly. In this paper, we focus on the design and implementation of the volumetric system models that we have developed in support of pinhole and parallel hole imaging. We provide illustrative experimental results based on a Jaszczak micro

phantom, a mouse bone scan, as well as transgenic and induced mouse models of amyloidosis. II. C ONIC D ETECTOR V IEW The purpose of the system model is to capture the physics of the data acquisition process. While effects such as attenuation and photon scattering could be included, we concentrate on encoding the sampling geometry of the system. Various approaches have been proposed. In some cases, the image voxels and the detector pixels are both represented by their geometric center points. This is common, for example, for line intersection models although sometimes subsampling is used to increase the model accuracy [6]. Collimation is taken into consideration by projecting the contour of each collimator aperture onto the detector and then restricting the calculations to the area thus defined. Another type of model uses solid angle based point spread functions. For that approach, the idea is to consider each image voxel a point source and then calculate the ratio of the area of each detector pixel illuminated thereby relative to the area of the corresponding emission sphere [7]. We refer to the literature for more detailed descriptions and for other alternatives. With reference to Fig. 1, our system models, which are precomputed and stored to disk prior to reconstruction, can be described as follows. Each detector pixel has a conic view of the image voxels. For pinhole collimation, the cone angle varies with the position of the detector pixel relative to the position of the pinhole. For parallel hole collimation, the cone angle is the same for all detector pixels. Both system models are based on the intersection volume of the cone with each image voxel. Rather than computing the exact intersection volume, we use an approximation. Let nc denote the normal vector corresponding to the conic view for a detector pixel being considered, and let nv denote a vector anchored at the apex of the cone and pointing toward the center of an image subvoxel. We then consider the subvoxel, and thus the voxel too, to be contained in the cone if the angle between the two vectors is less than or equal to θc , the cone angle. That is, nc · nv ≤ θc (2) acos |nc ||nv | Testing all image subvoxels to find the ones that are intersected by a cone is inefficient. We instead initialize the computation using a region-growing like search. The cone’s

0-7803-8701-5/04/$20.00 (C) 2004 IEEE

Detector column

Detector

Parallel hole collimators

nc

nc

θv

Pinhole aperture

θv Conic detector view

nv

nv

Conic detector view Sub−voxels

Fig. 1.

Sub−voxels

Conic view based volumetric system models for single pinhole imaging (left) and parallel hole imaging (right).

normal vector is used to find an image voxel that is guaranteed to be intersected. This voxel is then added to a queue. Subvoxel cone containment is tested for each voxel on the queue using Eq. (2). Voxels that fail the test are discarded. Voxels that pass the test, on the other hand, are added to the system matrix using the fraction of the voxel volume corresponding to the number of subvoxels intersected by the cone. Any of the twenty-six neighbors of the successful voxel that have not yet been processed are then added to the queue for subsequent consideration. When all detector pixels have been considered for all view angles, the computed intersection volumes are normalized to sum to one on a per voxel basis. Let the intersection volume for detector cone m with voxel n be written as vmn . Matrix element amn is then given by X amn = vmn / vmn . (3) m

When not using subsets, it is easily shown that this normalization ensures that the L1 -norms of vectors x(k+1) , x(k) , Ax(k) , and y are identical under the mild assumption that x(0) is non-negative. In other words, we achieve count preservation in that the total photon count associated with an image remains unchanged from iteration to iteration. The photon count associated with a forward projection of the image furthermore not only adds up to the same total, but equals that of the recorded projection data. III. C OMPUTATIONAL R ESULTS

The computational cost of the EM-ML algorithm is proportional to two times the number of nonzero elements in the system matrix both with respect to time, because each iteration requires a forward projection followed by a backprojection, and with respect to storage, because we need to store both the normalized volume intersection values as well as the corresponding voxel indices. The factors most likely to impact this number include the type of collimation, the size and number of projections, and the resolution of the reconstructed image. The computational work presented next is based on 60 projections, each 100 × 100 with a detector pitch of 1.5mm.

Data was acquired over the course of 30 minutes per data set. We made each image fill out the space defined by the orbit of the collimator for the given experiment but restricted the reconstruction to take place within the cylindrical field of view defined by the detector when it was smaller. The voxels were in all cases divided into 125 subvoxels. To reduce the time required to initialize the system matrix, we parallelized this part of the code using LAM-MPI based message passing. The reconstruction code itself was single threaded. The computations were performed using 16 commodity PCs, each equipped with 2.4GHz Pentium 4 CPUs and 2 Gbytes of memory. The computers, which run Linux, are interconnected via Gigabit ethernet. We first present results for a Jaszczak micro phantom filled with 99mTc. This data was acquired using a 0.5mm pinhole collimator. The cone angle was small, namely, 0.3 degrees. We reconstructed images for voxel resolutions of 0.6mm and 1.2mm. The former was of dimension 118 × 118 × 126 while the latter was 60 × 60 × 62. The two system matrices contained 150M and 57M non-zero elements and took respectively 9 and 4 minutes to initialize. During the queue based search for the voxels intersected by the detector cones, an average of 2000 voxels were processed per cone with a range of 1000 to 2600 for the 0.6mm image. The corresponding numbers were 900, 500, and 1200 for the 1.2mm image. These two sets of numbers should be compared with the 1.75M and 223K voxels that brute-force computations would have had to consider. The computational savings are thus substantial. One iteration of the EM-ML algorithm took 4 seconds for the 0.6mm image and 0.6 seconds for the 1.2mm image. Use of ordered subsets added a small but measurable amount of overhead. As an example, a 0.6mm reconstruction based on 8 ordered subsets took about 5 seconds per iteration. Likewise, a 1.2mm reconstruction increased to 2 seconds per iteration. Figure 2 shows a non-central slice from the two images, in both cases reconstructed using 16 iterations and 2 ordered subsets. We note that while the 1.6mm rods are difficult to see in the 1.2mm image, they are clearly seen in the 0.5mm image. The second set of computational results are for a 99mTc mouse bone scan. This data was acquired using 2mm parallel

0-7803-8701-5/04/$20.00 (C) 2004 IEEE

Fig. 2. Jaszczak micro phantom pinhole reconstructions: 0.6mm voxels (left) and 1.2mm voxels (right). Rod diameters are 1.2mm, 1.6mm, 2.4mm, 3.2mm, 4.0mm, and 4.8mm.

Fig. 3.

Bone scan parallel hole reconstruction.

hole collimation. The cone angle was relatively large, namely, 4.3 degrees. We made the voxel resolution equal to the detector pitch, i.e., 1.5mm, and reconstructed a 38 × 38 × 100 image. The system matrix contained 90M non-zero elements and took about 2 minutes to initialize. During the queue based search for the voxels intersected by the detector cones, an average of 900 voxels were processed per cone with a range of 150 to 1500. This is again substantially less work than that of a bruteforce computation which would have involved all 144K voxels. Each iteration of the EM-ML algorithm took about 2 seconds. Figure 3 shows a sagittal slice of the image produced after 8 iterations using 4 ordered subsets. Notice the brightness of the hind legs, the spine and the shoulders. IV. A MYLOID I MAGING R ESULTS We now present imaging results from two amyloid studies conducted using an early version of the microSPECT system. The single detector head of this prototype system was designed to produce 64×64 projections with a detector pitch of 1.25mm. Data was in both cases obtained for 60 views over the course of 30 minutes using 2mm parallel hole collimation. We here show 64 × 64 × 64 images that we reconstructed using 20 iterations of the EM-ML algorithm without ordered subsets. The first study involves a transgenic mouse model of systemic AA-amyloidosis [8]. This H2/huIL-6 mouse constitutively expresses the human interleukin-6 transgene rendering it in a state of chronic inflammation. Intravenous injection of

AEF, an amyloid enhancing factor consisting of an extract of spleen tissue from a mouse with amyloidosis, results in the rapid and systemic deposition of AA-amyloid in the spleen and liver and to a lesser degree the kidneys and the heart. Four mice were injected with 100µg AEF and for the purpose of control, four mice were administered an equivalent amount of saline solution. After seven weeks the mice received 600µCi of I-125 labeled SAP (a ubiquitous component of all types of amyloid), having two days prior been provided with 1% Lugol’s solution in their drinking water to block thyroid uptake of free radioiodine. The mice were sacrificed 24-hrs post SAP injection whereafter image data was acquired. Figure 4 shows an example of a reconstructed microSPECT image overlaid and co-registered with a CT image. The former represents a color map of the amyloid burden while the latter shows the anatomy of the diseased animal. The CT data was obtained using a highresolution MicroCAT II small animal x-ray CT system (ImTek, Inc., Knoxville, TN). These detailed 3D dual-modality views of systemic AA-amyloid disease in a mouse are, to the best of our knowledge, the first of their kind. We computed the biodistribution of the amyloid-bound I-125 SAP in these mice and compared that data with the equivalent I-125 SAP image counts. We obtained the former by measuring the specific activity associated with tissue samples gathered post mortem and the latter by first segmenting out the relevant organs in the CT images and then summing up the activity in the corresponding regions of the SPECT images. As indicated by Fig. 5, the two datasets are in excellent agreement. While a single scaling factor cannot be computed that will align all the I-125 SAP percent image counts with the corresponding percent injected dose per gram numbers, the same trends and organ activity ratios are observed. The second study deals with an induced mouse model of localized AL-amyloidosis. An amyloidoma composed of human AL-amyloid extract was introduced between the scapulae of a Balb/c mouse by administering a subcutaneous bolus of highly concentrated material. The amyloidoma is resolved naturally after two weeks, but the process can be accelerated dramatically following treatment with a fibril-reactive monoclonal antibody (mAb) designated 11-1F4 [9]. To examine the ability of the 11-1F4 antibody to bind AL amyloid in vivo, six mice received 50mg injections of human-derived AL amyloid. After a week, four mice were intravenously given an injection of I-125 labeled 11-1F4 mAb while two control mice received an I-125 labeled isotopematched control reagent. Following a three day period to permit the clearance of unbound 11-1F4 mAb, the mice were sacrificed and image data was acquired and processed as described above. Fig. 6 shows the accumulation of the radiolabeled 11-1F4 mAb at the site of the AL amyloidoma. Comparison of the image data with biodistribution results revealed excellent correlation. In both cases, significantly higher concentrations of the 111F4 mAb were observed in the amyloidoma compared with the heart, kidneys, spleen, and liver.

0-7803-8701-5/04/$20.00 (C) 2004 IEEE

Biodistribution Data 25

% ID/gm

20 15 10 5 0

Mouse1

Mouse2

Mouse3

Mouse4

Mouse3

Mouse4

I−125 SAP Data 60

% Image Counts

50 40 30 20 10 0

Mouse1

Mouse2

Fig. 5. Biodistribution data and I-125 SAP percent image counts. Black, gray, and white bars represent hepatic, spleenic, and pancreatic activity.

Fig. 4. MicroSPECT images of I-125 labeled SAP in a H2/huIL-6 mouse with AA-amyloidosis overlaid and co-registered with graylevel CT image (top) and skeleton rendered CT image (bottom). Bright object is splenic amyloid.

Fig. 6. MicroSPECT image of mouse induced with AL-amyloidosis overlaid on top of co-registered skeleton rendered CT image. Bright object to the right is the amyloidoma. Non-specific blood pool background to the left is heart and liver.

V. ACKNOWLEDGEMENTS

[3] A. Weisenberger, B. Kross, S. Majewski, V. Popov, M.F. Smith, B. Welch, R. Wojcik, S. Gleason, J. Goddard, M. Paulus, and S.R. Meikle, “Development and testing of a restraint free small animal SPECT imaging system with infrared based motion tracking,” IEEE Med. Imag. Conf., (Portland, OR), Oct. 2003. [4] L.A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging, vol. 1, pp. 113–122, 1982. [5] H. Hudson and R. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging, vol. 13, pp. 601– 609, 1994. [6] G. Zeng, G. Gullberg, B. Tsui, and J. Terry, “Three-dimensional iterative reconstruction algorithms with attenuation and geometric point response correction,” IEEE Trans. Nuc. Sci., vol. 2, pp. 693–702, 1991. [7] Z. Cho, J. Jones, and M. Singh, Foundations of Medical Imaging, Wiley, 1993. [8] A. Solomon, D.T. Weiss, M. Schell, R. Hrncic, C.L. Murphy, J. Wall, M.D. McGavin, H.J. Pan, G.W. Kabalka, and M.J. Paulus, “Transgenic mouse model of AA-amyloidosis,” Am. J. Pathol., vol. 154, pp. 1267–1272, 1999. [9] R. Hrncic, J. Wall, D. Wolfenbarger, C.L. Murphy, M. Schell, D. Weiss, and A. Solomon, “Antibody-mediated resolution of light chain-associated (AL) amyloid deposits,” Am. J. Pathol., vol. 157, pp. 1239–1246, 2000.

This work was supported by the NIH under grant number 1 R01 EB00789-01A2. The computer equipment was acquired as part of SInRG, a University of Tennessee grid infrastructure grant supported by the NSF under grant number EIA–9972889. Dr. T. Hoffman, US VA Biomolecular Imaging Center, Columbia, MO, provided the pinhole phantom data as well as the parallel hole bone scan data. R EFERENCES [1] J. Wall, M. Paulus, S. Kennel, S. Gleason, J. Baba, J. Gregor, D. Wolfenbarger, D. Weiss, and A. Solomon, “Imaging of primary (AL) amyloidosis with an amyloid reactive monoclonal antibody,” Xth Intl. Symp. on Amyloid and Amyloidosis, Tours, France, April 2004. [2] J. Wall, S. Kennel, M. Paulus, S. Gleason, J. Gregor, J. Baba, M. Schell, T. Richey, B. O’Nuallain, R. Donnell, P. Hawkins, D. Weiss, and A. Solomon, “Visualization and quantitation by dual modality highresolution microradiographic imaging of systemic amyloid deposits in a novel murine model of AA-amyloidosis,” Am. J. Pathology, submitted.

0-7803-8701-5/04/$20.00 (C) 2004 IEEE