initNormMix {mritc}

Get the Initial Estimate of the Parameters of a Normal Mixture Model

Description

Obtain initial estimation of proportions, means, and standard deviations of different components (tissue types for MRI) based on threshold values generated by Otsu's method implemented by a fast algorithm, or proportion of different components.

Usage

   initOtsu(y, m)
   initProp(y, prop)

Arguments

y

a vector of intensity values of an image.

m

number of classes (tissue types for MRI) minus 1.

prop

the initial estimate of proportion of different components.

Details

The exhaustive search part of the function for Otsu's algorithm is adapted from combn. For initProp, the threshold values are quantiles based on the initial estimate of proportion of different components.

Value

prop

a vector of initial estimate of the proportions of different components (tissue types for MRI).

mu

a vector of initial estimate of the means of different components (tissue types for MRI).

sigma

a vector of initial estimates of the standard deviations of different components (tissue types for MRI).

Note

For initOtsu, it supports any number of m. However, for MRI data, it can be slow if m is bigger than 3 even with the fast algorithm implemented, since the Otsu's algorithm uses an exhaustive search. But it should be fine with m equal to 2, which corresponds to the typical case in MRI classification with three major tissue types CSF, GM, and WM.

References

Nobuyuki Otsu (1979). A threshold selection method from gray-level histograms IEEE Transactions on Systems, Man and Cybernetics vol. 9, 62-66

Ping-Sung Liao, Tse-Sheng Chen and Pau-Choo Chung (2001) A Fast Algorithm for Multilevel Thresholding Journal of Information Science and Engineering vol. 17, 713-727

Examples

  #Example 1
  prop <- c(.3, .4, .3)
  mu <- c(40,  90, 130)
  sigma <- c(4, 8, 4)
  y <- floor(rnormmix(n=100000, prop, mu, sigma)[,1])
  initOtsu(y, 2)
  initProp(y, prop)

  #Example 2
  T1 <- readMRI(system.file("extdata/t1.rawb.gz", package="mritc"),
                c(91,109,91), format="rawb.gz")
  mask <-readMRI(system.file("extdata/mask.rawb.gz", package="mritc"),
                 c(91,109,91), format="rawb.gz")
  initOtsu(T1[mask==1], 2)
  initProp(T1[mask==1], c(0.17, 0.48, 0.35))