R/fuzzySpatialCMeansSegmentation.R

In stnava/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' fuzzySpatialCMeansSegmentation
#'
#' Fuzzy spatial c-means for image segmentation.
#'
#' Image segmentation using fuzzy spatial c-means as described in
#'
#' Chuang et al., Fuzzy c-means clustering with spatial information for image
#' segmentation.  CMIG: 30:9-15, 2006.
#'
#' @param image image to be segmented.
#' @param mask optional mask image.  Otherwise, the entire image is used.
#' @param numberOfClusters number of segmentation clusters
#' @param m fuzziness parameter (default = 2).
#' @param p membership importance parameter (default = 1).
#' @param q spatial constraint importance parameter (default = 1).
#' \code{q = 0} is equivalent to conventional fuzzy c-means.
#' @param radius neighborhood radius (scalar or array) for spatial
#' constraint.
#' @param maxNumberOfIterations iteration limit (default = 20).
#' @param convergenceThreshold Convergence between iterations is measured
#' using the Dice coefficient (default = 0.02).
#' @param verbose print progress.
#' @return list containing segmentation and probability images
#'
#' @author NJ Tustison
#'
#' @examples
#' image <- antsImageRead(getANTsRData("r16"))
#' mask <- getMask(image)
#' fuzzy <- fuzzySpatialCMeansSegmentation(image, mask)
#'
#' @export fuzzySpatialCMeansSegmentation

fuzzySpatialCMeansSegmentation <- function(
    image, mask = NULL, numberOfClusters = 4,
    m = 2, p = 1, q = 1, radius = 2, maxNumberOfIterations = 20, convergenceThreshold = 0.02,
    verbose = FALSE) {
  if (is.null(mask)) {
    mask <- antsImageClone(image) * 0 + 1
  }

  x <- image[mask != 0]

  v <- seq(from = 0, to = 1, length.out = numberOfClusters + 2)[2:(numberOfClusters + 1)]
  v <- v * (max(x) - min(x)) + min(x)
  cc <- length(v)

  if (verbose == TRUE) {
    cat("Initial cluster centers: ", v, "\n")
  }

  xx <- matrix(data = 0, nrow = cc, ncol = length(x))
  for (i in seq.int(cc))
  {
    xx[i, ] <- x
  }

  if (length(radius) == 1) {
    radius <- rep(radius, image@dimension)
  }

  segmentation <- antsImageClone(image) * 0
  probabilityImages <- NULL

  iter <- 0
  diceValue <- 0
  while (iter < maxNumberOfIterations && diceValue < 1.0 - convergenceThreshold) {
    # Update membership values

    xv <- matrix(data = 0, nrow = cc, ncol = length(x))
    for (k in seq.int(cc))
    {
      xv[k, ] <- abs(x - v[k])
    }

    u <- matrix(data = 0, nrow = nrow(xv), ncol = ncol(xv))
    for (i in seq.int(cc))
    {
      n <- xv[i, ]

      d <- n * 0
      for (k in seq.int(cc))
      {
        d <- d + (n / xv[k, ])^(2 / (m - 1))
      }
      u[i, ] <- 1 / d
    }
    u[is.nan(u)] <- 1


    # Update cluster centers

    v <- rowSums((u^m) * xx, na.rm = TRUE) / rowSums(u^m, na.rm = TRUE)

    if (verbose == TRUE) {
      cat("Updated cluster centers: ", v, "\n")
    }

    # Spatial function

    h <- matrix(data = 0, nrow = nrow(u), ncol = ncol(u))
    for (i in seq.int(cc))
    {
      uImage <- antsImageClone(image) * 0
      uImage[mask != 0] <- u[i, ]
      uNeighborhoods <- getNeighborhoodInMask(uImage, mask, radius)
      h[i, ] <- colSums(uNeighborhoods, na.rm = TRUE)
    }

    # u prime

    d <- rep(0, ncol(u))
    for (k in seq.int(cc))
    {
      d <- d + (u[k, ]^p) * (h[k, ]^q)
    }

    probabilityImages <- list()
    uprime <- matrix(data = 0, nrow = nrow(u), ncol = ncol(u))
    for (i in seq.int(cc))
    {
      uprime[i, ] <- (u[i, ]^p) * (h[i, ]^q) / d
      uprimeImage <- antsImageClone(image) * 0
      uprimeImage[mask != 0] <- uprime[i, ]
      probabilityImages[[i]] <- uprimeImage
    }

    tmpSegmentation <- antsImageClone(image) * 0
    tmpSegmentation[mask != 0] <- max.col(t(uprime))

    diceValue <- labelOverlapMeasures(segmentation, tmpSegmentation)$MeanOverlap[1]
    iter <- iter + 1

    if (verbose == TRUE) {
      cat("Iteration ", iter, " (out of ", maxNumberOfIterations, "):  ",
        "Dice overlap = ", diceValue, "\n",
        sep = ""
      )
    }
    segmentation <- tmpSegmentation
  }
  return(list(
    segmentationImage = segmentation,
    probabilityImages = probabilityImages
  ))
}