Joint Tumor Segmentation in PET-CT Images Using Co-Clustering and Fusion Based on Belief Functions

1). Term 1:

Similar to the original ECM, the first term of (7) denotes the restriction between clustering distortions and mass functions.

2). Term 2:

The second term of (7) is a spatial regularization defined in the framework of belief functions for enhancing local homogeneity of neighboring voxels. According to the spatial prior of a PET volume, matrix $M^{pt}={\left\{m_i^{pt}\right\}}_{i=1}^n$ that we want to learn can be regarded as a specific random field, matrix where each mass function $m_i^{pt}$ is a random vector in ${ℝ}^3$, and its distribution depends on the mass functions of adjacent voxels. Let $\Phi ={\{\Phi (i)\}}_{i=1}^n$ be a 3-D neighborhood system, where Φ(i) = {1, …, T} is the set of the T neighbors for voxel i. The corresponding mass functions of voxels in Φ(i) are $\left\{m_{i,1}^{pt},\dots ,m_{i,T}^{pt}\right\}$, while the feature vectors of these voxels are $\left\{X_{i,1}^{pt},\dots ,X_{i,T}^{pt}\right\}$. Then, based on the above assumption about M^pt, ${\sum }_{t\in \Phi (i)}\left[d_m^2\left(m_i^{pt},m_{i,t}^{pt}\right)\right]\left[d^2\left(X_i^{pt},X_{i,t}^{pt}\right)\right]$ of (7) quantifies the smoothness around voxel i, where $d^2\left(X_i^{pt},X_{i,t}^{pt}\right)$ denotes the distance between voxel i and its neighbor t in the feature space; while, $d_m^2\left(m_i^{pt},m_{i,t}^{pt}\right)$ measures the independence between $m_i^{pt}$ and $m_{i,t}^{pt}$, i.e., the inconsistency between the mass functions of voxel i and its neighbor t. Based on (4), we set $d_m^2\left(m_i^{pt},m_{i,t}^{pt}\right)=\left(m_i^{pt}-m_{i,t}^{pt}\right)J{\left(m_i^{pt}-m_{i,t}^{pt}\right)}^T$, where matrix J is defined by (5).

3). Term 3:

To describe image voxels, we include textural features as complementary information for voxel intensities. The challenge to this end is that a large amount of textural features can be extracted, but without prior knowledge concerning the most informative features. In addition, the quantified high-dimensional feature vectors may contain unreliable variables due to noise and limited resolution of PET images. To tackle this challenge, distance metric adaptation and/or feature selection (e.g., [61]) are desirable.

In our previous work [62], a supervised method has been proposed to learn a low-rank dissimilarity metric for improving the performance of distance-based classifiers on high-dimensional datasets containing unreliable and imprecise features. Distinct from this previous work, here our goal is to adapt distance metric for given data in an unsupervised learning protocol, so as to improve the performance of clustering algorithms. Therefore, we look for a matrix $D^{pt}\in {ℝ}^{p\times q}$ during clustering, under the constraint q ≪ p, by which the dissimilarity between any two feature vectors, say $X_1^{pt}$ and $X_2^{pt}$, can be quantified as

A desired matrix D^pt transforms the original feature space to a low-dimensional subspace, where discriminant input features play a more significant role than uninformative ones in calculating the dissimilarity. To find such a D^pt, the distortion between any $X_i^{pt}$ and cluster centroid ${\overline{V}}_j^{pt}$, i.e., $d^2\left(X_i^{pt},{\overline{V}}_j^{pt}\right)$ in the first term of (7), is represented by (8). Moreover, the spatial regularization is also used to adapt the distance metric, where $d^2\left(X_i^{pt},X_{i,t}^{pt}\right)$ is also quantified by (8). During the iterative minimization of (7), a large dissimilarity $d_m^2\left(m_i^{pt},m_{i,t}^{pt}\right)$ between $m_i^{pt}$ and $m_{i,t}^{pt}$ will indicate that D^pt obtained as current step is unfit. Then, it should be adjusted at the next step to reduce the dissimilarity between $X_i^{pt}$ and $X_{i,t}^{pt}$, so as to bring adjacent voxels closer.

To ensure the effectiveness of metric updating, the third term of (7) is defined as the sparsity regularization ||D^pt||_2,1 of matrix D^pt, i.e., $\mathcal{F}={\sum }_{i=1}^p\sqrt{{\sum }_{j=1}^q{\left(D_{i,j}^{pt}\right)}^2}$. It aims at selecting the most reliable input features to calculate distance d²(·, ·) in the feature space (which exists in the other three terms). By forcing rows of D^pt to be zero, the proposed method only selects reliable features for clustering, while the influence of unreliable features is controlled. Scalar λ is a tuning parameter that controls the influence of this regularization.

4). Term 4:

The last term of (7) is used to prevent the cost function being trivially solved with D^pt = 0, which collapses all the features vectors into a single point. Vectors ${\overline{X}}_{{\omega }_1}^{pt}$ and ${\overline{X}}_{{\omega }_2}^{pt}$ are two seeds for the positive tissue and the background, respectively, which can be predetermined automatically. We empirically set a constant parameter (i.e., = 1) for this term, mainly considering it serves to prevent D^pt = 0, in which case the penalty will be activated and the corresponding loss will get close to infinity. In addition, we experimentally found that this term is insensitive to changing parameters.

1). Initialization:

We initialize the mass functions (i.e. M^pt and M^ct) and the singleton cluster centers (i.e. V^pt and V^ct) for both PET and CT via the original ECM algorithm. Based on the initial mass functions that obtained in PET, a very limited number of tumor and background seeds are predefined. After that, the two feature vectors, i.e., ${\overline{X}}_{{\omega }_1}^{pt}$ (resp. ${\overline{X}}_{{\omega }_1}^{ct}$) and ${\overline{X}}_{{\omega }_2}^{pt}$ (resp. ${\overline{X}}_{{\omega }_2}^{ct}$), used in the last term of (7) are calculated as the barycenters of the tumor and background seeds, respectively. The output dimension, namely the number of columns in feature transformation matrix D^pt (resp. D^ct), is then determined by applying principle component analysis on all the feature vectors ${\left\{X_i^{pt}\right\}}_{i=1}^n$ (resp. ${\left\{X_i^{ct}\right\}}_{i=1}^n$). The initial D^pt (resp. D^ct) is constructed by the top 95% eigenvectors.

The optimization procedure alternates between three parts: the clustering in PET, the clustering in CT, and the fusion of M^pt and M^ct at current step. For the clustering in each monomodality, the optimization iterates between cluster assignment (i.e. M^pt or M^ct estimation) in the E-step, and both prototype determination (i.e. V^pt or V^ct estimation) and distance metric adaptation (i.e. D^pt or D^ct estimation) in the M-step.

2). Optimization in PET:

This procedure only relates to the minimization of the first term and the last term of (10), which is performed in an EM-like protocol.

3). Optimization in CT:

The adaptation of M^ct, V^ct, and D^ct, which relates to the last two terms of (10), is following the same way as that for M^pt, V^pt, and D^pt discussed above. To update M^ct at current step, M^pt obtained by Section III-C2 is utilized in (12).

4). Fusion of PET and CT:

It is worth noting that M^pt and M^ct obtained in Section III-C2 and III-C3 are in fact two independent pieces of evidence regarding the same target tumor. They are complementary, as PET and CT can provide, respectively, functional and anatomical information. In addition, they are also consistent due to the soft context term defined by (9). Therefore, in this step, M^pt is further updated by fusing it with M^ct via the Dempster’s rule (i.e. (3)). The updated M^pt is then used as the initialization for the optimization (i.e. Section III-C2) in the next iteration.

The whole optimization procedure will not terminate the alternation between the steps described in Section III-C2, III-C3, and III-C4, until the value of (10) has no significant change between two consecutive iterations.

1). Illustrative Results of Co-Segmentation:

An example to illustrate the co-segmentation performance of the proposed method is shown in Fig. 2, where each row represents a different slice in the axial plane of the same tumor. In the first and second column, contours delineated in PET (green) and CT (magenta), are compared respectively to the ground truth (blue); while, in the last column, all of them are overlaid in the fused images. As can be seen, segmentation in PET is in consistent with that in CT, while they are not identical. The 3-D tumor volumes in PET and CT are further shown in Fig. 3, where the green, magenta, and orange regions represent, respectively, true positive voxels, false positive voxels, and false negative voxels, as compared to the ground truth.

2). Co-Segmentation vs. Mono-Modal Segmentation:

To demonstrate the effectiveness of the proposed co-segmentation method, its performance was compared with that of monomodal segmentation in PET and CT. In our experiments, mono-modal segmentation was performed by ECM-MS [26], i.e., by setting γ in (10) to zero, and removing simultaneously the fusion procedure described in Section III-C4.

Two illustrative results are shown in Fig. 4, from which we can find that co-segmentation outperformed independent segmentation in PET and CT in both two cases. Correspondingly, the quantitative results on the 21 subjects are summarized in Table I, from which we can find that co-segmentation led to the best performance in terms of all the five metrics, compared to independent segmentation produced by our previous mono-modal method (i.e., ECM-MS [26]). More specifically, it is worth noting that our previous ECM-MS obtained good segmentation performance in segmenting PET images, which thanks to two main reasons: 1) PET images are high in contrast; 2) blur and heterogeneity inherent in them were effectively handled in the framework of belief functions, with the help of the proposed spatial regularization and distance metric adaptation procedure. The extension of ECM-MS, i.e., the proposed co-clustering algorithm, further improved the segmentation in PET images, as anatomical information from CT images was included as the complementary knowledge for more accurate tumor delineation. On the other hand, the segmentation in CT images was greatly improved by the co-segmentation strategy. It is mainly because the boundary between target and background is invisible in CT images, which was effectively tackled by co-segmentation via incorporating functional information from PET images.

3). Comparison With the State-of-the-Art Methods:

In this group of experiments, our proposed co-clustering method was compared with four mono-modal methods: 1) 3D-LARW [12], 2) the original ECM [53], 3) SECM [24], and 4) FCM-SW [23]. In addition, it was also compared with two co-segmentation methods: 5) a statistical method using hidden Markov tree model (denoted as HMT) [44], and 6) M²EDN [47], a CNN-based method implemented via using PET and CT as multi-channel inputs. For the supervised M²EDN method, the segmentation results were quantified by 10-fold cross validation. To train a reliable network in each iteration, the training data was augmented on-the-fly via random rotations and small distortions.

The quantitative results of all competing methods are summarized in Table II, from which we can have at least two observations. First, the co-segmentation methods (i.e., HMT, M²EDN, and our proposed method) yield overall better performance than the mono-modal methods (i.e., 3D-LARW, ECM, SECM, and FCM-SW), which demonstrates that the combination of anatomical CT information and functional PET information is beneficial for more accurate tumor segmentation. Second, our proposed co-clustering method outperforms the other two co-segmentation methods (i.e., HMT and M²EDN). This is mainly due to: 1) Capitalizing on the specific context term (i.e., (9)) and fusion strategy defined in the framework of BFT, our proposed method provides an effective way to improve the delineation of target tumor. 2) The performance of HMT and M²EDN may be hampered, respectively, by the difficulty in conjointly defining the distributions of PET and CT voxels, and the limited number of training samples for the construction of a deep neural network.

To be more comprehensive, the visual examples obtained by all competing methods are also presented in Fig. 5 for comparison. The first column of Fig. 5 presents the axial slices of four different tumors, where the first row is a slice corresponds to a large tumor, the second row represents a relatively small tumor, the third row represents a smaller tumor, and the last row represents a heterogenous tumor. The second column to the last column of Fig. 5 compare the contours delineated by different methods (green line) with those delineated by clinicians (blue line). We can observe that the contours delineated by the proposed method (the last column) are more consistent with the reference contours.

Furthermore, it is worth mentioning that the segmentation results presented in this paper are also comparable to previous PET-CT co-segmentation studies. For example, our method yields DSC = 0.87 ± 0.04, which is comparable to Song et al. [33] (DSC = 0.81 ± 0.08), Ju et al. [36] (DSC = 0.84±0.06), Cui et al. [37] (DSC = 0.87±0.04), and Salah et al. [44] (DSC = 0.89 ± 0.13). It potentially implies that our proposed co-clustering method is a competitive alternative compared with existing co-segmentation methods.

1). Effectiveness of textural features:

Our proposed method adopts textural features in addition to intensity-based features for the characterization of image voxels, and then in an unsupervised way, informative features are automatically selected to further determine a discriminant feature transformation for the co-clustering of PET and CT voxels. To evaluate its effectiveness, we alternatively removed textural features extracted based on GLSZM and GLCM, and compared the corresponding segmentation results with those obtained by using all extracted features. The quantitative comparison in terms of DSC, SEN, and PPV is shown in Fig. 6. We can observe that including textural features could effectively improve the performance of our method, e.g., DSC and PPV were increased, respectively, from 0.83 ± 0.06 to 0.87 ± 0.04 and 0.83 ±0.15 to 0.88 ± 0.06 in this group of experiments.

2). Effectiveness of information fusion based on Dempster’s rule:

To further refine the clustering results, the Dempster’s combination rule (3) was adopted in our method for the fusion of complementary knowledge from PET and CT images (i.e. Section III-C4), as well as the fusion of information from neighboring voxels (i.e. Section III-D). To evaluate the effectiveness of such fusion strategy, we excluded it from the proposed method. The corresponding segmentation results are then summarized in Table III, from which we can find that the fusion procedure could help to refine the final segmentations.

3). Sensitivity to parameters:

As an illustration, we set γ to 0.001, and orderly chose a η and a λ from {0.001, …, 0.003, 0.07, 0.01, 0.03, …, 0.07, 0.1, 0.2} and {1, …, 10}, respectively. Then, the proposed method was applied to segment a relatively large tumor (volume of 135.80 mL) and a relatively small tumor (volume of 7.60 mL). The obtained DSCs are summarized in Fig. 7. We can observe that the obtained DSCs as a function of λ is relatively stable in a wide range. On the other hand, for the relatively large tumor, our method had better performance with small η; while, on the contrary, large η is better for the relatively small tumor.