Semi-supervised discriminative classification with application to tumorous tissues segmentation of MR brain images

Abstract

Due to the large data size of 3D MR brain images and the blurry boundary of the pathological tissues, tumor segmentation work is difficult. This paper introduces a discriminative classification algorithm for semi-automated segmentation of brain tumorous tissues. The classifier uses interactive hints to obtain models to classify normal and tumor tissues. A non-parametric Bayesian Gaussian random field in the semi-supervised mode is implemented. Our approach uses both labeled data and a subset of unlabeled data sampling from 2D/3D images for training the model. Fast algorithm is also developed. Experiments show that our approach produces satisfactory segmentation results comparing to the manually labeled results by experts.

Appendices

Appendix A: Explanation of EBM

For the semi-supervised problem, we set labels of the unlabeled data to zeros initially. Thus, if t _i = 0, the probability P(t _i  = 0|y) ≡ λ. The factor λ makes the function P(t _i |y _i ) with respect to t _i be a probability, which means P(t _i =  1|y _i ) + P(t _i = −1|y _i ) + P(t _i = 0|y _i ) ≡ 1. As Fig. 7 shows, this model can be considered as a degenerated ordered category model (OCM) [71], where the variance of the probability P(t _i  = 0|y _i ) is infinite. We define the margin as the range where P(t _i  = 0|y) is larger than P(t _i  = 1|y) and P(t _i  = −1|y). In the margin the difference between P(t _i  = 1|y) and P(t _i  = −1|y) is smaller than the difference outside the margin. Therefore, the margin of EBM represents the more uncertain labels. The parameter λ controls this margin.

Fig. 7

Illustration of extended Bernoulli model

Moreover, Fig. 7b–d show the relationship between the prior and the posterior probability of latent variable. In GP and GRF, we can assume that each latent variable is also conditionally Gaussian:

$$ P(y_i\vert{{\mathbf{y}}}_{N-\{i\}}, {{\mathbf{X}}}_N) = N(\mu_i, \sigma_i){\mathop=\limits^\Delta} P(y_i) $$

(21)

where μ _i and σ _i are related to the input points and labels (see the graphical model in Fig. 1b).

As Fig. 7b, c shows, for t _i  = 1 and t _i  = −1, the mean and the variance of posterior P(y _i |t _i  = 1) and P(y _i |t _i  = −1) are related to the likelihood P(t _i |y _i ) and the prior P(y _i ). If μ _i is near zero, the posterior of latent variable y _i is affected by the label t _i . It will have a positive estimated y _i when the label is 1, and negative y _i when the label is −1.

If the label is t _i  = 0, due to the Bayesian formulation, we have P(t _i |y _i )P(y _i ) = P(y _i |t _i )P(t _i ). The probabilities P(t _i ) and P(t _i |y _i ) are both constants for t _i = 0, so the posterior probability P(y _i |t _i  = 0) only depends on the prior P(y _i ). If μ _i is still near zero, we will get a zero estimated y _i by maximizing the posterior probability. This is why we choose a graph regularization-based prior. As mentioned earlier, each covariance between two points of this prior is related to all the training data. Thus, if there are a small amount of labeled data in the training set, μ _i of an unlabeled x _i will be affected by the labeled data more than the one choosing the traditional prior. Then, μ _i is non-zero, which will lead to a non-zero estimated y _i (shown in Fig. 7 c).

Furthermore, by comparing Fig. 7c and d, we can see that the margins do not affect the estimation of the latent variable. The estimated y _i remains the same in spite of different margin models being imposed on the process y → t. However, any point whose latent variable y _i that falls inside the margin will be labeled zero, which makes it remains unlabeled. This kind of points does not contribute to the prediction function (see in the prediction phase). Therefore, the classification boundary will be changed.

Appendix B: Derivation of hyper-parameter estimation

The hyper-parameters are the standard deviations Θ = {σ _c , σ _p } of the exponential weight, which is shown in (19). We take σ _c as an example. The hyper-parameter is estimated by gradient search, which minimizes the negative logarithmic likelihood (10). According to Eqs. 5 and 8 by using the fact ∇Ψ = 0, and taking derivatives on both sides, we have

$$ \frac{\partial{{\mathbf{y}}}_N}{\partial \sigma_c} = ({{{\mathbf{I}}} + {{\mathbf{K}}}_N{\varvec{\Uppi}}_N})^{- 1}\frac{\partial{{\mathbf{K}}}_N} {\partial \sigma_c}(- {\varvec{\upalpha}}_N) $$

(22)

And according to Eq. 4, we have

$$ \begin{aligned} \frac{\partial {{\mathbf{K}}}_N^{- 1}}{\partial \sigma_c} &= - {{\mathbf{D}}}^{- \frac{1}{2}}\frac{\partial{{\mathbf{W}}}}{\partial \sigma_c} {{\mathbf{D}}}^{- \frac{1}{2}} \\ &\quad + {{\mathbf{D}}}^{- 1}\left(\frac{\partial \sqrt{{\mathbf{D}}}} {\partial \sigma_c}{{\mathbf{W}}}\sqrt{{\mathbf{D}}} + \sqrt{{\mathbf{D}}}{{\mathbf{W}}} \frac{\partial \sqrt{{\mathbf{D}}} }{\partial \sigma_c}\right){{\mathbf{D}}}^{- 1} \end{aligned} $$

(23)

where

$$ \begin{aligned} \frac{\partial{{\mathbf{W}}}_{ij}}{\partial \sigma_c} &= {{\mathbf{W}}}_{ij} {\left\|{{{\mathbf{x}}}_c(i) - {{\mathbf{x}}}_c(j)} \right\|^2} / {\sigma_c^3}\\ \frac{\partial \sqrt{{{\mathbf{D}}}_{ii}}}{\partial \sigma_c}& = \frac{1}{2\sqrt{{{\mathbf{ D}}}_{ii}}}\sum_k {{{\mathbf{W}}}_{ik} {\left\|{{{\mathbf{x}}}_c(i) - {{\mathbf{x}}}_c(k)} \right\|^2}} / {\sigma_c^3}\\ \end{aligned}$$

(24)

Furthermore, differences of K _N and ${{\mathbf{K}}}_N{\varvec{\Uppi}}_N$ are given by

$$ \frac{\partial {{\mathbf{K}}}_N}{\partial \sigma_c} = - {{\mathbf{K}}}_N\frac{\partial{{\mathbf{K}}}_N^{- 1}}{\partial \sigma_c}{{\mathbf{K}}}_N $$

(25)

and

$$ \frac{\partial{{\mathbf{K}}}_N{\varvec{\Uppi}}_N}{\partial \sigma_c} = \frac{\partial{{\mathbf{K}}}_N}{\partial \sigma_c} {\varvec{\Uppi}}_N + {{\mathbf{K}}}_N\frac{\partial {\varvec{\Uppi}}_N}{\partial \sigma_c} $$

(26)

where

$$ \begin{aligned} \frac{\partial {\varvec{\Uppi}}_N}{\partial \sigma_c}& = \frac{\partial {\varvec{\Uppi}}_N} {\partial {{\mathbf{y}}}_N}\frac{\partial{{\mathbf{y}}}_N}{\partial \sigma_c} \\ &= {\rm diag}\left(\frac{t_i^3\exp (t_i y_i )\left( {1 - \exp (t_i y_i )} \right)}{\left( {1 + \exp (t_i y_i )} \right)^3} \right) {\rm diag}{\left(\frac{\partial{{\mathbf{y}}}_N}{\partial \sigma_c }\right)} \end{aligned} $$

(27)

Therefore, the derivative of objective function can be given by (11).

Appendix C: Derivation of prediction function

For predicting new test points, we first distinguish between different sizes of covariance matrix K with a subscript, such that: (1) K _N is the covariance matrix of the N input training data, (2) K _N+1 is (N + 1) × (N + 1) covariance matrix of the vector (y ^T _N ,y _N+1)^T, where y _N+1 is the latent variable of a new given point in the test set. However, it is difficult to compute K _N+1 by K _N in an explicit expression, since K _N+1 itself depends on all the training data and each new point.³ We make use of ${\mathbf{k}}_i = W_{N+1,i} = \exp (-{\|{\mathbf{x}}_{N+1}-{\mathbf{x}}_{i} \|^2} / {2\sigma}^2 )$ as the covariance between a new given point and the ith training point. Therefore, K _N+1 can be given by

$$ {{\mathbf{K}}}_{N+1}= \left[\begin{array}{*{20}l} {{\mathbf{K}}}_N & \nu{{\mathbf{k}}}\\ \nu{{\mathbf{k}}}^{\rm T} & k_*\\\end{array}\right] $$

(28)

where ν is a scale factor to make k compatible with K _N . Note that the covariance matrix K _N depends on the global distance information, while the covariances of new point and the training points are only depend on the local distance information.

Then, K ⁻¹ _N+1 is given by

$$ {{\mathbf{K}}}_{N+1}^{-1}= \left[\begin{array}{*{20}l} {{\mathbf{M}}} & {{\mathbf{m}}}\\ {{\mathbf{m}}}^{\rm T} & \mu\\ \end{array}\right] $$

(29)

By using the partitioned inverse equations [61], we have

$$ \begin{aligned} \mu &= (k_* - \nu^2 {{\mathbf{k}}}^{\rm T} {{\mathbf{K}}}_N^{-1} {{\mathbf{k}}})^{-1}\\ {{\mathbf{m}}} &= -\mu \nu {{\mathbf{K}}}_N^{-1} {{\mathbf{k}}} \\ {{\mathbf{M}}} &= {{\mathbf{K}}}_N^{-1} + \frac{1}{\mu} {{\mathbf{m}}} {{\mathbf{m}}}^{\rm T}\\ \end{aligned} $$

(30)

Since the optimal ${\hat{{\mathbf{y}}}}_N$ has been estimated, we only need to minimize (15) with respect to y _N+1: μ y _N+1 + m ^T y _N = 0, which leads to

$$ {\hat{y}}_{N + 1} = {{\mathbf{k}}}^{\rm T}{{{\mathbf{K}}}_N}^{- 1}{\hat{{\mathbf{y}}}}_N = {{\mathbf{k}}}^{\rm T}({{\mathbf{I - S}}}){\hat{{\mathbf{y}}}}_N. $$

(31)

Here the scale factor ν can be omitted.

Appendix D: Explanation of the decomposition

The authors in [42] proved that the weight matrix W can be approximately decomposed as

$$ {\hat{{\mathbf{W}}}} = {{\mathbf{V}}} {\varvec{\Uplambda}}_{{\mathbf{S}}} {{\mathbf{V}}}^{\rm T} $$

(32)

where

$$ {{\mathbf{V}}} = \left[\begin{array}{*{20}l} {{\mathbf{A}}} \\ {{\mathbf{B}}} \\\end{array} \right] {{\mathbf{A}}}^{-1/2} {{\mathbf{U}}}_{{\mathbf{S}}} {\varvec{\Uplambda}}_{{\mathbf{S}}}^{-1/2} $$

(33)

A is the sub-block of weight matrix generated by the sampling points, B is the weight matrix between the sampling points and the rest. We assume that W is denote by ${{\mathbf{W}}} = \left[ \begin{array}{*{20}l}{{\mathbf{A}}} & {{\mathbf{B}}} \\ {{\mathbf{B}}}^{\rm T} & {{\mathbf{C}}}\\\end{array} \right],$ where C is the weight matrix of the rest points. The matrix U _S and ${\varvec{\Uplambda}}_{{\mathbf{S}}}$ is the eigenvectors and values of the matrix ${{\mathbf{S}}} = {{\mathbf{A}}} + {{\mathbf{A}}}^{-1/2} {{\mathbf{B}}} {{\mathbf{B}}}^{\rm T} {{\mathbf{A}}}^{-1/2} = {{\mathbf{U}}}_{{\mathbf{S}}} {\varvec{\Uplambda}}_{{\mathbf{S}}} {{\mathbf{U}}}_ {{\mathbf{S}}}^{\rm T},$ where A ^−1/2 denotes the symmetric positive definite square root⁴ of A. Then we can find that the constraint V V ^T = I is satisfied automatically and the approximated weight matrix is given by

$$ {\hat{{\mathbf{W}}}} = \left[\begin{array}{*{20}l} {{\mathbf{A}}} & {{\mathbf{B}}} \\ {{\mathbf{B}}}^{\rm T} & {{\mathbf{B^{\rm T} A^{-1}} \mathbf{B}}} \\ \end{array}\right] $$

(34)

The difference of C and B ^T A ⁻¹ B is just the Schur complement. In addition, for the purpose of our accelerated algorithm, we need to approximate the matrix ${\hat{{\mathbf{D}}}}^{- \frac{1}{2}} {\hat{{\mathbf{W}}}} {\hat{{\mathbf{D}}}}^{- \frac{1}{2}}.$ This has been exploited by [42], which replace the matrix A and B as

$$ \begin{aligned} {{\mathbf{A}}}_{ij} &\leftarrow\frac{{{\mathbf{A}}}_{ij}}{\sqrt{{\hat{{\mathbf{ d}}}_i} {\hat{{\mathbf{d}}}}_j}} \quad(i,j = 1,\ldots,n)\\ {{\mathbf{B}}}_{ij}& \leftarrow \frac{{{\mathbf{B}}}_{ij}}{\sqrt{\hat{{\mathbf{d}}}}_i {\hat{{\mathbf{d}}}}_{j+m}}\quad (i = 1,\ldots,n;j=1,\ldots,m) \end{aligned} $$

(35)

The vector ${\hat{{\mathbf{d}}}}$ can be evaluated by

$$ {\hat{{\mathbf{d}}}} = {\hat{{\mathbf{W}}}} {{\mathbf{1}}} = \left[\begin{array}{*{20}l} {{\mathbf{A}}}{{\mathbf{1}}}_M + {{\mathbf{B}}}{{\mathbf{1}}}_{N-M} \\ {{\mathbf{B}}}^T {{\mathbf{1}}}_M + {{\mathbf{B^T A^{-1}} \mathbf{B}}} {{\mathbf{1}}}_{N-M} \\ \end{array}\right] $$

(36)

where 1 is the column vector of ones. Thus, we can rewrite the decomposition of ${\hat{{\mathbf{D}}}}^{- \frac{1}{2}} {\hat{{\mathbf{W}}}} {\hat{{\mathbf{D}}}}^{- \frac{1}{2}}$ as ${{\mathbf{U}}}{\varvec{\Uplambda}}{{\mathbf{U}}}^{\rm T},$ where ${\varvec{\Uplambda}}$ is an M × M matrix.