Leader-follower clustering algorithm for automatic segmentation of cardiac PET studies J. M. Mateos-Pérez, Student Member, IEEE, J. J. Vaquero, Senior Member, IEEE, C. García-Villalba, L. Cussó, M. Desco Member, IEEE Abstract– Quantitative analysis of dynamic cardiac PET studies provides useful diagnostic information regarding different heart diseases. To correctly perform the kinetic analysis, tracer time-activity curves (TAC) must be precisely extracted from the chosen tissue on the imaging data. In this work we present an alternative clustering algorithm for segmentation of dynamic studies using a leader-follower clustering approach, in which the number of clusters is unknown.
I. INTRODUCTION
D
ynamic molecular imaging studies provide quantitative information on the spatial distribution of different tracers, as well as functional information regarding the different tissues [1]. The conventional kinetic parameters estimation approach involves a segmentation of the image in order to extract the time-activity curves (TAC) for each region. This can be done manually, but it is a process very timeconsuming, error-prone, non-repeatable and user-dependent. Therefore, several automatic methods have been proposed in the literature [2-8]. Clustering methods that take into account spatial information have also been developed [9]. The number of classes specified in the segmentation is a critical parameter for obtaining good results [10]. In this work we present an alternative clustering algorithm for segmentation of dynamic studies using a leader-follower clustering algorithm, in which the number of clusters is unknown and the different voxels are grouped by TAC similarity regardless of their spatial position. Pearson’s correlation score is used as the main metric for the algorithm. A simulated study and a cardiac swine 13NH3 study have been used to test the method.
Manuscript received November 14, 2011. This work was supported in part by the CENIT-AMIT Ingenio 2010, Ministerio de Ciencia e Innovación, RETIC-RECAVA, Ministerio de Sanidad y Consumo. The authors would also like to thanks Dr. Michael W. Dae and Jose Manuel Udias for providing the data for this work. J. M. Mateos-Pérez is with CIBERSAM, Madrid, Spain (e-mail:
[email protected]. Telephone: +34 91 426 50 67). J. J. Vaquero is with the Departamento de Bioingeniería e Ingeniería Aeroespacial, Universidad Carlos III de Madrid, Spain (e-mail:
[email protected]). C. García-Villalba is with the Unidad de Medicina y Cirugía Experimental, Hospital General Universitario Gregorio Marañón, Madrid, Spain (e-mail:
[email protected]). L. Cussó is with the Unidad de Medicina y Cirugía Experimental, Hospital General Universitario Gregorio Marañón, Madrid, Spain (e-mail:
[email protected]). M. Desco is with the Unidad de Medicina y Cirugía Experimental, Hospital General Universitario Gregorio Marañón, Madrid, Spain and the Departamento de Bioingeniería e Ingeniería Aeroespacial, Universidad Carlos III de Madrid, Spain (e-mail:
[email protected]).
II. MATERIALS AND METHODS A. Leader-follower clustering One of the most popular clustering algorithms is k-means. This algorithm finds the k mean vectors that describe a given dataset. Starting from random or pre-defined centroids, the algorithms explores the dataset iteratively until no further modification is made on the final clusters, or until a certain threshold is reached [11]. Value for k is critical for a good segmentation [10]: more clusters imply more different regions to be segmented. As is the case in medical imaging, sometimes it is not possible to know how many different regions are present, due to pathology, noise, spill-over and other imaging effects. Therefore, a low value for k will group together regions that should be separated. It is possible to implement a clustering algorithm that does not depend on setting an initial cluster number but uses a sensitivity threshold. The pseudocode for this algorithm is written as follows [11]: begin initialize η,θ w1 ← x do accept new x j ← arg min ||x – wj’|| if ||w – wj|| < θ then wj ← wj + ηx else add new w ← x w ← w/||w|| until no more patterns return w1, w2,… end
In this algorithm, THETA is a sensitivity threshold that can be used to make the algorithm more restrictive (more clusters will be created) or less restrictive (fewer clusters). B. Implementation Let G(TAC1,...,TACn) be a group of n voxels defined by n TACs, and let G be the mean TAC from the voxels belonging to that particular group and ρG j the mean correlation score for all the voxels inside that group, with a mean amplitude μ and a standard deviation σ . A voxel with a TACi will be included in a given group Gj if: 1.
The voxel is similar to the voxels already present in that group ρ (G j , TACi ) > ρ (Gk , TACi ), ∀j ≠ k and
their correlation score is greater than a given threshold, ρt , ρ (G j , TACi ) > ρt . 2. The amplitude for this particular TAC is higher than μ − σ . If these requirements are met, that voxel will be included into group j and the values of G j will be updated accordingly. The first condition checks that the voxel belongs to that particular group, and the second condition force the group to accept only less noisy voxels as time passes. If the voxel does not belong to any known group, a new one will be created. The algorithm was implemented using the Java programming language as a plug-in for the PMOD software package (PMOD Technologies Ltd., Zurich, Switzerland). All tests were carried out on a desktop workstation powered by an Intel Core 2 Quad CPU with 4 GB of RAM. The COLT v1.2.0 (CERN – European Organization for Nuclear Research, Geneva, Switzerland) library was used for computational purposes. C. Simulated study A simulated study using a cylinder with two differentiated regions was used to test the algorithm. The activity for each cylinder is defined by the following equations: C (t ) hot _ cylinder =
Fig. 2: frame from the pig study. The ventricles can be seen very clearly inside the myocardium due to their lack of activity. The lungs can be seen due to some residual uptake.
5s frames, 2 x 15s frames, 3 x 60s frames, 2 x 300 s frames). The tracer injections were performed as a bolus using an automatic injector. Images were reconstructed using a FBP algorithm at 2 mm/pixel with 4.8mm cut-off and no postfiltering and scatter correction with convolution substraction. A manual segmentation of both ventricles and the myocardium was obtained by a trained technician at our institution.
k1 k Ca + Ca (4 − 1 )e − k2 t , k2 k2
k1 = 0.01min −1 , k2 = 0.05min −1 C (t )cold _ cylinder =
k1 k Ca − 1 Ca e − k2 t , k2 k2
(1.1)
k1 = 0.05 min −1 , k2 = 0.01min −1 Where Ca = 12.4 μCi/ml. The matrix size is 175 x 175 x 61 voxels, with a voxel size of 1mm x 1mm x 2mm. 20 frames of 60 seconds each were simulated. FBP was used as reconstruction method.
Fig. 3: activities for segmented regions. Correlation has been calculated against the theoretical activities Fig. 1: a frame from the simulated study.
D. Swine study A cardiac swine 13NH3 study was also used to test the algorithm. 740 MBq were administered for each study and image acquisition was done using a matrix size of 128 x 128 x 47, with a voxel size of 2.34 x 2.34 x 3.27 mm. Twenty five dynamic frames were acquired for a total of 900 seconds (18 x
IV. DISCUSSION III. RESULTS A. Simulated study The results for the automatic segmentation and TAC extraction for the simulated study can be seen on Fig. 3. Note the high correlation between the extracted TAC and the theoretical curves for the different regions. A visual reference for the segmented areas can be seen on Fig. 4.
Fig. 4: (left) reference image for simulated study. (Center) Hot and (right) cold cylinder segmented by the algorithm (θ = 0.5).
B. Swine study The results for the automatic segmentation and TAC extraction for the swine study can be seen on Fig. 5. Fig. 6 shows the segmented areas along with the manual segmentation for comparison purposes. The time employed in the segmentation was 13.26 s.
Fig. 5: TACs automatically extracted from swine study and correlation with manually obtained TAC.
Fig. 6: (left) right ventricle, (center) left ventricle and (right) myocardium as segmented by the presented tool. The overlapped dotted line represents the ROI drawn by manual segmentation.
A precise algorithm for segmentation and TAC extraction is a fundamental first step when performing any kind of kinetic analysis. Due to the noisy nature of dynamic PET sequences, it is important to group together as many voxels belonging to the same region as possible in order to get a noise-free curve. In this paper we have presented a clustering algorithm which uses shape similarity as its metric. The fact that this clustering algorithm does not depend on a random prior, as k-means does, allows for repeatability. There is only one parameter that needs to be set by the user, θ, which will implicitly define the final number of clusters. The algorithm offers a visual reference that allows the user to assess the final regions segmented by the algorithm. Though the final visual result is quite noisy, it can be observed the high agreement between the extracted TAC and the theoretical curve (simulated study) and the manual segmentation (real study). The algorithm groups voxels by TAC similarity, independently of their location within the volume; a better result can potentially be obtained by implementing spatial constrictions. It has to be taken into account that a high θ will cause many clusters to be created. The visual result provides aid to assess whether the value has been set too high or too low. REFERENCES [1]
R. Weissleder and U. Mahmood, “Molecular Imaging,” Radiology, vol. 219, no. 2, pp. 316-333, 2001.
[2]
R. Di Paola et al., “Handling of Dynamic Sequences in Nuclear Medicine,” IEEE Transactions on Nuclear Science, vol. 29, no. 4, pp. 1310-1321, 1982.
[3]
D. C. Barber, “The use of principal components in the quantitative analysis of gamma camera dynamic studies.,” Physics in medicine and biology, vol. 25, no. 2, pp. 283-92, Mar. 1980.
[4]
A. Sitek, G. T. Gullberg, and R. H. Huesman, “Correction for ambiguous solutions in factor analysis using a penalized least squares objective.,” IEEE Transactions on Medical Imaging, vol. 21, no. 3, pp. 216-25, Mar. 2002.
[5]
R. Klein, R. S. Beanlands, A. Adler, and R. DeKemp, “Model-based factor analysis of dynamic sequences of cardiac positron emission tomography,” in Nuclear Science Symposium Conference Record, 2008. NSS ’08. IEEE, 2008, pp. 5198-5202.
[6]
H. M. Wu et al., “Factor analysis for extraction of blood timeactivity curves in dynamic FDG-PET studies.,” Journal of nuclear medicine : official publication, Society of Nuclear Medicine, vol. 36, no. 9, pp. 1714-22, Sep. 1995.
[7]
X. Zheng, G. Tian, S.-cheng Huang, and D. Feng, “A Hybrid Clustering Method for ROI Delineation in Small Animal Dynamic PET Images : Application to the Automatic Estimation of FDG Input Functions,” IEEE Transactions on Information Technology in Biomedicine, 2010.
[8]
H. Becker, “A Survey of Correlation Clustering,” in Advanced Topics in Computational Learning Theory, 2005, pp. 1-10.
[9]
J. Kim, W. Cai, D. Feng, and S. Eberl, “Segmentation of VOI from multidimensional dynamic PET images by integrating spatial and
temporal features.,” IEEE Transactions on Information Technology in Biomedicine, vol. 10, no. 4, pp. 637-46, Oct. 2006. [10]
K.-P. Wong, D. Feng, S. R. Meikle, and M. J. Fulham, “Segmentation of dynamic PET images using cluster analysis,” IEEE Transactions on Nuclear Science, vol. 49, no. 1, pp. 200-207, Feb. 2002.
[11]
R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification. Wiley-Interscience, 2001, p. 654.
