A hybrid region growing tumour segmentation method for low contrast and high noise Nuclear Medicine (NM) images by combining a novel non-linear diffusion filter and global gradient measure (HNDF-GGM-RG).

Poor spatial resolution and low signal-to-noise ratio (SNR) along with the finite image sampling constraint make lesion segmentation on Nuclear Medicine (NM) images (e.g., PET-Positron Emission Tomography) a challenging task. Since the size, signal-to-background ratio (SBR) and SNR of lesion vary within and between patients, performance of the conventional segmentation methods are not consistent against statistical fluctuations. To overcome these limitations, a hybrid region growing segmentation method is proposed combining non-linear diffusion filter and global gradient measure (HNDF-GGM-RG). The performance of the algorithm is validated on PET images and compared with the 40%-fixed threshold and a state-of-the-art active contour (AC) methods. Segmented volume, dice similarity coefficient (DSC) and percentage classification error (% CE) were used as the quantitative figures of merit (FOM) using the torso NEMA phantom that contains six different sizes of spheres. A 2:1 SBR was created between the spheres and background and the phantom was scanned with a Siemens TrueV PET-CT scanner. 40T method is SNR dependent and overestimates the volumes ( ≈ 4.5  times ) . AC volumes match with the true volumes only for the largest three spheres. On the other hand, the proposed HNDF-GGM-RG volumes match closely with the true volumes irrespective of the size and SNR. Average DSC of 0.32 and 0.66 and % CE of 700% and 160% were achieved by the 40T and AC methods respectively. Conversely, average DSC and %CE are 0.70 and 60% for HNDF-GGM-RG and less dependent on SNR. Since two-sample t-test indicates that the performance of AC and HNDF-GGM-RG are statistically significant for the smallest three spheres and similar for the rest, HNDF-GGM-RG can be applied where the size, SBR and SNR are subject to change either due to alterations in the radiotracer uptake because of treatment or uptake variability of different radiotracers because of differences in their molecular pathways.

Introduction

Positron Emission Tomography (PET), a Nuclear Medicine (NM) based medical imaging modality, provides functional images of the tumour with quantitative information of the radiotracer uptake. Robust, precise, and reproducible segmentation of the tumour on PET images is vital for machine learning based accurate diagnosis, treatment planning and response assessment. However, poor spatial resolution due to the partial volume effect (PVE) resulting from finite resolution of the PET camera and low signal-to-noise ratio (SNR) make delineation of the tumour on PET images a challenging task. Functional volume segmentation of clinical images conventionally relies on manual delineation of region of interest (ROI) by expert radiologist either on PET images directly by going through each 2D slice to obtain precise boundary information or using co-registered anatomical images (CT or MRI) (Hogenauer et al., 2016; Trigonis et al., 2014). The accuracy of the manual ROI delineation on the PET images is very much dependent on the intensity window chosen for visualization (Chen et al., 1998). The clinician also need to mentally reconstruct a structure in 3D since only 2D image slices can be seen which is time consuming and tedious work (Shareef et al., 1999). Irrespective of the method employed, manual delineation of ROIs is always laborious, highly operator dependent, requires significant knowledge of local anatomy and may be expected to produce significant intra and inter observer variability (Krak et al., 2005; Maroy et al., 2008) which makes manual delineation less reproducible (Foster et al., 2014).

To overcome the limitations of manual delineation method and make the segmentation method robust, precise, and reproducible, a number of computer-aided segmentation techniques have been proposed for PET images, that can usually be classified as semiautomatic or fully automatic methods (Foster et al., 2014; Li et al., 2008; Zaidi and El Naqa, 2010). Since the intensity in PET images directly correlates with the radiotracer uptake, intensity threshold based segmentation method is widely used due to its simplicity and predictability. However, due to poor spatial resolution, the histogram of PET images are not bimodal and hence determination of an optimal threshold is a challenging task. Inaccurate threshold may result in either overestimation or underestimation of the tumor volume. Threshold intensity (IThreshold) is generally defined as a percentage of the maximum intensity (IMax) within the roughly delineated volume and is given bywhere varies from 0 to 100%. A fixed value of is usually used for threshold based segmentation. A value of 40% for is proposed by Erdi et al. (1997) for lesions larger than 4 ml with a signal to background ratio greater than five. Several other researchers have also been using 40% threshold (40T) for tumour segmentation on PET images (Hong et al., 2007; Nestle et al., 2005; Schinagl et al., 2007; Wanet et al., 2011). However, the same fixed threshold value may not be optimal across different situations where statistical fluctuations are expected either due to treatment or progression of the disease. To overcome these challenges, several adaptive threshold based methods have been proposed over the years (Daisne et al., 2003; Jentzen et al., 2007; Schaefer et al., 2008). One of the main drawback of the adaptive threshold method is that the parameters need to be calibrated for each individual scanner and reconstruction protocol and indirectly dependent on the maximum intensity, IMax. Gradient based segmentation method is another alternative to the threshold based method (Geets et al., 2007; Wanet et al., 2011). However, it has found limited applicability in PET image segmentation because it requires substantial amount of preprocessing to get rid of the noise. Several region growing methods have also been proposed for automatic PET image segmentation method (Day et al., 2009; Li et al., 2008; Tan et al., 2017). The criteria of inclusion of voxel in the region for these region growing methods generally depends on the statistics of the surrounding voxels and thus the delineated ROIs vary with the statistical fluctuations. On the other hand, the performance of these methods is very much dependent on the selection of the seed point.

This paper proposes a fully automatic and easily implementable hybrid region growing tumour segmentation method combining a novel non-linear diffusion filter and global gradient measure (HNDF-GGM-RG) primarily based on the works of Tamal (2017) and Hojjatoleslami and Kittler (1998). The method overcomes the limitations of the existing region growing methods and does not require calibration like adaptive threshold based method. A stopping criteria for the proposed region growing method is automatically selected based on the global gradient information of the region of the smoothed images rather that the local neighborhood gradient information. The performance of the method is compared with the most widely investigated 40% fixed threshold and a state-of-the-art Chan and Vese model of active contour (AC) methods (Chan and Vese, 1999; Kass et al., 1987) using lesions with different size and SNR.

Materials and methods

Non-linear diffusion filter (NDF)

The stopping criteria of the diffusion process or diffusion coefficient is generally determined by the gradient intensity information (Perona and Malik, 1990). However, for low resolution images (e. g., NM images), only the gradient intensity information is not sufficient to define the edge specifically in the presence of high noise. Since the orientation of the gradient vector is systematic on the boundary and random for noise, it has been showed that utilization of the gradient vector orientation coherence within a fixed size window along with the gradient intensity can provide superior results compared to the conventional edge based diffusion filtering method (Tamal, 2017). In this case the diffusion coefficient is calculated as (Tamal, 2017)where is the gradient magnitude of the original image smoothed with a Gaussian kernel of width,. k is a parameter that is used to control the level of smoothing. is the sum of the cosine angles of gradient vectors of the smoothed image within a predefined window and is calculated for each 2D slice to form a 3D volume of the cosine angles. In such case, is represented in 2D by and determined as indicates the location of the current pixel and is the index of the pixel within the window with and specifying the size of the window.

Peripheral contrast measure based region growing (PCM-RG)

The central idea of the peripheral contrast measure based region growing (PCM-RG) method is that the region is grown by pixel aggregation by using similarity and discontinuity measures (Hojjatoleslami and Kittler, 1998). In this method, once an arbitrary pixel is selected as a seed point, amongst all the boundary pixels of the seed point, the pixel with the highest grey level value can only be included in the current region. In this way, the pixels with the higher grey levels are included first in the regions followed by the absorption of the boundary pixels with the monotonically lower and lower gray levels. During the inclusion of each single pixel, current boundary (CB) and internal boundary (IB) are determined. CB is the set of pixels adjacent to the current region and belongs to the background (Figure 1(a)). On the other hand, IB is defined as the boundary produced by the set of connected outermost pixels or edge pixels and belongs to the current region (Figure 1(a)). Both CB and IB evolve as the current region grows.

Figure 1

(a) Region growing method with current boundary (CB), internal boundary (IB) and current region. (b) a 26-connected voxel neighborhood system.

Hybrid non-linear diffusion filter and global gradient measure based region growing (HNDF-GGM-RG)

The above region growing method cannot provide accurate segmented ROI when the signal to background ratio (SBR) or contrast is low and noise is high as in the case of PET images. The method is also required to be modified and adapted for 3D medical images. Instead of calculating PCM directly from the original image which is unreliable, a method is proposed to use the smoothed images filtered using both the gradient vector orientation and gradient magnitude information. A global gradient measure based stopping criteria is then proposed to apply on the filtered images. The method will be called as HNDF-GGM-RG and is defined aswhere are the voxels on the edge of ROI (IB) and are the voxels on the boundary of the ROI (CB) respectively of the filtered image. The method requires two passes to provide the final segmented results. In the first pass, the region growing process is run to include a predefined number of voxels (VInitial) which is determined by the roughly delineated ROI cossntaining the tumour. Voxel containing the maximum intensity, IMax within the ROI is selected as the seed point for the region growing method. For 3D a 26-connected voxel neighborhood system is implemented to find the boundary voxels. In a 26-connected system, each voxel is neighbor to other voxel if they share same faces, edges or corners. In the system, each voxel with coordinates is connected to a voxel with coordinate . The 26-connected system is graphically depicted in Figure 1(b). The same 26-connected voxels are used to find the global gradient for each stage of the region growing process. As each individual voxel gets included in the region, the current boundary and the internal boundary of the ROI evolves (Figure 2). For the inclusion of each individual voxel, HNDF-GGM is calculated using Eq. (4). Upon completion of the region growing covering the whole roughly delineated ROI, the final stopping criteria (i.e., the number of voxels to be included in the region) is defined by the number of voxels (VFinal) that corresponds to the maximum HNDF-GGM (Figure 2). Once VFinal is determined, the region growing method is run again (second pass) starting from the same seed point with the stopping criteria being defined by the VFinal. Since the HNDF-GGM uses the difference between two neighboring boundaries rather than that of two neighboring voxels, it is less sensitive to noise compared to the other gradient magnitude based methods that only rely on gradient measure in a particular direction (e.g., x, y and z directions). Moreover, the proposed method does not use the information of the image statistics and hence is robust to change in size, contrast and SNR.

Figure 2

Illustration of Boundary evolving process of the proposed HNDF-GGM-RG method. Light blue, light grey, dark blue, dark grey and black colours represent original ROI, background, current boundary (CB), current region and internal boundary (IB) respectively. (a) True ROI (light blue) in a background (light grey). (b) Region growing method starting from a seed point (dark grey) and its boundary (dark blue). (c) to (h) show how the current boundary (dark blue) and edge or the internal boundary (black) of the current region evolves as the region (dark grey + black) grows. For every region starting from only one voxel or seed voxel until all voxels within the roughly delineated area (VInitial) are included, HNDF-GGM is calculated for each step using Eq. (2). The graph shows the HNDF-GGM against the number of voxels. Voxel number corresponding to the maximum HNDF-GGM (red filled circle) is the final number of voxel (VFinal) that should be included in the final segmented region starting from the same seed point.

The HNDF-GGM-RG algorithm for PET images is summarized below:

B = HNDF-GGM-RG (I, ROI).

Input:

I: PET image;

ROI: a rough region of interest manually delineated to enclose the tumor.

Intermediate Output:

NDF-I: filtered image after applying filter on the original PET image

First pass:

Determine the total number of voxel (VInitial) from roughly delineated ROI to include in the region;

Define the voxel (VMax) as seed point that contains maximum intensity, IMax.

Iteration:

For k = 1: 1: VInitial

Obtain global gradient measure for inclusion of each voxel: HNDF-GGM = HNDF-GGM-RG (NDF-I, ROI,VMax);

End

Determine final voxel number (VFinal) corresponding to the highest HNDF-GGM.

Second pass:

Start region growing at VMax

Iteration:

For k = 1: 1: VFinal

Update the segmented image: B = HNDF-GGM-RG (NDF-I, ROI,VMax);

End

Output:

B: binary image with segmented region containing VFinal number of voxels.

Validation

The torso NEMA phantom containing six spheres with diameter of 10, 13, 17, 22, 28 and 37 mm (0.52, 1.15, 2.57, 5.58, 11.49 and 26.52 cm3 volume) was filled with 1668 KBq/ml 18F solutions (Figure 3). The uniform background was filled with 838 KBq/ml 18F solution to yield 2:1 SBR or contrast between the homogeneous hot spheres and the colder uniform background. Low contrast of 2:1 was chosen for validation because for higher contrast levels the performance of most the segmentation methods converges.

Figure 3

(a) Torso NEMA phantom with six fillable spheres and a cold insert in the middle. (b) The dimensions of the shperes.

The phantom was scanned on the TrueV PET-CT scanner (Siemens, USA) for 120 min in a 3D mode which provides 109 image planes or slices covering a 21.6 cm axial FOV (field of view). The 3D images were then reconstructed using a 256 × 256 × 109 matrix with voxel dimensions of 2.67 × 2.67 × 2.00 mm. Five different scan durations (900, 1200, 2000 and 4000 s corresponding to 15, 20, 33.3 and 66.6 min) to represent different levels of SNR were reconstructed using OSEM reconstruction algorithm with 4 iterations and 21 subsets. To reconstruct five different non-overlapping and overlapping realizations for all durations the starting time of each static frame were shifted. A 4-mm FWHM (full width at half maximum) Gaussian filter was applied on all the reconstructed images after applying radioactivity decay correction. The boundaries calculated from the known diameter and position of each sphere was then used to estimate the true volume of interest (VOITrue).

Validation measures

The performance (robustness, accuracy and reproducibility) of the HNDF-GGM-RG method was compared with the most widely investigated 40% fixed threshold and active contour (AC) methods for different statistical situations (size and SNR). For the purpose of the validation, along with the segmented volumes, Dice Similarity Coefficient (DSC) and Classification Error (CE) were also analyzed for each segmentation method. DSC provides quantitative measures of spatial overlap index with the VOITrue. It can be used to evaluate segmentation accuracy. DSC is given bywhere ∩ is the intersection and VOISeg is the segmented volume. DSC value 0 indicates complete non-overlap and 1 indicates complete match or overlap between two volumes.

Percentage classification error (CE) is defined aswhere PCE (positive classification errors) refers to the background voxels that are classified as voxels belonging to the sphere. On the other hand, NCE (negative classification errors) refers to the voxels within the sphere belonging to the background. Higher CE value is indicative of poor segmentation accuracy.

Results

Representative segmentation results of one realization for 40T, AC and HNDF-GGM-RG along with the true volumes for 13 and 28 mm sphere for all acquisition durations are shown in Figure 4. Average volumes of five realizations along with their standard deviations for all six spheres and acquisition durations are shown in a bar graph in Figure 5.

Figure 4

Segmentation results of 40T, AC and HNDF-GGM-RG on one slice of the NEMA phantom. Top row: 13 mm diameter sphere. Bottom row: 28 mm diameter sphere for 900 s (a, e), 1200 s (b, f), 2000 s (c, g) and 4000 s (d, h) acquisition durations representing different SNR.

Figure 5

Mean segmented volumes (cm3) of five realizations for all six spheres and acquisition durations along with the VOITrue. The error bar represents standard deviation of the means. Each of the panel (a to e) represents mean and standard deviation of the segmented volumes for one sphere size starting from the smallest.

For 10 mm sphere, the average segmented volumes across different acquisition durations for 40T, AC and HNDF-GGM-RG methods are 15.04, 7.03 and 0.93 cm3 generating average error of 2792%, 1252% and 78% respectively (in comparison to VOITrue of 0.52 cm3). For 13 mm sphere, 40T method segments the same roughly delineated ROI irrespective of acquisition durations (average volume of 19.70 cm3 with 1612.74% error). The mean volume for the AC method for this sphere is 6.58 cm3 with 472% error. The performance of HNDF-GGM-RG method significantly improves for 13 mm sphere with the average sphere volume of 1.73 cm3 with 50% error. The absolute average errors of for the 40T method are 995%, 490%, 277% and 188% with average segmented volumes of 28.14, 32.95, 43.39 and 76.39 cm3 for 17 mm, 22 mm, 28 mm and 37 mm spheres respectively. Since 40T method is not able to accurately segment the lesions from background for low contrast, the final segmented ROI is very much dependent on the roughly delineated ROI.

AC method provides mean segmented volume of 5.85, 5.71, 10.04 and 23.49 cm3 with 127%, 18%, 13% and 12% error respectively. For the NDF-GGM-RG method, the error values are 46%, 46%, 31% and 22% for the same spheres with average segmented volumes of 1.39, 3.01, 7.86 and 20.54 cm3. The standard deviations of volumes for the 40T method are higher compared to that of the AC and HNDF-GGM-RG method indicating that the volumes delineated by the 40T method are less reproducible specially for 22 mm, 28 mm and 37 mm spheres. One the other hand, for 13 mm and 17 mm spheres, the 40T method delineated similar ROIs as the roughly delineated ROIs and therefore, the standard deviations for these two spheres are much lower but with high error.

The average percentage errors across different sizes and SNR are approximately 131% and 45% for AC and HNDF-GGM-RG methods respectively with the error for the AC method is very much dependent on SNR if the size of the tumour is small . However, the 40T method is dependent on the noise irrespective of the size of the sphere.

Average dice similarity coefficient (DSC) of five realizations is shown for all three methods in Figure 6. DSC is always the highest for the HNDF-GGM-RG method for the smallest three spheres compared to the 40T and AC. For the 10 mm diameter sphere, DSC is very low for all methods. For the 13 mm and 17 mm spheres, average DSC for HNDF-GGM-RG method is 0.49 and 0.68 respectively. For the same two spheres, DSC for the 40T method is 0.11 and 0.15 and for the AC is 0.27 and 0.57 with the AC method being very much dependent on SNR (0.36–0.82 for 17 mm sphere). For the biggest three spheres (22 mm, 28 mm and 37 mm), the average DSC are, 0.65, 0.77 and 0.87 for the HNDF-GGM-RG method. 0.31, 0.46 and 0.56 for the 40T method. 0.72, 0.83 and 0.90 are the DSCs for the AC method. Though the DSCs are similar for both HNDF-GGM-RG and AC for the biggest three spheres, the differences between the average DSC for the proposed method between acquisition durations are negligible indicating lower dependency on the noise.

Figure 6

Mean dice similarity coefficient (DSC) of five realizations for all six spheres and acquisition durations. The error bar represents standard deviation of the means. Each of the panel (a to e) represents mean and standard deviation of the DSC for one sphere size starting from the smallest.

The average percentage classification error (CE) of five realization of all the six spheres are shown in Figure 7. CE decreases with the increase of the size of the spheres for all three methods. . CE values for the biggest five spheres (13 mm–37 mm) are 137%, 50%, 55%, 38% and 23% for the HNDF-GGM-RG method compared to 1689%, 1145%, 537%, 274%, and 193% of 40T and 515%, 179%, 61%, 33% and 19% of AC. Since, the 40T method always segment ROIs much bigger than the original volume, percentage CEs are always higher for this method.

Figure 7

Mean percentage classification error (CE) of five realizations for all six spheres and acquisition durations. The error bar represents standard deviation of the means. Each of the panel (a to e) represents mean and standard deviation of the CE for one sphere size starting from the smallest.

The CEs for the 40T method are dependent on the SNR and the standard deviation for each case is very high. The CE of the AC method is only dependent on the SNR for the smallest three spheres. On the other hand, the differences in CE for the HNDF-GGM-RG method are less than 15% across acquisition durations.

The one-way ANOVA test performed on all spheres of all three methods indicates that the methods are significantly different with . Two-Sample t-Test with AC and HNDF-GGM-RG also shows that these methods are significantly different for the smallest three spheres. However, for the 22 mm, 28 mm and 37 mm sphere, these two methods are not significantly different with the value of p being 0.14, 0.67 and 0.17 respectively.

Discussion

Accurate segmentation of suspected lesion in PET is very important for further analysis for diagnosis, treatment planning and response assessment of the disease. Threshold based method is widely used because of its simplicity to implement and predictability. However, dependency on the size and noise makes it less accurate specially when the contrasts and counts are low. This paper proposes a novel hybrid global gradient measure based region growing method for low contrast and high noise NM images. The performance of the proposed methods was compared with the most widely investigated 40T threshold and state-of-the art active contour (AC) methods for different sizes and acquisition durations using Torso NEMA phantom.

For the smallest three spheres, the 40T method segment the same roughly delineated ROI. These spheres are subject to high level of partial volume effects (PVE) arising from limited resolution of PET camera and finite image sampling (Soret et al., 2007). Because of that, the value of the IMax is much lower and thus 40% threshold is too low compared to the optimum threshold value. The correction for PVE is vital when the diameter of the lesion is 2–3 smaller the than resolution of the scanner typically measured as full width half maximum (FWHM) (Bettinardi et al., 2014). The typical resolution of the TrueV PET-CT scanner used for this study is approximately 6 mm FWHM and therefore, the segmented volumes using the 40T method are grossly mismatched. For the biggest three spheres, where PVE corrections are not essential, volumes of the segmented ROIs by the 40T method increase with the increase of acquisition duration. DSC decreases and percentage CE increase with the acquisition duration. This is because, noise is reduced with the increase of acquisition duration and at 4000 s the noise is the lowest. 40% threshold is not optimum in such cases and the percentage threshold needs to be increased to avoid the gross overestimation of the volumes. This investigation also indicates that the 40T method is not only dependent on the size but also on the noise since IMax is also dependent on the size and noise (Akamatsu et al., 2015). Previous study has confirmed the effects of size on the threshold based segmentation method (Erdi et al., 1997). However, the effects of noise on the threshold based segmentation method was not explicitly addressed in the literature.

For the proposed HNDF-GGM-RG and AC, the accuracy (volume, reproducibility, DSC and CE) of the segmented volumes also becomes better with the increase of the sphere size. However, the differences in accuracy are much lower compared to the 40T method and are much less dependent on the noise. For the smallest two spheres with 13 mm and 17 mm diameter, the accuracy of the segmented volumes is significantly improved using the HNDF-GGM-RG method compared to both 40T and AC. This is because, the proposed method does not dependent on the maximum intensity value IMax and conventional detection of the edge both of which vary with the size and noise or acquisition duration of the image.

Overall, the proposed HNDF-GGM-RG segmentation method segments volume with less error, high DSC and low CE compared to the 40% threshold across different sizes and acquisition durations and for the smallest three spheres compared to the AC method. The method also does not require any user intervention and easy to implement. It is also insensitive to the location of the seed point.

Conclusion

A novel fully automatic hybrid global gradient measure based region growing (HNDF-GGM-RG) segmentation method is presented in this study and the accuracy of the segmentation method was compared with the most widely investigated 40% fixed threshold and AC methods across different sizes and SNR for low contrast. The HNDF-GGM-RG method provides more accurate results and is robust to the changes in size and noise. The method can be applied where the size of the lesion and uptake can change due to treatment. The method can also be used for different radiotracers where variable uptake are expected due to differences in molecular pathways of the tracers. The method is also fully automatic and does not require to select the seed point. The higher accuracy, reproducibility and robustness of the segmented volumes by the proposed HNDF-GGM-RG method indicate that it can be implemented for routine clinical purposes.

Declarations

Author contribution statement

M. Tamal: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This study was funded by the Deanship of Scientific Research (DSR), Imam Abdulrahman Bin Faisal University, Dammam, Saudi Arabia (grant number 2018-033-Eng).

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.