Model-based computed tomography image estimation: partitioning approach

ABSTRACT

There is a growing interest to get a fully MR based radiotherapy. The most important development needed is to obtain improved bone tissue estimation. The existing model-based methods perform poorly on bone tissues. This paper was aimed at obtaining improved bone tissue estimation. Skew-Gaussian mixture model and Gaussian mixture model were proposed to investigate CT image estimation from MR images by partitioning the data into two major tissue types. The performance of the proposed models was evaluated using the leave-one-out cross-validation method on real data. In comparison with the existing model-based approaches, the model-based partitioning approach outperformed in bone tissue estimation, especially in dense bone tissue estimation.

Keywords:

• Computed tomography
• magnetic resonance imaging
• CT image estimation
• skew-Gaussian mixture model
• Gaussian mixture model

1. Introduction

Magnetic resonance (MR) imaging and computed tomography (CT) are the most widely used diagnostic imaging technologies in medicine. They are used to obtain more detailed cross-sectional images of the human body. CT uses ionizing radiation to record a pattern of radiodensities to obtain cross-sectional images. The ionizing radiations are attenuated as they pass through the tissues of patients. The amount of attenuations depends on the tissue types. The differences in attenuation between adjacent tissues create contrast on CT images. Tissues with higher (or lower) attenuation appear brighter (or darker) on grayscale CT images. As a result, air, soft, and bone tissues appear as darkest, darker, and white on grayscale CT images. Therefore, the CT image is excellent for identifying and assessing the structures of bone tissues. On the other hand, exposing a patient to ionizing radiation in CT imaging may have a risk for radiation-related cancer.

MR imaging is remarkably different from CT. It does not depend on ionizing radiations. MR imaging relies on the absorption and emission of radio waves from tissue protons exposed to a strong magnetic field. Thus, MR imaging is safer than CT imaging. The relative MR signal intensity differences between adjacent anatomic structures determine tissue contrast on MR images. In comparison with CT images, MR images show much better soft tissue contrast and noticeably improves the delineation of tumors. However, MR images are poor in depicting bone tissues. The reason is that bone, air, and rapidly flowing blood appear black on grayscale MR images.

A new innovation MR-only based radiotherapy enhances tissue contouring and precision in soft tissue therapy setup. It also improves biological information at treatment planning and avoids registration errors, which are errors due to the transformation of different images of the same scene into one coordinate system, between CT and MR images. Moreover, MR-only based radiotherapy is a cost-effective approach as it reduces redundant imaging. However, co-registered CT and MR image complement each other due to better bone tissue imaging of CT [8,10,20]. MR images are not directly applicable for attenuation correction. Attenuation refers to the loss of signals due to absorption or scattering out of the signals in the body. CT images are vital for attenuation correction in positron emission tomography (PET) imaging. This is due to the direct relation between CT image intensities and PET attenuation coefficients. However, the CT scanner is not available in recently combined PET/MR imaging scanner. Therefore, MR-only based radiotherapy and the combined PET/MR imaging scanner can be successful if one can obtain a reliable CT image estimation. As a result, we need to develop a reliable CT image estimation method from MR images.

Huynh et al. [14] used a learning-based method to estimate CT image from the MR image. A patch of CT image is estimated directly from a given MR image patch using the structured random forest. The robustness of the estimation has been evaluated using a new ensemble model. Nie et al. [30] proposed a 3D deep learning-based method for patch-wise estimation of CT images from MR images. The neural network generates structured output and it preserves the neighborhood information in the estimated CT image. Arabi et al. [1] suggested a two-step atlas-based algorithm to estimate CT image from MR image sequences. The estimation is mainly concerned with pinpointing of bone tissues.

Johansson et al. [16] used a Gaussian mixture model (GMM) to obtain CT substitute from MR images without taking spatial dependence between neighboring voxels into account. Johansson et al. [17] investigated the uncertainty associated to the voxel-wise estimation of CT images. By considering spatial dependence between neighboring voxels, Kuljus et al. [22] extended the work of Johansson et al. [16] by using hidden Markov model (HMM) and Markov random field (MRF) model. Kuljus et al. [22] compared the estimation quality of GMM, HMM, and MRF. In terms of mean absolute error, HMM outperformed the other models and it was computationally robust than MRF. However, it had a weaker estimation quality on dense bone tissues. Even though MRF had superior performance on bone tissue estimation, it was computationally expensive.

The main aim of this article is to further investigate the voxel-wise estimation of CT images from MR images by partitioning the data into non-bone and bone tissues. According to Johansson et al. [17] and Kuljus et al. [22], the estimation of CT image was poor on air and bone tissues. It is this result that motivated the partitioning of the data into non-bone and bone tissues in order to further explore the estimation of CT images. Even though there is no clear-cut CT image intensity boundary between these tissue types, Waterstram-Rich and Gilmore [33] and Washington and Leaver [34] provide informative threshold delimiting these tissues. Waterstram-Rich and Gilmore [33] used 150 Hounsfield units (HU) as a lower limit for bone tissues. On the other hand, Washington and Leaver [34] utilized 200 HU as an approximate delimiting value of the tissues.

The partitioning of the data may introduce skewness. Consequently, there is a need to relax the normality assumption used in [16,17,22]. Azzalini [4] proposed a univariate skew-normal model that relaxed the normality assumption by incorporating a skewness parameter in the distributional assumption. Azzalini and Dalla Valle [5] extended the univariate skew-normal to a multivariate skew-normal. A multivariate skew-normal is a tractable extension of a multivariate normal distribution with an extra parameter to regulate skewness. Lin et al. [27] introduced a univariate skew-normal mixture model in order to deal with population heterogeneity and skewness. Lin [25] extended the univariate skew-normal mixture model to a multivariate skew-normal mixture model which is an alternative to the most widely used multivariate Gaussian mixture model.

Kahrari et al. [19] developed a multivariate skew-normal-Cauchy distribution and represented it as a shape mixture of the multivariate skew-normal distribution. Kahrari et al. [18] modified the multivariate skew-normal-Cauchy distribution and the modified version becomes a shape mixture of a special case of the fundamental skew-normal distribution developed by Arellano-Valle and Genton [3] with a univariate half-normal mixing distribution. The class of scale mixtures of skew-normal-Cauchy distributions has been represented as a shape mixture of the class of scale mixtures of skew-normal distributions with a univariate half-normal mixing distribution [18]. Jamalizadeh and Lin [15] presented the scale-shape mixtures of skew-normal distributions for modeling asymmetric data. Arellano-Valle et al. [2] established a flexible class of multivariate distributions obtained by both scale and shape mixtures of multivariate skew-normal distributions. Based on the skew-t-normal distribution [13], Tamandi et al. [32] introduced the shape mixtures of the skew-t-normal distributions that contain one additional shape parameter to regulate skewness and kurtosis. The shape mixtures of the skew-t-normal distributions are a flexible extension of the skew-t-normal distribution. Lin et al. [26] developed a multivariate extension of the skew-t-normal distribution that is obtained as a scale mixture of the multivariate skew-normal distribution introduced by Azzalini and Dalla Valle [5]. Based on the multivariate skew-t-normal distribution, Lin et al. [26] introduced a robust probabilistic mixture model which is composed of a weighted sum of a finite number of different multivariate skew-t-normal densities. The flexible mixture model based on the multivariate skew-t-normal distribution includes mixtures of normal, t and skew-normal distributions as special cases. Cabral et al. [9] proposed mixture models which consist of members of skew-normal independent distributions (the skew-normal, the skew-t, the skew-slash and the skew-contaminated normal) and the mixture models are developed using the multivariate skew-normal distribution in [5].

The most common approach to estimate the parameters of these skew-mixture models is the EM algorithm [11]. However, its M-step for the recent skew-mixture models is computationally intractable. Alternatively, we use an EM-type algorithm to estimate these skew-mixture models. That is, we make further assumptions at the M-step of the EM algorithm. For instance, expectation conditional maximization (ECM) algorithm [29], which replaces the M-step of EM with several computationally simple conditional maximization steps, can be enough to estimate the parameters of some mixture models. In mixtures of multivariate skew-t-normal distributions and shape mixtures of skew-t-normal distributions, expectation conditional maximization either (ECME) algorithm [28] was exploited to estimate its parameters by replacing some conditional maximization-steps of ECM with the conditional maximization likelihood step that maximizes the correspondingly constrained actual-likelihood function. Cabral et al. [9] developed an EM-type algorithm that removed some obstacles (for instance, Monte Carlo integration) during the parameter estimation process in mixture models of skew-normal independent distributions. In this article, we used a mixture model based on the multivariate skew-normal distribution in [5] and developed EM-algorithm for its parameter estimation. That is, further assumption is not required at the M-step of EM algorithm to estimate the parameters of the skew-Gaussian mixture model.

In this work, we use skew-Gaussian mixture model (SGMM) and Gaussian mixture model (GMM) to further explore the estimation of CT images by partitioning the data into two major tissue types: non-bone and bone tissues. Non-bone tissues consist of subclasses such as white matter, blood, water, fat, gray matter, air, etc. while bone tissues consist of subclasses such as cranium, mandible, frontal bone, nasal bone, orbital bones, cortical bone, cancellous bone, etc. These facts motivated us to apply mixture models on these tissue types. SGMM involves a weighted sum of the joint skew-normal distributions of a CT image intensity and its corresponding intensities of MR images. The number of skew-normal distributions in the mixture depends on the number of underlying tissue types or clusters. Latent variables that represent the underlying tissue types are utilized during the parameter estimation process of the model through incomplete data assumption in EM-algorithm framework [11]. Voxel-wise point estimator of CT image was obtained as a weighted sum of the conditional expected value of a CT image intensity given its corresponding intensities of MR images and the underlying tissue type. The probability that an underlying tissue type is determined based on the intensities of MR images was used as a weight of the conditional expected value.

In summary, this study is concerned with comparing the CT image estimation performance of the partitioning approach with HMM, MRF, and GMM on bone tissues. The models HMM, MRF, and GMM are trained on the full data (data that are not partitioned into non-bone and bone tissues). We are also interested to compare the predictive quality of the partitioning approaches, SGMM and GMM${}^{\ast }$ (GMM applied to each partition), on bone tissues.

This article is organized as follows. The second section describes data acquisition and demonstrates statistical methods. The third section presents the results obtained and the final section discusses the implication of the results.

2. Statistical methodology

In this section, we describe the data, formulate SGMM, and develop the parameter estimation method. We also demonstrate the CT image estimation method and present evaluation method of the estimation. For the remaining models (GMM, HMM and MRF), we refer to Kuljus et al. [22] and the references therein.

2.1. Data acquisition

CT and MR images were obtained from the head of five patients. Four MR images were acquired from each patient using two dual echo ultrashort echo-time sequences with flip angles of 10 degrees and 30 degrees. The ultrashort echo-time sequences sampled a first echo (free induction decay) and a second echo (gradient echo) from the same excitation with an echo time of 0.07 and 3.76 ms. MR image of a patient was reconstructed to $192\times 192\times 192$ matrix. An entry in the matrix represents a signal intensity corresponding to a three-dimensional tissue (voxel) with size $1.33\times 1.33\times 1.33{\mathrm{mm}}^3$. One CT image of a patient was acquired using gradient echo Lightspeed with 2.5 mm slice thickness. The acquired CT image was reconstructed with an in-plane resolution of $0.78\times 0.78{\mathrm{mm}}^2$. One binary mask (an image with voxel value 1 (or 0) representing the region of interest (or the surrounding air)) was also developed to demarcate the head of a patient from its surrounding air. The main use of the binary mask is to exclude the surrounding air from the acquired CT and MR images. For each patient, the binary mask, the CT image, and the four MR images were co-registered and resampled to the same resolution (voxel-to-voxel correspondence and set to the same voxel dimension) using linear interpolation. For further technical details, we refer to Johansson et al. [16]. Voxel values of the CT image, the binary mask, and the four MR images were organized into six columns to obtain data for a patient. The organized data of each patient were column stacked and the surrounding air removed to obtain data for model fitting. Figure 1 shows a slice data for a given patient.

Figure 1. Binary mask (a), CT image (b) and MR images (c)–(f).

2.2. Data partitioning

This subsection describes data partitioning during model training. It also demonstrates how MR images of new patients are utilized during CT image prediction.

2.2.1. Data partition: model training

CT image intensity threshold was utilized to partition the data into two major tissue types. Using 150 HU CT image intensity as a limit, 50 HU (is selected to take the delimiting value provided by Washington and Leaver [34]) overlap was allowed during the parameter estimation process. The overlap of the major tissue types was motivated in order to minimize the effect of fuzzy boundary of the tissue types. Accordingly, CT image intensities in (− 1024 HU, 200 HU) and ( 100 HU, 3071 HU] were assumed to represent non-bone and bone tissues. The minimum number of voxels for non-bone and bone tissues are 6,214,160 and 1,292,068.

2.2.2. MR images partition: CT image estimation

We are interested to predict CT images from MR images of new patients. Since we only have MR images of the new patients, there is a need to partition the MR images of the new patients into non-bone and bone tissues. There is poor bone tissue information on MR images and following that we estimate CT images for the new patients using a model trained on the full data, which is the data not partitioned into non-bone and bone tissues. The CT image intensity threshold and the estimated CT images are used to obtain MR data corresponding to the two major tissue types (non-bone and bone tissues). The trained models on the major tissue types are utilized to obtain the desired CT images of the new patients.

2.3. Statistical model: mixture of multivariate skew-normal model

Let $Y_{i1}$ and ${\mathbf{Y}}_{i2}=(Y_{i2},Y_{i3},\dots ,Y_{id}{)}^{\prime }$ represent voxel i of CT image and its corresponding MR images. In our real data, $d=5$. A d-dimensional random vector ${\mathbf{Y}}_i=(Y_{i1},{\mathbf{Y}}_{i2}^{\prime }{)}^{\prime }$ is assumed to follow a multivariate skew-normal distribution $\mathcal{S}\mathcal{N}({\mathbf{y}}_i|\mathit{\eta },\mathbf{\Sigma },\mathit{\lambda })$ with a d-dimensional location parameter vector $\mathit{\eta }$, a $d\times d$-dimensional positive definite dispersion matrix $\mathbf{\Sigma }$, and a d-dimensional skewness parameter vector $\mathit{\lambda }$. Its density can be given by (1) $f\left({\mathbf{y}}_i|\mathit{\eta },\mathbf{\Sigma },\mathit{\lambda }\right)=2\mathcal{N}\left({\mathbf{y}}_i|\mathit{\eta },\mathbf{\Sigma }\right)\mathrm{\Phi }\left({\mathit{\lambda }}^{\prime }{\mathbf{\Sigma }}^{-1/2}\left({\mathbf{y}}_i-\mathit{\eta }\right)\right),$(1) where ${\mathbf{\Sigma }}^{-1/2}{\mathbf{\Sigma }}^{-1/2}={\mathbf{\Sigma }}^{-1}$, $\mathrm{\Phi }(\cdot )$ is a univariate standard normal distribution function, $i=1,2,3,\dots ,n$, and n is the number of voxels. According to Lachos et al. [24], the stochastic representation of ${\mathbf{Y}}_i$ may be given by (2) ${\mathbf{Y}}_i=\mathit{\eta }+{\mathbf{\Sigma }}^{1/2}\mathit{\delta }{U}_i+{\mathbf{\Sigma }}^{1/2}{\left(\mathbf{1}-\mathit{\delta }{\mathit{\delta }}^{\prime }\right)}^{1/2}{\mathbf{V}}_i,$(2) where (3) $U_i\sim \mathcal{H}\mathcal{N}\left(u_i|0,1,\left(0,\mathrm{\infty }\right)\right),{\mathbf{V}}_i\sim \mathcal{N}\left({\mathbf{v}}_i|\mathbf{0},\mathbf{1}\right),\mathrm{and}\mathit{\delta }=\frac{\mathit{\lambda }}{\sqrt{1+{\mathit{\lambda }}^{\prime }\mathit{\lambda }}}.$(3) In Equation (2(2) ${\mathbf{Y}}_i=\mathit{\eta }+{\mathbf{\Sigma }}^{1/2}\mathit{\delta }{U}_i+{\mathbf{\Sigma }}^{1/2}{\left(\mathbf{1}-\mathit{\delta }{\mathit{\delta }}^{\prime }\right)}^{1/2}{\mathbf{V}}_i,$(2) ), ${\mathbf{V}}_i$ and $U_i$ are assumed to be independent. The notations $\mathbf{0}$, $\mathbf{1}$, and $\mathcal{H}\mathcal{N}(u_i|0,1,(0,\mathrm{\infty }\left)\right)$ in Equation (3(3) $U_i\sim \mathcal{H}\mathcal{N}\left(u_i|0,1,\left(0,\mathrm{\infty }\right)\right),{\mathbf{V}}_i\sim \mathcal{N}\left({\mathbf{v}}_i|\mathbf{0},\mathbf{1}\right),\mathrm{and}\mathit{\delta }=\frac{\mathit{\lambda }}{\sqrt{1+{\mathit{\lambda }}^{\prime }\mathit{\lambda }}}.$(3) ) represent d-dimensional zero vector, $d\times d$-dimensional identity matrix, and half-normal distribution, respectively. Using Equation (2(2) ${\mathbf{Y}}_i=\mathit{\eta }+{\mathbf{\Sigma }}^{1/2}\mathit{\delta }{U}_i+{\mathbf{\Sigma }}^{1/2}{\left(\mathbf{1}-\mathit{\delta }{\mathit{\delta }}^{\prime }\right)}^{1/2}{\mathbf{V}}_i,$(2) ), a hierarchical model can be given by (4) $\begin{array}{rl}{\mathbf{Y}}_i|U_i& =u_i\sim \mathcal{N}\left({\mathbf{y}}_i|\mathit{\eta }+u_i\mathit{\xi },\mathbf{\Omega }\right),\\ & U_i\sim \mathcal{H}\mathcal{N}\left(u_i|0,1,\left(0,\mathrm{\infty }\right)\right),\end{array}$(4) where $\mathit{\delta }=\frac{\mathit{\lambda }}{\sqrt{1+{\mathit{\lambda }}^{\prime }\mathit{\lambda }}},\mathit{\xi }={\mathbf{\Sigma }}^{1/2}\mathit{\delta },\mathrm{and}\mathbf{\Omega }={\mathbf{\Sigma }}^{1/2}\left(\mathbf{1}-\mathit{\delta }{\mathit{\delta }}^{\prime }\right){\mathbf{\Sigma }}^{1/2}=\mathbf{\Sigma }-\mathit{\xi }{\mathit{\xi }}^{\prime }.$ Let $Z_i$ be a categorical random variable representing the underlying tissue types at voxel i. Define an indicator variable: $\begin{array}{r}{Z}_{ik}=1_{\left(Z_i=s\right)}=\left\{\begin{array}{ll}1& \mathrm{if}k=s,\\ 0& \mathrm{otherwise},\end{array}\right.\end{array}$ where $k=1,2,\dots ,K$. The definition implies that $P(Z_{ik}=1)=P(Z_i=k)$. Let $P(Z_i=k)$= ${\pi }_k$, which represents the weight that the $i\mathrm{th}$ observation belongs to a tissue class k. To incorporate tissue heterogeneity into the statistical modeling, ${\mathbf{Y}}_i|Z_{ik}=1$ is assumed to follow a multivariate skew-normal distribution $\mathcal{S}\mathcal{N}({\mathbf{y}}_i|{\mathit{\eta }}_k,{\mathbf{\Sigma }}_k,{\mathit{\lambda }}_k)$. This means that ${\mathbf{Y}}_i$ follows a mixture of multivariate skew-normal distributions. Its density may be given by $f\left({\mathbf{y}}_i|\mathbf{\Psi }\right)=\sum _{k=1}^K{\pi }_k\mathcal{S}\mathcal{N}\left({\mathbf{y}}_i|{\mathit{\eta }}_k,{\mathbf{\Sigma }}_k,{\mathit{\lambda }}_k\right),$ where $\begin{array}{rl}& {\pi }_k\ge 0,\sum _{k=1}^K{\pi }_k=1,{\mathit{\psi }}_k=\left\{{\pi }_k,{\mathit{\eta }}_k,{\mathbf{\Sigma }}_k,{\mathit{\lambda }}_k\right\},\mathbf{\Psi }=\left\{{\mathit{\psi }}_1,{\mathit{\psi }}_2,\dots ,{\mathit{\psi }}_K\right\},\\ & k=1,2,\dots ,K.\end{array}$ The unknown parameters in $\mathbf{\Psi }$ are estimated from the independent observations ${\mathbf{y}}_i$.

2.4. Parameter estimation method

The log-likelihood function of the data $\mathbf{y}=({\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_n{)}^{\prime }$ can be given by $\log f\left(\mathbf{y}|\mathbf{\Psi }\right)=\sum _{i=1}^n\log \left\{\sum _{k=1}^K{\pi }_k\mathcal{S}\mathcal{N}\left({\mathbf{y}}_i|{\mathit{\eta }}_k,{\mathbf{\Sigma }}_k,{\mathit{\lambda }}_k\right)\right\}.$ In general, there is no explicit analytical solution for $\arg \underset{\mathbf{\Psi }}{max}\log f\left(\mathbf{y}|\mathbf{\Psi }\right)$. However, iterative maximizing procedure under the idea of incomplete data via EM-algorithm can be used to obtain an optimal estimate of the parameters. Let ${\mathbf{Z}}_i$ be a K-dimensional column vector of $Z_{ik}$. Its realization is a K-dimensional vector consisting 1 at only one location and 0 at the remaining locations. The latent random vector ${\mathbf{Z}}_i$ follows a multinomial distribution with one trial and $P(Z_{ik}=1)={\pi }_k$. Using the indicator variable $Z_{ik}$, the hierarchical model (4(4) $\begin{array}{rl}{\mathbf{Y}}_i|U_i& =u_i\sim \mathcal{N}\left({\mathbf{y}}_i|\mathit{\eta }+u_i\mathit{\xi },\mathbf{\Omega }\right),\\ & U_i\sim \mathcal{H}\mathcal{N}\left(u_i|0,1,\left(0,\mathrm{\infty }\right)\right),\end{array}$(4) ) can be extended to $\begin{array}{rl}{\mathbf{Y}}_i|U_i=u_i,Z_{ik}& =1\sim \mathcal{N}\left({\mathbf{y}}_i|{\mathit{\eta }}_k+u_i{\mathit{\xi }}_k,{\mathbf{\Omega }}_k\right),\\ U_i|Z_{ik}& =1\sim \mathcal{H}\mathcal{N}\left(u_i|0,1,\left(0,\mathrm{\infty }\right)\right),\end{array}$ where ${\mathit{\delta }}_k=\frac{{\mathit{\lambda }}_k}{\sqrt{1+{\mathit{\lambda }}_k^{\prime }{\mathit{\lambda }}_k}},{\mathit{\xi }}_k={\mathbf{\Sigma }}_k^{1/2}{\mathit{\delta }}_k,\mathrm{and}{\mathbf{\Omega }}_k={\mathbf{\Sigma }}_k-{\mathit{\xi }}_k{\mathit{\xi }}_k^{\prime }.$ The observed data $\mathbf{y}$ is assumed to be incomplete data. It is augmented with the latent matrix $\mathbf{z}=({\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_n{)}^{\prime }$ and the latent vector $\mathbf{u}=(u_1,u_2,\dots ,u_n{)}^{\prime }$ to form a complete dataset $(\mathbf{y},\mathbf{u},\mathbf{z})$ in EM-algorithm framework. Assuming that $({\mathbf{Y}}_i,U_i,{\mathbf{Z}}_i)$ is independent of $({\mathbf{Y}}_j,U_j,{\mathbf{Z}}_j)$ for every $i\ne j$, the complete-data log-likelihood is given by (5) $\begin{array}{rl}\log f\left(\mathbf{y},\mathbf{u},\mathbf{z}|\mathbf{\Theta }\right)& =\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\left\{-\frac{1}{2}{\left({\mathbf{y}}_i-{\mathit{\eta }}_k-{\mathit{\xi }}_k{u}_i\right)}^{\prime }{\mathbf{\Omega }}_k^{-1}\left({\mathbf{y}}_i-{\mathit{\eta }}_k-{\mathit{\xi }}_k{u}_i\right)\right\}\\ & +\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\left\{\log {\pi }_k-\frac{1}{2}\log |{\mathbf{\Omega }}_k|+const\right\},\end{array}$(5) where ${\mathit{\theta }}_k=\left\{{\pi }_k,{\mathit{\eta }}_k,{\mathit{\xi }}_k,{\mathbf{\Omega }}_k\right\},\mathbf{\Theta }=\left\{{\mathit{\theta }}_1,{\mathit{\theta }}_2,\dots ,{\mathit{\theta }}_K\right\},k=1,2,\dots ,K,$ and const is a constant function of parameters.

The E-step of EM-algorithm involves computing the expected value of the complete-data log-likelihood given $\mathbf{Y}$ and the current estimate ${\mathbf{\Theta }}^{old}$ of $\mathbf{\Theta }$. It may be given by the Q-function: (6) $\mathcal{Q}\left(\mathbf{\Theta },{\mathbf{\Theta }}^{\mathrm{old}}\right)=E\left[\log f\left(\mathbf{Y},\mathbf{U},\mathbf{Z}|\mathbf{\Theta }\right)|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right].$(6) Using Equation (5(5) $\begin{array}{rl}\log f\left(\mathbf{y},\mathbf{u},\mathbf{z}|\mathbf{\Theta }\right)& =\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\left\{-\frac{1}{2}{\left({\mathbf{y}}_i-{\mathit{\eta }}_k-{\mathit{\xi }}_k{u}_i\right)}^{\prime }{\mathbf{\Omega }}_k^{-1}\left({\mathbf{y}}_i-{\mathit{\eta }}_k-{\mathit{\xi }}_k{u}_i\right)\right\}\\ & +\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\left\{\log {\pi }_k-\frac{1}{2}\log |{\mathbf{\Omega }}_k|+const\right\},\end{array}$(5) ), the expected value in Equation (6(6) $\mathcal{Q}\left(\mathbf{\Theta },{\mathbf{\Theta }}^{\mathrm{old}}\right)=E\left[\log f\left(\mathbf{Y},\mathbf{U},\mathbf{Z}|\mathbf{\Theta }\right)|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right].$(6) ) involves computing $E\left[Z_{ik}|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right],E\left[Z_{ik}{U}_i|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right],\mathrm{and}E\left[Z_{ik}{U}_i^2|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right].$ The expected value $E[Z_{ik}|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}]$ can be given by $\begin{array}{rl}E\left[Z_{ik}|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right]& =\frac{{\pi }_k^{\mathrm{old}}\mathcal{S}\mathcal{N}\left({\mathbf{y}}_i|{\mathit{\eta }}_k^{\mathrm{old}},{\mathbf{\Sigma }}_k^{\mathrm{old}},{\mathit{\lambda }}_k^{\mathrm{old}}\right)}{\sum _{j=1}^K{\pi }_j^{\mathrm{old}}\mathcal{S}\mathcal{N}\left({\mathbf{y}}_i|{\mathit{\eta }}_j^{\mathrm{old}},{\mathbf{\Sigma }}_j^{\mathrm{old}},{\mathit{\lambda }}_j^{\mathrm{old}}\right)}\\ & ={\gamma }_{ik}^{\mathrm{old}},\end{array}$ where ${\gamma }_{ik}$ is the responsibility that the component k of the mixture takes for explaining the observation ${\mathbf{y}}_i$. Let ${\vartheta }_{ik}=E[Z_{ik}{U}_i|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}]$, and ${\psi }_{ik}=E[Z_{ik}{U}_i^2|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}]$. Then ${\vartheta }_{ik}$ can be simplified as follows: (7) $\begin{array}{rl}{\vartheta }_{ik}& =E\left[Z_{ik}E\left[U_i|\mathbf{Y},Z_{ik},{\mathbf{\Theta }}^{\mathrm{old}}\right]|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right]\\ & =E\left[U_i|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right]E\left[Z_{ik}|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right]\\ & ={\gamma }_{ik}^{\mathrm{old}}E\left[U_i|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right].\end{array}$(7) Using similar procedure, ${\psi }_{ik}$ may be given by (8) ${\psi }_{ik}={\gamma }_{ik}^{\mathrm{old}}E\left[U_i^2|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right].$(8)

If Equations (1(1) $f\left({\mathbf{y}}_i|\mathit{\eta },\mathbf{\Sigma },\mathit{\lambda }\right)=2\mathcal{N}\left({\mathbf{y}}_i|\mathit{\eta },\mathbf{\Sigma }\right)\mathrm{\Phi }\left({\mathit{\lambda }}^{\prime }{\mathbf{\Sigma }}^{-1/2}\left({\mathbf{y}}_i-\mathit{\eta }\right)\right),$(1) ) and (4(4) $\begin{array}{rl}{\mathbf{Y}}_i|U_i& =u_i\sim \mathcal{N}\left({\mathbf{y}}_i|\mathit{\eta }+u_i\mathit{\xi },\mathbf{\Omega }\right),\\ & U_i\sim \mathcal{H}\mathcal{N}\left(u_i|0,1,\left(0,\mathrm{\infty }\right)\right),\end{array}$(4) ) are satisfied, then using inverse matrix adjustment formula in [31] and matrix determinant lemma in [12], $U_i|{\mathbf{Y}}_i={\mathbf{y}}_i\sim \mathcal{T}\mathcal{N}\left(u_i|\frac{{\mathit{\xi }}^{\prime }{\mathbf{\Omega }}^{-1}\left({\mathbf{y}}_i-\mathit{\eta }\right)}{1+{\mathit{\xi }}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{\xi }},\frac{1}{1+{\mathit{\xi }}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{\xi }},\left(0,\mathrm{\infty }\right)\right),$ where $\mathcal{T}\mathcal{N}(u|\mu ,{\sigma }^2,(0,\mathrm{\infty }\left)\right)$ is a truncated normal with location parameter μ, scale parameter σ and support $(0,\mathrm{\infty })$. Based on the truncated normal probability distribution, (9) $\begin{array}{r}E\left[U_i|{\mathbf{Y}}_i={\mathbf{y}}_i,U_i>0\right]=\frac{1}{\beta }\left[\alpha +\frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}\right],\end{array}$(9) (10) $\begin{array}{rl}\mathrm{var}\left[U_i|{\mathbf{Y}}_i={\mathbf{y}}_i,U_i>0\right]& =\frac{1}{{\beta }^2}\left[1-\alpha \frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}-{\left(\frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}\right)}^2\right],\\ E\left[U_i^2|{\mathbf{Y}}_i={\mathbf{y}}_i,U_i>0\right]& =\frac{1}{{\beta }^2}\left[1+\alpha \frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}+{\alpha }^2\right],\end{array}$(10) where $\alpha =\frac{{\mathit{\xi }}^{\prime }{\mathbf{\Omega }}^{-1}\left({\mathbf{y}}_i-\mathit{\eta }\right)}{\sqrt{1+{\mathit{\xi }}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{\xi }}},\beta =\sqrt{1+{\mathit{\xi }}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{\xi }},$ and $\phi \left(\cdot \right)$ is a univariate standard normal density. The expected values $E\left[U_i|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right]\mathrm{and}E\left[U_i^2|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right]$ in Equations (7(7) $\begin{array}{rl}{\vartheta }_{ik}& =E\left[Z_{ik}E\left[U_i|\mathbf{Y},Z_{ik},{\mathbf{\Theta }}^{\mathrm{old}}\right]|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right]\\ & =E\left[U_i|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right]E\left[Z_{ik}|\mathbf{Y},{\mathbf{\Theta }}^{\mathrm{old}}\right]\\ & ={\gamma }_{ik}^{\mathrm{old}}E\left[U_i|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right].\end{array}$(7) ) and (8(8) ${\psi }_{ik}={\gamma }_{ik}^{\mathrm{old}}E\left[U_i^2|\mathbf{Y},Z_{ik}=1,{\mathbf{\Theta }}^{\mathrm{old}}\right].$(8) ) are obtained from Equations (9(9) $\begin{array}{r}E\left[U_i|{\mathbf{Y}}_i={\mathbf{y}}_i,U_i>0\right]=\frac{1}{\beta }\left[\alpha +\frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}\right],\end{array}$(9) ) and (10(10) $\begin{array}{rl}\mathrm{var}\left[U_i|{\mathbf{Y}}_i={\mathbf{y}}_i,U_i>0\right]& =\frac{1}{{\beta }^2}\left[1-\alpha \frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}-{\left(\frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}\right)}^2\right],\\ E\left[U_i^2|{\mathbf{Y}}_i={\mathbf{y}}_i,U_i>0\right]& =\frac{1}{{\beta }^2}\left[1+\alpha \frac{\phi \left(\alpha \right)}{\mathrm{\Phi }\left(\alpha \right)}+{\alpha }^2\right],\end{array}$(10) ) by replacing $\mathit{\eta }$, $\mathit{\xi }$, and $\mathbf{\Omega }$ with their corresponding estimates ${\mathit{\eta }}_k^{old}$, ${\mathit{\xi }}_k^{old}$, and ${\mathbf{\Omega }}_k^{old}$. The M-step of the EM-algorithm is given by ${\mathbf{\Theta }}^{\mathrm{new}}=\arg \underset{\mathbf{\Theta }}{max}\mathcal{Q}\left(\mathbf{\Theta },{\mathbf{\Theta }}^{\mathrm{old}}\right),$ and it is available in closed form. Using EM-algorithm, the estimates of the parameter $\mathbf{\Theta }$ can be updated as shown in Algorithm 1.

At mth iteration of the E-step, we need to compute the following expressions: $\begin{array}{rl}{\mathbf{\Sigma }}_k^{\left(m\right)}& ={\mathbf{\Omega }}_k^{\left(m\right)}+{\mathit{\xi }}_k^{\left(m\right)}{\mathit{\xi }}_k^{{\left(m\right)}^{\prime }};\\ {\mathit{\lambda }}_k^{\left(m\right)}& =\frac{{\left[{\mathbf{\Sigma }}_k^{\left(m\right)}\right]}^{-1/2}{\mathit{\xi }}_k^{\left(m\right)}}{\sqrt{1-{\mathit{\xi }}_k^{{\left(m\right)}^{\prime }}{\left[{\mathbf{\Sigma }}_k^{\left(m\right)}\right]}^{-1}{\mathit{\xi }}_k^{\left(m\right)}}};\\ {\alpha }_{ik}^{\left(m\right)}& =\frac{{\mathit{\xi }}_k^{{\left(m\right)}^{\prime }}{\left[{\mathbf{\Omega }}_k^{\left(m\right)}\right]}^{-1}\left({\mathbf{y}}_i-{\mathit{\eta }}_k^{\left(m\right)}\right)}{{\beta }_k^{\left(m\right)}};\\ {\beta }_k^{\left(m\right)}& =\sqrt{1+{\mathit{\xi }}_k^{{\left(m\right)}^{\prime }}{\left[{\mathbf{\Omega }}_k^{\left(m\right)}\right]}^{-1}{\mathit{\xi }}_k^{\left(m\right)}}.\end{array}$

In general, EM-algorithm converges to a local optimum. As a result, different initial values for the parameters are utilized during the estimation process to select the optimal estimates. A clustering method K-means has been employed to initialize the location parameters, the mixing coefficients, and the dispersion matrices. This clustering method is used to partition the observations into clusters in which each observation belongs to the cluster with the nearest mean. In this study, the clusters may represent tissue types or mixture of tissue types such as white matter, blood, water, fat, gray matter, air, cortical bone, cancellous bone, etc. We randomly initialize the remaining parameters of SGMM. Log-likelihood value cannot be computed analytically for MRF [22] and therefore, we cannot use it for selecting the optimal estimates. Following that we used mean squared error as criterion for selecting optimal estimates in SGMM in order to utilize the same criterion for the models used in this work. Two steps are employed during parameter estimation process. For a given number of tissue types, we estimated the parameters of the models and repeated this step to select the optimal parameter estimates using mean squared error. We selected the number of classes or tissue types through cross-validation using mean squared error. The stopping criterion for the convergence of the parameter estimation process is $\underset{ik}{max}\left|{\gamma }_{ik}^{\left(m+1\right)}-{\gamma }_{ik}^{\left(m\right)}\right|,$ with an upper limit $5\times {10}^{-5}$, where m denotes the iteration number of EM-algorithm.

2.5. Estimation of CT images

Let a d-dimensional vectors $\mathit{\eta }$, $\mathit{\nu }={\mathbf{\Sigma }}^{-1/2}\mathit{\lambda }$, and a $d\times d$ dispersion matrix $\mathbf{\Sigma }$ be partitioned as follows: ${\mathbf{Y}}_i=\begin{array}{}\end{array}\left[\begin{array}{c}{Y}_{i1}\\ {\mathbf{Y}}_{i2}\end{array}\right],\mathit{\eta }=\begin{array}{}\end{array}\left[\begin{array}{c}{\eta }_1\\ {\mathit{\eta }}_2\end{array}\right],\mathit{\nu }=\begin{array}{}\end{array}\left[\begin{array}{c}{\nu }_1\\ {\mathit{\nu }}_2\end{array}\right],\mathrm{and}\mathbf{\Sigma }=\begin{array}{}\end{array}\left[\begin{array}{cc}{\mathrm{\Sigma }}_{11}& {\mathbf{\Sigma }}_{12}\\ {\mathbf{\Sigma }}_{21}& {\mathbf{\Sigma }}_{22}\end{array}\right].$ The dimension of $Y_{i1}$, ${\eta }_1$, ${\nu }_1$, and ${\mathrm{\Sigma }}_{11}$ is a $1\times 1$. The random variable $Y_{i1}$ represents the $i^{th}$ voxel in CT image, and the random vector ${\mathbf{Y}}_{i2}$ denotes the corresponding voxel in MR images. If we assume that ${\mathbf{Y}}_i$ follows SGMM: (11) ${\mathbf{Y}}_i\sim \mathcal{S}\mathcal{N}\left({\mathbf{y}}_i|\mathit{\eta },\mathbf{\Sigma },\mathit{\lambda }\right),$(11) then ${\mathbf{Y}}_{i2}\sim \mathcal{S}\mathcal{N}\left({\mathbf{y}}_{i2}|{\mathit{\eta }}_2,{\mathbf{\Sigma }}_{22},\mathit{\tau }\right),$ where $\mathit{\tau }={\mathbf{\Sigma }}_{22}^{1/2}\left(\frac{{\mathbf{\Sigma }}_{22}^{-1}{\mathbf{\Sigma }}_{21}{\nu }_1+{\mathit{\nu }}_2}{\sqrt{1+{\nu }_1{\mathrm{\Sigma }}_{11}^c{\nu }_1}}\right)\mathrm{with}{\mathrm{\Sigma }}_{11}^c={\mathrm{\Sigma }}_{11}-{\mathbf{\Sigma }}_{12}{\mathbf{\Sigma }}_{22}^{-1}{\mathbf{\Sigma }}_{21}.$ According to Khounsiavash et al. [21], if Equation (11(11) ${\mathbf{Y}}_i\sim \mathcal{S}\mathcal{N}\left({\mathbf{y}}_i|\mathit{\eta },\mathbf{\Sigma },\mathit{\lambda }\right),$(11) ) holds, then the probability density function of $Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}$ can be given by (12) $f\left(y_{i1}|{\mathbf{y}}_{i2}\right)=\mathcal{N}\left(y_{i1}|{\eta }_1^c,{\mathrm{\Sigma }}_{11}^c\right)\frac{\mathrm{\Phi }\left({\mathit{\lambda }}^{\prime }{\mathbf{\Sigma }}^{-1/2}\left({\mathbf{y}}_i-\mathit{\eta }\right)\right)}{\mathrm{\Phi }\left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)},$(12) where ${\eta }_1^c={\eta }_1+{\mathbf{\Sigma }}_{12}{\mathbf{\Sigma }}_{22}^{-1}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\mathrm{and}{\mathrm{\Sigma }}_{11}^c={\mathrm{\Sigma }}_{11}-{\mathbf{\Sigma }}_{12}{\mathbf{\Sigma }}_{22}^{-1}{\mathbf{\Sigma }}_{21}.$

Proposition 2.1

Using Equation (12(12) $f\left(y_{i1}|{\mathbf{y}}_{i2}\right)=\mathcal{N}\left(y_{i1}|{\eta }_1^c,{\mathrm{\Sigma }}_{11}^c\right)\frac{\mathrm{\Phi }\left({\mathit{\lambda }}^{\prime }{\mathbf{\Sigma }}^{-1/2}\left({\mathbf{y}}_i-\mathit{\eta }\right)\right)}{\mathrm{\Phi }\left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)},$(12) ), the expected value of $Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}$ can be given by (13) $E\left[Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}\right]={\eta }_1^c+\frac{{\mathrm{\Sigma }}_{11}^c{\nu }_1}{\sqrt{1+{\nu }_1{\mathrm{\Sigma }}_{11}^c{\nu }_1}}\frac{\phi \left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}{\mathrm{\Phi }\left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}.$(13)

Proof.

Let $\kappa \left({\mathbf{y}}_{i2}\right)=\mathrm{\Phi }\left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\right({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\left)\right)$ and $B={\mathit{\nu }}_2^{\prime }({\mathbf{y}}_{i2}-{\mathit{\eta }}_2)-{\nu }_1{\eta }_1$. The expected value of $Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}$ can be given by $\begin{array}{rl}E\left[Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}\right]& ={\int }_{\mathbb{R}}{y}_{i1}f\left(y_{i1}\mid {\mathbf{y}}_{i2}\right)\mathrm{d}{y}_{i1}\\ & =\frac{1}{\kappa \left({\mathbf{y}}_{i2}\right)}{\int }_{\mathbb{R}}{y}_{i1}\mathcal{N}\left(y_{i1}\mid {\eta }_1^c,{\mathrm{\Sigma }}_{11}^c\right)\mathrm{\Phi }\left({\mathit{\lambda }}^{\prime }{\mathbf{\Sigma }}^{-1/2}\left({\mathbf{y}}_i-\mathit{\eta }\right)\right)\mathrm{d}{y}_{i1}\\ & =\frac{1}{\kappa \left({\mathbf{y}}_{i2}\right)}{\int }_{\mathbb{R}}{\int }_0^{\mathrm{\infty }}{y}_{i1}\mathcal{N}\left(y_{i1}\mid {\eta }_1^c,{\mathrm{\Sigma }}_{11}^c\right)\mathcal{N}\left(x\mid {\nu }_1{y}_{i1}+B,1\right)\mathrm{d}x\mathrm{d}{y}_{i1}\\ & =\frac{1}{\kappa \left({\mathbf{y}}_{i2}\right)}{\int }_0^{\mathrm{\infty }}\mathcal{N}\left(x\mid B+{\nu }_1{\eta }_1^c,1+{\nu }_1{\mathrm{\Sigma }}_{11}^c{\nu }_1\right)E[Y_{i1}|X=x]\mathrm{d}x,\end{array}$ where $Y_{i1}|X=x\sim \mathcal{N}\left(x\mid {\eta }_1^c+\mathrm{\Lambda }{\nu }_1\left(x-B-{\nu }_1{\eta }_1^c\right),\mathrm{\Lambda }\right)\mathrm{with}\mathrm{\Lambda }={\left(\frac{1}{{\mathrm{\Sigma }}_{11}^c}+{\nu }_1^2\right)}^{-1}.$ The result in Equation (13(13) $E\left[Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}\right]={\eta }_1^c+\frac{{\mathrm{\Sigma }}_{11}^c{\nu }_1}{\sqrt{1+{\nu }_1{\mathrm{\Sigma }}_{11}^c{\nu }_1}}\frac{\phi \left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}{\mathrm{\Phi }\left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}.$(13) ) follows from well known properties of the truncated normal distribution and inverse matrix adjustment formula.

Using Equation (13(13) $E\left[Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}\right]={\eta }_1^c+\frac{{\mathrm{\Sigma }}_{11}^c{\nu }_1}{\sqrt{1+{\nu }_1{\mathrm{\Sigma }}_{11}^c{\nu }_1}}\frac{\phi \left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}{\mathrm{\Phi }\left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}.$(13) ), we can obtain the point estimator of $Y_{i1}$ by $\begin{array}{rl}E\left[Y_{i1}|{\mathbf{Y}}_{i2}\right]& =E\left[E\left[Y_{i1}|{\mathbf{Y}}_{i2},Z_i\right]|{\mathbf{Y}}_{i2}\right]\\ & =\sum _{k=1}^{K}P\left(Z_i=k|{\mathbf{Y}}_{i2}\right)E\left[Y_{i1}|{\mathbf{Y}}_{i2},Z_i=k\right].\end{array}$ In this framework, the latent variable $Z_i$ represents the underlying tissue classes. The weight $P(Z_i=k|{\mathbf{Y}}_{i2})$ can be computed using Bayes' theorem. The expected value $E[Y_{i1}|{\mathbf{Y}}_{i2},Z_i=k]$ is obtained from Equation (13(13) $E\left[Y_{i1}|{\mathbf{Y}}_{i2}={\mathbf{y}}_{i2}\right]={\eta }_1^c+\frac{{\mathrm{\Sigma }}_{11}^c{\nu }_1}{\sqrt{1+{\nu }_1{\mathrm{\Sigma }}_{11}^c{\nu }_1}}\frac{\phi \left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}{\mathrm{\Phi }\left({\mathit{\tau }}^{\prime }{\mathbf{\Sigma }}_{22}^{-1/2}\left({\mathbf{y}}_{i2}-{\mathit{\eta }}_2\right)\right)}.$(13) ) by indexing the parameters with k.

2.6. Model validation

The main focus of this work is to investigate CT image predictive quality of SGMM and GMM by partitioning the data into two major tissue types. It is also aimed at comparing the estimation quality of the partitioning approach with GMM, HMM, and MRF on the bone tissues.

We used leave-one-out cross-validation method to compare the predictive quality of the models. That is, we keep one head to validate the models and use the remaining heads for training the models. Using the threshold CT image intensity, CT image in a validation dataset is partitioned into non-bone and bone tissues. Let $Y_i^{\left(t\right)}$ and ${\hat{Y}}_i^{\left(t\right)}$ be CT image and its corresponding estimated CT image intensities in a given tissue region t. One can use the square and the absolute loss functions to assess the estimation cost. They can be given by ${\left({\hat{Y}}_i^{\left(t\right)}-Y_i^{\left(t\right)}\right)}^2\mathrm{and}|{\hat{Y}}_i^{\left(t\right)}-Y_i^{\left(t\right)}|,$ where t=1,2. Since the mean square error heavily weights the outliers, the mean absolute error (MAE) can be employed as a main tool to evaluate the estimation performance of the models. The MAE can be given by $MA{E}_t=\frac{1}{n_t}\sum _{i=1}^{n_t}|{\hat{Y}}_i^{\left(t\right)}-Y_i^{\left(t\right)}|,$ where $n_t$ is the number of voxels in partition t.

The peak signal-to-noise ratio (PSNR) can also be used to quantify the overall quality of the estimation. It takes square loss function into account through the mean squared error. The PSNR may be given by $PSNR=10{\log }_{10}\left(\frac{n{M}^2}{\sum _{i=1}^n{\left({\hat{Y}}_i-Y_i\right)}^2}\right),$ where $Y_i$ and ${\hat{Y}}_i$ represent CT image and the estimated CT image intensities at voxel i, respectively, and M is the maximal intensity in CT image. We exploited patient-wise leave-one-out cross-validation and mean absolute error to compare the estimation quality of the models. One can also use peak signal-to-noise ratio to compare the estimation performance of the models. The better model has lower MAE and higher PSNR. Since MAE and PSNR are crude estimation quality measures, we utilized smoothed residual and absolute residual plots to further evaluate the estimation quality of the models through the tissues of the heads. A moving average over non-overlapping windows in CT image intensities can be used as a main tool to investigate and identify the model that outperforms in bone tissue estimation. Over the non-overlapping windows on $Y_i$ with a window size of 20 HU, the averages for $Y_i$, ${\hat{Y}}_i-Y_i$, and $|{\hat{Y}}_i-Y_i|$ over the windows can be computed to obtain the smoothed residual and absolute residual plots.

In addition to the above model performance evaluation methods, we assessed the performance of the models using Bland-Altman plot [6,7], which is a graphical method to compare two measurements by plotting the differences between the two measurements ${\hat{Y}}_i-Y_i$ against their averages $({\hat{Y}}_i+Y_i)/2$. If one of the two measurements is a reference measurement, then the differences can be plotted against the reference measurement [23] and that coincides with the smoothed residual plot via cross-validation. The main advantage of the Bland-Altman plot is that it reveals a relationship between the differences and the magnitude of the measurements, to examine possible systematic bias, and outliers.

We can also complement the MAE based performance evaluation of the methods by Wilcoxon signed-rank test, which is a non-parametric alternative to the paired Student's t-test. To compare SGMM to every one of the remaining methods, the paired method-method differences in MAEs were examined by using the one-sided Wilcoxon signed rank test, with the null hypothesis that SGMM is the same as one of the other methods with respect to the MAE values and the research hypothesis that SGMM is better than the method being compared.

3. Results

A CT image intensity 100 HU was utilized as a delimiting value during CT image estimation. The optimal parameter estimates were received for K=8 in GMM and K=5 for both HMM and MRF. In the case of SGMM and GMM${}^{\ast }$, we have received the optimal parameter estimates for K=6 for both major tissue types. Table 1 demonstrates a summary of mean absolute errors for the bone tissues and p-values of the Wilcoxon signed-rank test.

Table 1. Mean absolute errors and p-values for bone tissues.

	Model
Head	SGMM	GMM ${}^{\ast }$	HMM	MRF	GMM
1	316.25	315.16	324.94	307.95	314.41
2	349.52	348.05	360.03	328.89	365.12
3	303.58	301.78	331.16	322.48	328.01
4	272.21	269.33	296.04	280.63	292.78
5	350.41	349.28	366.13	357.22	359.56
Average	318.39	316.72	335.66	319.43	331.98
p-Value	–	0.985	0.030	0.500	0.053

The rows of the table represent validation datasets. The table presents mean absolute errors received for the models. In terms of MAE, it is apparent that for each head the partitioning approach (SGMM and GMM${}^{\ast }$) and MRF had better performance than HMM and GMM. Moreover, the differences of the average MAEs show that the partitioning approach and MRF achieved better results than HMM and GMM. For instance, SGMM outperforms HMM with the method-to-method difference of average MAEs $-17.27$ and the standard deviation of the method-to-method differences of MAEs 8.23. The p-values of Wilcoxon signed-rank test show that SGMM outperforms HMM and GMM significantly on bone tissues.

The results in Table 2 show that the partitioning approach has noticeable outperformance on dense bone tissues (approximately with CT image intensities greater than 900 HU according to Washington and Leaver [34]) as compared to the remaining models. Even though MRF is computationally expensive, it had better performance than HMM and GMM on dense bone tissues. The p-values of Wilcoxon signed-rank test reveal that SGMM has significantly better performance than the remaining methods except for GMM${}^{\ast }$ on dense bone tissues.

Table 2. Mean absolute errors and p-values for dense bone tissues.

	Model
Head	SGMM	GMM ${}^{\ast }$	HMM	MRF	GMM
1	365.68	361.27	418.99	401.43	407.85
2	406.09	403.65	458.23	394.89	492.96
3	331.85	328.44	407.01	404.80	386.72
4	236.89	232.15	290.42	295.46	290.77
5	436.42	433.84	509.20	488.10	505.49
Average	355.39	351.87	416.77	396.94	416.76
p-Value	–	0.985	0.030	0.053	0.030

Table 3 demonstrates the estimation quality of the models on non-bone tissues. The best result was received for HMM. However, there is already a good contrast between soft tissues and air on MR images. The remaining models had similar behavior on non-bone tissues.

Table 3. Mean absolute errors and p-values for non-bone tissues.

	Model
Head	SGMM	GMM ${}^{\ast }$	HMM	MRF	GMM
1	111.96	112.67	97.75	106.15	114.10
2	116.56	117.53	99.69	103.60	116.12
3	124.65	125.46	94.93	117.31	122.00
4	116.61	117.19	101.54	112.52	111.76
5	122.08	122.42	98.87	126.51	118.19
Average	118.37	119.05	98.55	113.22	116.43
p-Value	–	0.030	0.985	0.947	0.947

Table 4 presents the overall summary of CT image estimation quality. In comparison to the other models, we received better result for HMM. The reason is that HMM outperformed the other models on non-bone tissues. However, this is not the main interest in this work.

Table 4. Combined mean absolute errors and p-values.

	Model
Head	SGMM	GMM ${}^{\ast }$	HMM	MRF	GMM
1	144.15	144.58	133.56	137.95	145.67
2	151.44	152.04	138.68	137.34	153.40
3	161.00	161.28	142.92	158.99	163.85
4	146.71	146.62	139.17	145.05	146.78
5	160.53	160.62	143.89	165.36	158.85
Average	152.77	153.03	139.65	148.94	153.71
p-Value	–	0.069	0.985	0.911	0.140

Table 5 demonstrates the prediction accuracy of the models in terms of PSNR. The results show that the models had similar behavior. On the bone tissues, PSNR has shown that the models have similar estimation performance.

Table 5. Evaluation of the methods based on PSNR and p-values.

	Model
Head	SGMM	GMM ${}^{\ast }$	HMM	MRF	GMM
1	19.92	19.90	20.29	20.30	20.29
2	19.56	19.50	20.23	20.40	20.15
3	19.75	19.71	20.13	20.11	19.69
4	19.76	19.76	20.43	20.46	20.69
5	19.59	19.58	19.78	19.63	19.66
Average	19.72	19.69	20.17	20.18	20.09
p-Value	–	0.978	0.030	0.030	0.053

On the basis of average, mean absolute error was utilized to compare CT image prediction accuracy of the models. Smoothed absolute residual plots were also employed to assess the estimation quality of the models. In comparison to MAE, smoothed absolute residual plots are powerful to evaluate the estimation performance of the models through the tissues of the heads. Figure 2 presents smoothed absolute residual plots for the models. The absolute residuals were averaged over non-overlapping windows in CT image intensities with window size 20 HU.

Figure 2. Smoothed absolute residual plot for the five patients.

It is clear from Figure 2 that none of the models outperformed throughout the tissues of the head. However, SGMM and GMM${}^{\ast }$ outperformed the other models on dense bone tissues. Figure 3 shows smoothed residual plot. The deviations of observed CT image intensities from the estimated CT image intensities were exploited to obtain average residuals over non-overlapping windows in CT image intensities with window size 20 HU. It is evident from the plot that bone tissues were underestimated. However, the partitioning approach has improved underestimation of bone tissues.

Figure 3. Smoothed residual plot for the five patients.

The Bland-Altman plots of the methods are shown in Figure 4. It shows the bias of the methods in estimating the CT images. The bias is higher for higher CT image values. However, the partitioning approach has improved the bias in bone tissue estimation as compared to the remaining methods. The Bland-Altman plots also show higher variability on lower CT image values. In comparison to Figure 3, we observe that the higher variability on lower CT image values in Figure 4 is due the estimated CT image values. Notice, however, that this work is mainly concerned with improving bone tissue estimation.

Figure 4. Bland-Altman plot for the five patients.

The partitioning approach had better performance in tracing bone tissues. Bone tissues appear as bright white on CT images, see Figure 5. The figure shows slices of CT image and its corresponding predicted slices for a representative patient. The top right portion of the images in the first row of the figure clearly shows that the partitioning approaches are better in identifying bone tissues.

Figure 5. The first column (a) presents the original CT image slices and the remaining columns (b)–(f) show the corresponding predicted slices of CT image.

We presented the prediction errors in Figure 5, which corresponds to the predicted slices in Figure 6. It can be seen that the images of the prediction errors corresponding to the bone tissues appear darker for the partitioning approach: GMM${}^{\ast }$ and SGMM.

Figure 6. The first column (a) presents the original CT image slices and the remaining columns (b)–(f) show the prediction errors for each model.

4. Discussion

Statistical model for voxel-wise CT image estimation from MR images have been presented and evaluated using cross-validation on five datasets. Kuljus et al. [22] and Johansson et al. [16] have used voxel-wise estimation approach to study CT image estimations. Existing works suggest that the estimation quality on the bone tissues is poor. This study was aimed at probing the estimation of CT images by partitioning the data into two major tissue types: non-bone and bone tissues. Specifically, the focus of the study was to obtain improved bone tissue estimations. We proposed SGMM to relax the distributional assumption utilized by Kuljus et al. [22] and Johansson et al. [16]. The model was motivated to take the asymmetrical distribution that could arise from the partitioning of the data into account. We also used GMM in addition to SGMM to examine the effect of the distributional nature of the partitioned data on the estimation process. The study was also aimed at comparing CT image estimation performance of SGMM, HMM, GMM, and MRF on the bone tissues. In MRF and HMM, spatial dependence between neighboring voxels has been taken into account.

EM-algorithm was developed for SGMM parameter estimation. EM-gradient algorithm was used to estimate the parameters of MRF while EM-algorithm was utilized to estimate the parameters of GMM and HMM. For MRF, the updates of the parameters at M-step of EM-algorithm were difficult to obtain in an explicit form and thereby, a gradient based optimization was utilized during parameter estimation process. Moreover, Gibbs sampling was used at E-step of the algorithm. Thus, the estimation process was expensive in MRF. Unlike the other models, log-likelihood function in MRF involves Gibbs field and it is not computable. This means that log-likelihood based selection of optimal parameter estimate is not feasible analytically. Consequently, mean squared error was employed in selecting the optimal parameter estimates of the models.

Table 4 demonstrated that HMM outperformed the other models. This was essentially due to its best estimation performance on non-bone tissues, see Table 3. According to Karlsson et al. [20], the key task in CT image estimation from MR images is to obtain an improved bone tissue estimation. Table 1 revealed that GMM${}^{\ast }$ (GMM trained on each major tissue type), SGMM, and MRF had better prediction accuracy than HMM and GMM on the bone tissues. In addition, Table 2 shows that GMM${}^{\ast }$ and SGMM had noticeable outperformance on dense bone tissues. Based on the mean absolute error, HMM and GMM perform similarly on the bone tissues.

To provide statistical evidence for the CT image estimation performance of SGMM in comparison with the remaining methods, the Wilcoxon signed-rank tests have been carried out. The tests have shown that SGMM has better performance on bone tissues than HMM and GMM. In addition, the Wilcoxon signed-rank tests provided an evidence that SGMM has significant outperformance than MRF, HMM, and GMM on dense bone tissues.

The skewness assumption allowed to recognize skewness in the partitions of the data. The estimates of the skewness parameters demonstrated that the partitions have skewness property and the nature of skewness was dependent on the subtissue types. The estimates of the skewness parameters ranged from $-1.54$ to 3.03. In this particular application, however, the skewness assumption did not improve CT image estimation quality as compared to the symmetric assumption in GMM${}^{\ast }$, see Tables 1–5.

Figure 2 shows that the models performed better on soft tissues, that is tissues having CT image intensity closer to 0 HU. On the other hand, the models had weaker prediction accuracy on the two extremes that is on air and bone tissues. This pattern of residual plots has been observed in [16,22]. Moreover, the figure shows that none of the models outperformed throughout the tissues of the head. Tables 1–4 demonstrate that the predictive accuracy of the models was dependent on the heads. That is the results were not uniform over the heads. This might have a problem in real applications and needs a further investigation.

The bias of the methods in estimation of the CT images was manifested in the smoothed residual and Bland-Altman plots in Figures 3 and 4. Though the bias is higher on bone tissues, it was improved by our approach as compared to the remaining methods. This was the main concern of this work. The Bland-Altman plots also demonstrated higher variability of the estimation on lower CT image values.

5. Conclusions

In this study, we examined CT image estimation by partitioning the data into two major tissue types. This partitioning approach is an efficient way to get a good quality CT image substitute with improved estimation of bone tissues. Moreover, the SGMM and the developed algorithm to estimate its parameters is general and can be applied to other applications.

Acknowledgments

The authors thank Adam Johansson and Thomas Asklund for providing us data and David Bolin for providing us Mathlab code for MRF. Moreover, the authors thank Kristi Kuljus for HMM results. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at High Performance Computing Center North (HPC2N).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work is supported by the Swedish Research Council grant (reg. no. 340-2013-5342).