Automated and optimal detection of 3D articular cartilage using undecimated wavelets in MRI

Abstract

Articular cartilage is a thin and curvilinear structure whose image may get corrupted with noise during acquisition process. Image segmentation of this structure depends on expertise of an operator. In this paper, we propose an automated cartilage detection technique using 3D histogram and wavelet multiscale singularity analysis. 3D undecimated wavelet transform is implemented on MRI volume to obtain wavelet coefficients, which are used to determine an adaptive threshold for a given local resolution and a global threshold value. Local threshold helps to segment foreground cartilage edge details and is obtained using maximum likelihood estimate of the 3D histogram for a selected confidence level of intensity histogram. A global threshold is used to optimize coefficients using wavelet multiresolution singularity. The final wavelet coefficients are used to obtain a 3D model of cartilage tissue. The proposed method has been tested and validated using MRI and phantom datasets of articular cartilage. Quantitative analysis has been performed using mean square error (MSE), signal-to-noise ratio (SNR) and volumetric estimation of the datasets for different confidence and noise levels. The proposed method displays reduction in MSE for both denoised and noisy MRI volumes at different standard deviations of noise with overall improvement in SNR.

Appendices

Appendix 1: Rayleigh intensity interval estimation

To comput confidence interval for unknown mean $\mu$ for confidence level $\alpha 1$ for Rayleigh distribution, we assume the distribution to be sampled with respect to normal distribution.

$$\begin{aligned} -t_{1 - \alpha 1/2} \le \frac{ \mu _\mathrm{r} - \mu }{S/\sqrt{N}} \le t_{1 - \alpha 1/2} \end{aligned}$$

(17)

where sample mean $\mu _\mathrm{r}$ is obtained with unbiased estimator $\hat{\sigma }$, which is obtained using MLE of Rayleigh probability distribution given by Eq. (5) [26, 27]. For N independent samples, the likelihood equation for the distribution is given as [26],

$$\begin{aligned} L(P) = \prod ^{N}_{i}\left[ \frac{x}{ \sigma ^{2}} e^{- \frac{(x )^{2}}{2\sigma ^{2}}} \right] = \left( \frac{x}{\sigma ^{2}}\right) ^{N}\prod ^{N}_{i} e^{- \frac{(x )^{2}}{2\sigma ^{2}}} \end{aligned}$$

(18)

To obtain MLE for the distribution for unknown variance, take partial derivative with respect to $\sigma ^{2}$ of the log-likelihood equation. The MLE for Rayleigh is given as,

$$\begin{aligned} \hat{\sigma ^{2}}= q = \frac{1}{2N}\sum _{i = 1}^{N}x_{i}^{2} \end{aligned}$$

(19)

Mean computed using this estimate is then given as,

$$\begin{aligned} \mu _\mathrm{r} =\hat{\sigma }\sqrt{\frac{\pi }{2}} \end{aligned}$$

(20)

To compute confidence limits using mean, rearrange the terms in the above Eq. (18)

$$\begin{aligned} \mu _\mathrm{r}\pm t_{1 - \alpha 1/2}\left( \frac{S}{\sqrt{N}}\right) \end{aligned}$$

(21)

where $t_{1 - \alpha 1/2}$ can be obtained from the t-distribution table under given confidence level $\alpha 1$ for $N-1$ degrees of freedom.

Appendix 2: Gaussian intensity interval estimation

To compute confidence interval for Gaussian distribution with N independent observations for unknown mean, we use sampling for the given distribution from the normal distribution. The sampling distribution follows a t-distribution with $N-1$ degree of freedom and is given as follows [27],

$$\begin{aligned} T =\frac{ \hat{X} - \mu }{S/\sqrt{N}} \end{aligned}$$

(22)

we determine the $t_{(1 - \alpha 1/2 ),(N - 1)}$ value from the t-table for $N - 1$ degree of freedom and is given as follows [27],

$$\begin{aligned} -t_{1 - \alpha 1/2} \le \frac{ \hat{X} - \mu }{S/\sqrt{N}} \le t_{1 - \alpha 1/2} \end{aligned}$$

(23)

where $\hat{X}$ is also the MLE for Gaussian distribution. Probability distribution for Gaussian distribution is given in Sect. 2 [27],

The likelihood equation for Gauss distribution with N independent samples is [36],

$$\begin{aligned} L(P) = \prod _{i=1}^{N} \left[ \frac{1}{\sqrt{2\pi \sigma ^{2}}} e\left( - \frac{x_{i} - \mu }{2\sigma ^{2}}\right) \right] \end{aligned}$$

(24)

The MLE of the above equation with respect to mean is given as follows [36],

$$\begin{aligned} \hat{X} = \mu = \frac{1}{N}\sum _{i = 1}^{N}x_{i} \end{aligned}$$

(25)

Appendix 3: Wavelet interscale difference

Wavelet coefficients are obtained by convolution of wavelet function on 2D image function I(x, y) as [23]

$$\begin{aligned} w_{j} = \frac{1}{\sqrt{MN}} \sum _{x=0}^{M-1}\sum _{y=0}^{N-1}I(x,y)*\psi _{j,m,n}^{i}(x,y), i = \{H,V,D\}\nonumber \\ \end{aligned}$$

(26)

But image function I(x, y) follows a Rice probability distribution function given in Eq. (4) [25].

Thus wavelet coefficients obtained from the image are assumed to have a Rice distribution function, given as

$$\begin{aligned} P(w_{j}|I,\sigma ) = \frac{w_{j}}{\sigma ^{2}}\exp \left( -\frac{w_{j}^{2} + A_\mathrm{sig}^{2}}{2\sigma ^{2}}\right) I_{0}\frac{w_{j}A_\mathrm{sig}}{\sigma ^{2}} \end{aligned}$$

(27)

We know in wavelet domain the modulus of sum of wavelet coefficients can be given as [32],

$$\begin{aligned} N_{j}I(x_{0}) =\acute{A}s^{\alpha + 1} \approx P(w_{j}|I,\sigma ) \end{aligned}$$

(28)

And the interscale difference between wavelet coefficients at two scale is given as,

$$\begin{aligned} N_{j + 1}I(x_{0}) - N_{j}I(x_{0}) =\gamma \end{aligned}$$

(29)

But since the wavelet coefficients are governed by Rice probability distribution, the interscale difference is now given as,

$$\begin{aligned} \grave{\gamma } =\gamma \sigma ^{2}= & {} \left[ w_{j+1}exp\left( -\frac{w_{j + 1}^{2} + A_\mathrm{sig}^{2}}{2\sigma ^{2}}\right) I_{0}\frac{w_{j + 1}A_\mathrm{sig}}{\sigma ^{2}} \right. \nonumber \\&\left. -\,w_{j}exp(-\frac{w_{j}^{2} + A_\mathrm{sig}^{2}}{2\sigma ^{2}}) I_{0}\frac{w_{j}A_\mathrm{sig}}{\sigma ^{2}} \right] \end{aligned}$$

(30)