progress in mesh based spatio-temporal reconstruction

0 downloads 0 Views 479KB Size Report
tomographic reconstruction. In addition, it can also provide a unified framework for tomographic reconstruction of organs. (e.g., the heart) that undergo non-rigid ...
PROGRESS IN MESH BASED SPATIO-TEMPORAL RECONSTRUCTION Jovan G. Brankov, Ricard Delgado Gonzalo, Yongyi Yang, Mingwu Jin, and Miles N. Wernick Department of Electrical and Computer Engineering Illinois Institute of Technology, Chicago, IL 60616, USA

ABSTRACT In this paper, we present a new methodology for calculation of a 2D projection operator for emission tomography using a content-adaptive mesh model (CAMM). A CAMM is an efficient image representation based on adaptive sampling and linear interpolation, wherein non-uniform image samples are placed most densely in regions having fine detail. We have studied CAMM in recent years and shown that a CAMM is an efficient tool for data representation and tomographic reconstruction. In addition, it can also provide a unified framework for tomographic reconstruction of organs (e.g., the heart) that undergo non-rigid deformation. In this work we develop a projection operator model suitable for a CAMM representation such that it accounts for several major degradation factors in data acquisition, namely object attenuation and depth-dependent blur in detector-collimator response. The projection operator is calculated using a raytracing algorithm. We tested the developed projection operator by using Monte Carlo simulation for single photon emission tomography (SPECT). The methodology presented here can also be extended to transmission tomography.

I.

INTRODUCTION

In tomography it is customary to represent volumetric data using uniform sampling and voxel basis functions. In [1] we first proposed an alternative approach, in which the image to be reconstructed is represented by non-uniform sampling and mesh model basis functions (see Figure 1). Such an image representation is called a content-adaptive mesh model (CAMM). This model involves sampling of the image domain on a non-uniform grid (mesh nodes), followed by partitioning into a collection of non-overlapping patches (mesh elements). The image function is then determined over each mesh element (e.g. triangles in 2D and tetrahedrons in 3D) by interpolation from the image samples (nodal values). With a CAMM representation, image samples are placed most densely in regions having fine detail. In recent years, we have studied and developed efficient procedures for mesh generation. We developed a mesh generation procedure in [1] and [2] for 2D images, which was later extended for 3D image in [3] and for vector-valued (e.g., color) images in [4]. We first presented the use of a CAMM (and mesh in general) for emission tomography in [5]-[7]. There we showed that when compared to several commonly used pixel-based methods for image reconstruction, the CAMM reconstruction approach achieved the best performance in terms of defect detection measured by channelized Hotelling observer (CHO), signal to noise ratio (SNR), and computation time. Subsequently, in [8] and [9] we developed and demonstrated the promise of use of a deformable CAMM for gated cardiac single photon emission tomography (SPECT) imaging. In a deformable CAMM, the mesh structure is allowed to deform over time, making it well suited for representing organs exhibiting nonrigid motion (see Figure 2). In our previous development, we did not develop a methodology for direct calculation of a projection operator specific for a mesh representation. Instead, we resorted to a two-step approach whereby the projection operator is first calculated for a pixel representation, then combined with an interpolation operator that connects a mesh model to the pixel representation. In this paper, we develop a methodology for a direct 2D projection operator calculation in a mesh model. This is an extension of our recent work in [10]. In particular, this projection model allows for incorporation of both non-uniform attenuation due to the object absorption and the depth-dependent blur in the detector-collimator spatial response in SPECT. It is worth noting that the methodology presented here can be easily extended to transmission tomography or positron emission tomography (PET).

Research supported by NHLBI grant HL065425 Computational Imaging VI, edited by Charles A. Bouman, Eric L. Miller, Ilya Pollak, Proc. of SPIE-IS&T Electronic Imaging, SPIE Vol. 6814, 68140V, © 2008 SPIE-IS&T · 0277-786X/08/$18

SPIE-IS&T/ Vol. 6814 68140V-1 2008 SPIE Digital Library -- Subscriber Archive Copy

Figure 1. Reconstructed emission image (left), attenuation image (center), and the mesh structure (right) for representing both images.

60 Motion trajectories

55

Frame 16

50 45

Frame 1

Y

40 35 30 25 20 15 60 50 40 30 20 X

2

4

6

8

10

12

14

16

Frame No.

Figure 2. Deformable mesh is allowed to deform from frame to frame. To facilitate the projection operator calculation which incorporates non-uniform attenuation, a single CAMM should be used such that it accurately represents both the emission image and the attenuation image. For this we make use of our mesh generation procedure developed for vector-valued images [4]. This was first demonstrated for multi-modality image tomographic reconstruction in [11], in which a single mesh structure was generated for co-registered MRI and SPECT images where the anatomical structure obtained from an MRI image served as a spatial prior for SPECT image reconstruction. In the literature a great many methods have been developed for improving the quality of reconstructed images in tomography (see [12] or [13] for a review). Most of these methods are pixel-based, i.e., the image is represented and computed directly in a pixel basis. Some methods based on object modeling have also been described in the literature (e.g.,[14]-[18]). These methods typically assume a priori a geometric model of the object being imaged. For example, generalized cylinder models were used in [14] for image reconstruction from incomplete projections; parametric surface models were investigated in [15][16]. In addition, spherically-symmetric basis functions were used in [19][20] or point cloud representation in [21]. These methods have a similar philosophy to our proposed approach in that a model is used to combat the ill-posed nature of the reconstruction problem. In [22] a hybrid-grid approach was proposed where fine-grid pixels were used in the heart region and coarse-grid pixels were used in the rest of the cardiac image, aiming to reduce

SPIE-IS&T/ Vol. 6814 68140V-2

computation yet preserving accuracy where it matters. While it has a similar spirit to our proposed approach, the use of a CAMM in our approach allows for truly content-adaptive and systematic placement of the image samples.

II. MESH REPRESENTATION OF EMISSION AND ATTENUATION MAP Let f (p) denote an image function over a domain D ⊂ R 2 where p = ( x , y ) ∈ D . Mesh modeling starts with dividing this domain into a total of T non-overlapping mesh elements Tl , where l = 1, 2,..., T . The image function f (p) is represented over each element as an affine function, that is, a first degree polynomial of two variables. To do so each mesh element is defined by set of tree nodal positions { p1l = ( x1l , y1l ), p l2 = ( x2l , y 2l ), p l3 = ( x3l , y3l ) } and a vector consisting T

of three nodal values fl = ⎡⎣ f (p1l ), f (p l2 ), f (p l3 )⎤⎦ . The approximate values of the function f (p) can be now calculated as: T fˆ ( x, y ) = [ Al , Bl , Cl ][ x , y ,1] , ( x, y ) ∈ Tl

where Al , Bl and Cl are the interpolation coefficients which can be computed from nodal values and nodal positions as shown below: l

[ Al , Bl , Cl ]T

⎛ x1 ⎜ = ⎜ x2l ⎜ xl ⎝ 3

−1

y1l 1⎞ ⎟ y 2l 1⎟ fl = al fl y

l 3

(1)

1⎟⎠

where we define al as interpolation matrix over Tl , which we will use later. As pointed out earlier, we use a common mesh structure to represent both the attenuation map and emission image, and thus: T µˆ ( x, y ) = ⎡⎣ Alµ , Blµ , Clµ ⎤⎦ [ x, y ,1] , ( x, y ) ∈ Tl (2) where Alµ , Blµ and Clµ are the interpolation coefficients for the attenuation map. These coefficients can be computed as above. Note that the interpolation matrix al is the same for both the emission image and attenuation map. The example in Fig. 1 shows the CAMM structure generated by using the procedure described in [1] and [4] for the emission image fˆ ( x , y ) and absorption map µˆ( x , y ) . In our experiments, the emission map was obtained from a smoothed version of the filtered backprojection (FBP) reconstruction.

III. MESH TOMOGRAPHY MODEL In [3], we showed that the mesh tomographic model can be expressed as follows: E [g] = Hfm

(3)

where fm denotes a vector formed from the nodal values of the mesh structure, E [⋅] is the expectation operator and g denotes the collection of projection bins g j , i.e., measured data. Here the projection matrix H is used to model both the attenuation and collimator response. The methodology to calculate this projection matrix is presented next. A. Projection model with non-uniform attenuation Consider projection bin g j for a given angular position θ of the detector array (see Figure 3). Without loss of generality, we present our methodology in a rotated xy coordinate system.

SPIE-IS&T/ Vol. 6814 68140V-3

Detector Array

g,

(x,y) S

rczy(j)

Figure 3. Integration ray for gj and intersection of the ray with a mesh element. The projection can be calculated by the following equation: E[ g j ] =

d detector





ddetector

f r ( x , y0 ) ⋅ e



µr ( x , y ) dy

y0

(4)

dy0

−∞

where the f r ( x , y ) and µr ( x, y ) are the functions of the emission image and the attenuation map, respectively, expressed in the rotated reference system, and d detector is the distance at which the detector array is located. Next, by taking advantage of the fact that the same mesh structure is used for both the emission image and absorption map, we can split the integration into a sum of sub-integrals, where each sub-integral is defined over each mesh element Tl (see Figure 3). That is, E[ g j ] =

∑ ∫ l

where

ray ( j )

denotes

the

line

⎡ ⎢ ⎢ ⎢Tl ∩ ray ( j ) ⎢⎣

that



f r ( x, y0 ) ⋅ e

goes

from

Ω j ,l = ray ( j ) ∩ ⎡⎣ y > IPc ( j , l )⎤⎦ is the segment of ray ( j )

− ∑ ⎢⎢

⎤ IPc ( j ,l )



k ⎢ Tk ∩Ω j ,l ⎣

µr ( x , y ) dy ⎥⎥ −

infinity

⎥ ⎦

to



y0

the

µr ( x , y ) dy

⎤ ⎥ dy0 ⎥ ⎥ ⎥⎦

detector

(5)

bin g j ,

perpendicularly,

between the detector face and the triangle Tl , and

IPc ( j , l ) = ( xˆc , j ,,l , yˆ c , j ,l ) is the closest intersection point to the detector between the triangle Tl and ray ( j ) . Equation (5) is based on the fact that the contribution of each triangle to a certain detector bin is attenuated by the path from the given triangle to the detector face. Further, with dense mesh elements (i.e., small in area), Eq. (5) can be approximated by assuming that the emission takes place in the center of the intersection between the ray and the triangle. This approximation is commonly used in pixel based projection models [25][26]. Thus, Eq. (5) can be approximated as follows:

SPIE-IS&T/ Vol. 6814 68140V-4

E[ g j ] ≈

∑ l

⎡ ⎤ ⎡ − ∑ ⎢ ∫ µr ( x , y ) dy ⎥⎥ ⎢ {k |Tk ∩Ω j ,l ≠∅} ⎢⎣⎢Tk ∩ray ( j) ⎦⎥ ⎢e ⎢ ⎣



⋅e

1 2



µ r ( x , y ) dy

Tl ∩ ray ( l )



Tl ∩ray ( j )

⎤ ⎥ f r ( x , y0 )dy0 ⎥ ⎥ ⎦

(6)

Now, using mesh modeling and linear interpolation over mesh elements Eq. (5) can be compactly rewritten as follows: E[ g j ] =

⎡ −{ ∑ ∑l ⎢⎢e ⎣

}

M j ,k µˆ k

⋅e

k |Tk ∩Ω j ,l ≠∅

1 − M j ,l µˆ l 2

⋅ M j ,l fˆl

⎤ ⎥ ⎥⎦

(7)

where µˆ l and fˆl are vectors formed by the nodal values of the attenuation and emission map for the triangle Tl , and M j ,l is defined as follows:

M j ,l = d j ,l

⎛ ⎜ ⎜ ⎜ ⎜ ⎝

( xˆ ( yˆ

c , j ,l c , j ,l

+ xˆ f , j ,l ) 2 ⎞

T



+ yˆ f , j ,l ) 2 ⎟ al

(8)

⎟ ⎟ ⎠

1

B. Projection model with non-uniform attenuation compensation and collimator response In this section we describe how to include the collimator response in the projection model. A detailed theoretical analysis of parallel hole collimators can be found in [27]. It is generally accepted [13] that the point spread function (PSF) can be modeled as a Gaussian curve where the full with at a half maximum (FWHM) is linearly dependent on the distance (in our experiments, σ c (d ) = 2.354 (σ 0 + d ⋅ σ d ) ). In our approach, we will approximate this by casting a cone of rays (see Figure 4) from each detector bin and computing the integral over each ray, followed by summation with an appropriate weighing factor. d

j) 'ZsC S

I I I

j

/ // Surface

PSF for differeist distances from the detector

ofthe

Detector Array

cc. /

Focal Point

Figure 4. Cone of rays that approximates the parallel hole collimator response.

SPIE-IS&T/ Vol. 6814 68140V-5

Since at the face of the detector the PSF FWHM is equivalent to 2.354σ 0 , the inclination of the rays is adjusted by moving the focal point of the rays further from the detector face by d focal = d detector + σ 0 σ d . Next the rays of the cone cast at angle αi are defined as follows: ⎛

i



N rays

αi = arctan ⎜⎜ 3σ d

⎞ ⎟ + 1 ⎟⎠

(9)

and their associated weights are defined as follows: −

e

wi =

1+ 2

1 ⎛ 3i ⎞ ⎜ ⎟ 2 ⎜⎝ N rays +1 ⎟⎠

∑e

N rays



2

1 ⎛ 3b ⎞ ⎜ ⎟ 2 ⎜⎝ N rays +1 ⎟⎠

(10)

2

b =1

In summary the forward projection model can be described as follows: E[ g j ] =

∑ w ∑e

N rays





{

}

k∈ k |Tk ∩Ω j ,l ,i ≠∅

i

i =− N rays

M j ,k ,i µˆ k

e

1 − M j ,l ,i µˆ l 2

M j ,l ,i fˆl

(11)

l

where M j ,l ,i now incorporates the dependence on which ray of the cone is being acquired ( αi ). From Eqs. (7) and (11) it is evident that a linear relationship exists between the projected data and nodal values. Thus, the projection bin g j and the nodal values f ( pq ) , q = 1 N , can be described as follows:



E[ g j ] =

∑w N

j ,q

⋅ f (p q )

q =1

where w j ,q are the elements of the matrix H .

IV. SIMULATION RESULTS AND DISCUSSIONS In our evaluation study, the 4D NURBS-based cardiac-torso (NCAT) 2.0 phantom [28] as used to generate the emission and attenuation distribution. The SIMIND Monte Carlo package [29] was used to simulate gated SPECT imaging with Tc99m labeled sestamibi as the imaging agent. The SPECT system modeled was the Picker Prism3000 with a low-energy high-resolution (LEHR) collimator. The projections were 64 by 64 bins with a pixel size of 0.634 cm. For a circular camera rotation of a 28.5 cm radius, 64 projection sets were collected. The simulation was performed at a highlevel (100 million) of total counts, which can be considered noiseless in order to test the accuracy of the developed projection model. The scattered counts were about 34% of the total counts and for the image reconstruction simulations, the scattered counts have been subtracted from the measured data to mimic ideal scatter compensation. In reconstruction we used the expectation-maximization (EM) algorithm [30], assuming Poisson likelihood. We show reconstructed images in Figure 5 for the phantom (shown in top-left). For comparison, we also show images obtained by filter-back projection (FBP) (top,-right), and by a pixel-based ML-EM algorithm with attenuation and collimator response correction [25][26] (bottom-left). By using noiseless projection data, our goal was to demonstrate comparable accuracy of our projection operator and that of pixel-based one. In Figure 6 we show a plot of the peak signal-to-noise ratio ( PSNR ) vs the iteration number. The PSNR is calculated as: ⎛



N max(f ) 2 ⎟ dB 2 ⎜ ⎟ ⎜ ⎟ f − fˆ

PSNR = 10log ⎜ ⎝



where fˆ is the estimated image, f is the original phantom image, and N is the number of pixels.

SPIE-IS&T/ Vol. 6814 68140V-6

(12)

Figure 5. Reconstructed images of the phantom (top-left) by different methods: FBP reconstruction (topright), pixel based ML-EM reconstruction (bottom-left), and mesh based ML-EM reconstruction (bottomright). PSNR . number of Itereflons 20 9.5

19

8.5

18

17

6.5

16

/

15.5

—*——.—-

15

—0—!

5

Mesh EM with scatter correction and attenuation compensation Mesh EM with scatter correction, deblurring and attenuation compensation Pixel EM with scatter correction and attenuation compensation Pixel EM with scatter correction, deblurring and attenuation compensation ! 10 15 20 25 30 35 40 45 50 number of iterations

Figure 6. Peak signal-to-noise ratio vs. iterations for the FBP, pixel ML-EM and mesh ML-EM.

SPIE-IS&T/ Vol. 6814 68140V-7

V. CONCLUSIONS AND FUTURE WORK We presented a methodology to compute a projection operator for a mesh representation model with several major data degradation factors incorporated, namely object attenuation and depth-dependent blur in detector-collimator spatial response. The projection operator was derived based on a 2D vector-valued content-adaptive mesh model for both attenuation and emission images. The methodology presented here is based on SPECT but can be easily extended to other modalities. We tested our derived projection operator by reconstructing images from a Monte Carlo simulation of SPECT. In future work plan to further evaluate this approach using noisy projection data, and extend it to 3D volumetric data.

VI. REFERENCES [1] Y. Yang, M. N. Wernick, and J. G. Brankov, “A fast algorithm for accurate content-adaptive mesh generation,” Proceedings of the IEEE International Conference on Image Processing (ICIP-01), vol. 3, pp. 868-871, 2001. [2] Yongyi Yang, Miles N. Wernick, and Jovan Brankov, “A fast algorithm for accurate content-adaptive mesh generation,” IEEE Transactions on Image Processing, vol. 12, pp. 866-881, 2003. [3] Jovan G. Brankov, Yongyi Yang, and Miles N. Wernick, “Content-adaptive 3D mesh modeling for representation of volumetric images,” Proceedings of the IEEE International Conference on Image Processing, vol. 3, pp 849-852, 2002. [4] Jovan G. Brankov, Yongyi Yang, and Miles N. Wernick, “Accurate mesh representation of vector-valued (color) images,” Proceedings of the IEEE International Conference on Image Processing, 2003. [5] Jovan G. Brankov, Yongyi Yang, and Miles N. Wernick, “Tomographic image reconstruction using content-adaptive mesh modeling,” Proceedings of the IEEE International Conference on Image Processing (ICIP-01), vol. 1, pp. 690693, 2001. [6] Jovan G. Brankov, Yongyi Yang, and Miles N. Wernick, “Content-adaptive mesh modeling for fully-3D tomographic image reconstruction,” Proceedings of the IEEE International Conference on Image Processing, vol. 2, pp. 621-624, 2002. [7] Jovan Brankov, Yongyi Yang, and Miles N. Wernick, “Tomographic image reconstruction based on a contentadaptive mesh model,” IEEE Transactions on Medical Imaging, vol. 23, pp. 202-212, 2004. [8] Jovan G. Brankov, Y. Yang, M. V. Narayanan, M. N. Wernick “Motion-compensated 4D processing of gatet SPECT perfusion studies,” Conference Record of the 2002 IEEE Nuclear Science Symposium & Medical Imaging Conference, vol. 3, pp. 1380-1384, 2002. [9] Jovan G. Brankov, Yongyi Yang, and Miles N. Wernick, Spatio-temporal processing of gated SPECT images using deformable mesh modeling,” Medical Physics , vol. 32, no. 9, pp. 2839-2849, 2005. [10] Ricard Delgado Gonzalo and Jovan G. Brankov, “Mesh model based projection operator for emission Tomography,” Conference Record of the 2007 IEEE Nuclear Science Symposium & Medical Imaging Conference, 2007. [11] Jovan G. Brankov, Yongyi Yang, Richard M. Leahy, and Miles N. Wernick, “Multi-modality tomographic image reconstruction using mesh modeling,” Proceedings of the IEEE/NIH International Symposium on Biomedical Imaging, pp. 405-408, 2002. [12] R. Leahy and J. Qi, "Statistical approaches in quantitative positron emission tomography," Statistics and Computing, vol. 10, pp. 147-165, 2000. [13] M. N. Wernick et al., Emission Tomography, Elsevier, 2004. [14] Y. Bresler, J. A. Fessler, and A. Macovski, “A Bayesian approach to reconstruction form incomplete projections of a multiple object 3D domains,” IEEE Trans. Patt. Anal. Mach. Intell. vol. 11, pp. 840-858, 1989. [15] G. S. Cunningham, K. M. Hanson, and X. L. Battle, “Three dimensional reconstruction from low-count SPECT data using deformable models,” IEEE. Nuclear Science Symposium, Volume 2, pp. 1469-1473, 1997. [16] K. M. Hanson, G. S. Cunningham, G. R. Jennings, and D. R. Wolf, “Tomographic reconstruction based on flexible geometric models,” Proc. of IEEE Int. Conf. on Image Processing, vol. 2, pp. 145-147, 1994. [17] D. J. Rossi and A. S. Wilsky, “Reconstruction from projections based on detection and estimation of objects,” IEEE Trans. Acoust. Speech, Signal Processing, vol. 32, pp. 886-906, 1984. [18] K. Shmueli, W. R. Brody, and A. Macovski, “Estimation of blood vessels boundaries in X-ray images,” SPIE Conf. Digital Radiology, Sept. 14-16, 1981. [19] R. M. Lewitt, “Alternatives to voxels for image representation in iterative reconstruction algorithms,” Phys. Med. Biol., vol. 37, pp.705-716, 1992.

SPIE-IS&T/ Vol. 6814 68140V-8

[20] S. Matej and R. M. Lewitt, “Efficient 3D grids for image reconstruction using spherically-symmetric volume elements,” IEEE Trans. Nucl. Sci., vol. 42, pp.1361-1370, 1994. [21] A. Sitek, R. Husman, and G. Gillberg, “Tomographic reconstruction using an adaptive tetrahedral mesh defined by a point cloud,” IEEE Transactions on Medical Imaging, vol. 25, pp. 1172-1179, 2006. [22] Y. Zhang, J. A. Fessler, N. H. Clinthorne, W. L. Rogers, “A hybrid-grid parameterization method for SPECT reconstruction,” J. Nuc. Med. (Abs. Book), 36(5):172, May 1995. [23] R. Floyd and L. Steinberg, “An adaptive algorithm for spatial gray scale,” SID Int. Sym. Digest of Tech. Papers, 1975. [24] Dudgeon, D. E. and Mersereau, R. M., Multidimensional Digital Signal Processing: Prentice-Hall, 1984. [25] E. S. Chornoboy et al., IEEE Trans. Med. Imag. vol. 9, no. 1, pp. 99–110, Mar. 1990. [26] A. W. McCarty et al., IEEE Trans. Med. Imag. vol. 10, no. 3, pp. 426–436, Sep. 1991. [27] C. Metz et al., Phys .Med. Biol. vol 25. pp. 1059-1070,1980. [28] W. P. Segars, Development of a new dynamic NURBS-based cardiac-torso (NCAT) phantom, Ph.D. thesis, The University of North Carolina, 2001. [29] M. Ljungberg and S-V. Srand, Comp. Meth. Prog. Biomed., vol. 29, pp. 257–272, 1989. [30] K. Lange and R. Carson, Journal of Computer Assisted Tomography, vol. 8, pp. 306-316, 1984.

SPIE-IS&T/ Vol. 6814 68140V-9