MEDITOMO: a high performance software package for 3D SPECT imaging

Abstract

We describe a parallel software library, named MEDITOMO, designed for processing MEDIcal images obtained by SPECT (single photon emission computed tomography) TOMOgraphic systems. MEDITOMO is the core library of the PSE (problem solving environment) MEDIGRID, oriented to medical imaging analysis, which the authors are currently developing. MEDIGRID is employed in a Grid-computing infrastructure involving clinical departments and research institutes. The algorithms of MEDITOMO are the standard ones that are usually applied in the SPECT analysis, i.e. the conjugate gradient and the expectation maximization. The main contribution of this work concerns the introduction of the total variation seminorm as the edge-preserving regularization in both algorithms and the development of the parallel software library. Experiments carried out on synthetic and clinical data are shown.

Keywords:

• medical imaging
• inverse problems
• parallel algorithms
• edge-preserving regularization
• mathematical software

2000 AMS Subject Classification :

• 65R30
• 65R32
• 65Y05
• 65Y20
• 65F22

Acknowledgements

We thank Prof. Mario Bertero for his suggestions. We thank Andreas Robert Formiconi and Piero Calvini who gave us the prototype software and collaborated with the development of MEDITOMO. We thank the Careggi Hospital in Florence, which provided the data we used in the numerical experiments, and also thank the help of its doctors in nuclear medicine. The work is partially supported by the national MIUR-PRIN (2004–2006, 2006–2008) projects and by the MIUR-PON S.Co.P.E. project (2006–2008). The software library has been developed by the Department of Physiopathology of the University of Florence and the Department of Physics of the University of Genoa.

Notes

^†The likelihood function is replaced by the logarithm of the likelihood function, whose maxima obviously coincide with the maxima of the likelihood function.

^‡For the Bayes’ Law.

^†The values of J, M, L, a, b, and c depend on the acquisition system.

^†We use the same notation adopted for the discretization of the integral operator describing the blurred Radon transform introduced in Section 2. Indeed, even if the operators are different, the discrete problems arising in the case of the Gaussian or the Poisson distribution of the noise have the same computational characteristics. More important, the corresponding algorithms have the same computational kernels, is the matrix–vector product involving the matrix K.

^†This approach takes into account the specific structure of the matrix K. Due to its high sparsity (between 5 and 10%), the matrix K is stored in a look-up table 8. The look-up table is given in closed form and we can access only each row of K that corresponds to a given angular projection over the whole dataset.

*Each module has a name of the form XXtomo_YYY_ZZ, where

•	XX denotes the model (3d or 2d+1);
•	YYY denotes the geometry of the acquisition system (parallel par or fan fan);
•	ZZ denotes the algorithm (conjugate gradient cg, expectation maximization emos, or fixed point fp+tv).

^† In-vivo data were kindly provided by the Careggi Hospital in Florence.