dmri3d: Read MATLAB-Format Data from Mri3d.exe

#' @title Read MATLAB-Format data from mri3d.exe
#' @description The MATLAB program \code{mri3d.exe} is used for MRI segmentation at
#' the University Hospital of Zuerich. With \code{File/Save/Measurements}, a MATLAB
#' file (default extension \code{.mat}) with segmentation slices and 3D-mesh data
#' can be exported  and imported into R by function \code{readMri3d}.
#' The sampling of slices can be optionally decimated.
#' @param filename Name of file exported with from mri3d.exe
#' @param dmin Smallest permissible length of an edge. Used in function
#' \code{simplify.owin} in package \code{spatstat}. No reduction when \code{dmin==NULL}.
#' @return List with the following items:
#' \code{vars}: parameter settings used for mri3d.exe, and record data: examination date,
#' protocol name, examination name, filename, directory; and from 3D analysis:
#' vol3d (Volume) and area3d (surface area).
#'
#' \code{slices}: data frame with slice number, segment number, x and y
#'
#' \code{slicearea}: data frame with slice, segment and area
#'
#' \code{indices}: 3D integer array of indices of hull triangles; can be NULL
#'
#' \code{vertices}: 3D numeric array of vertex positions; can be NULL
#'
#' @examples
#' filename = system.file("extdata", "rec1.mat", package = "Dmri3d")
#' mri = readMri3d(filename)
#' str(mri)
#' \dontrun{
#' library(rgl)
#' open3d()
#' wire3d(with(mri,tmesh3d(vertices[[1]], indices[[1]], homogeneous=FALSE)))
#' wire3d(with(mri,tmesh3d(vertices[[2]], indices[[2]], homogeneous=FALSE)))
#' wire3d(with(mri,tmesh3d(vertices[[3]], indices[[3]], homogeneous=FALSE)))
#' library(Rvcg)
# Rvcg functions need 4D homogeneous coordinates
#' plotSimplified = function(mri,index){
#'   open3d()
#'   mesh3d = list(vb=rbind(mri$vertices[[index]],
#'                          rep(1,ncol(mri$vertices[[index]]))),
#'                 it=mri$indices[[index]])
#'   class(mesh3d) = "mesh3d"
#'   mesh3dA = vcgQEdecim(mesh3d, percent = 0.1, normcheck = TRUE)
#'   mesh3dB = vcgSmooth(mesh3d, iteration=1)
#'   wire3d(mesh3dA)
#' }
#' plotSimplified(mri,1)
#' plotSimplified(mri,2)
#' plotSimplified(mri,3)
#' }
#' @import spatstat
#' @import R.matlab
#' @import assertthat
#' @export
readMri3d = function(filename, dmin = NULL) {
#  if (FALSE) {
#    filename = system.file("extdata", "voloutofbounds.mat",
#                           package = "Dmri3d")
#    suppressPackageStartupMessages(library(spatstat))
#    suppressPackageStartupMessages(library(R.matlab))
#    suppressPackageStartupMessages(library(assertthat))
#    dmin = NULL
#  }
  if (!file.exists(filename))
    stop(paste(filename, "does not exist"))
  #
  if (!grep("\\.mat$", filename))
    stop(paste(filename, "must have extension <mat>"))
  mat = readMat(filename, drop = TRUE)
  if (names(mat)[1] != "mri3d.data")
    stop(paste(filename, "is not a valid mri3d data file"))
  mat = mat$mri3d.data
  #unknownStruct = paste(basename(filename), " has an currently unhandled structure")
  exeVersion = try(unlist(mat["version", ,]), silent = TRUE)
  # Bug in matlab  (jumps into sdata)
  if (class(exeVersion) == "try-error") {
    # no version info
    exeVersion = NA
  } else if (is.null(exeVersion)) {
    exeVersion = unlist(mat[4, ,]) # try mod
    if (!is.character(exeVersion))
      exeVersion = unlist(mat[3, ,]) # try sdata
    if (!is.character(exeVersion))
      exeVersion = NA
  }
  
  expectedNames = c("img", "var", "sdata")
  matNames = dimnames(mat)[[1]]
  if (!all(expectedNames %in% matNames))
    stop("Required top level names ", expectedNames, " not found:", matNames)
  
  img = mat["img", 1, 1][[1]]
  var = mat["var", 1, 1][[1]][,,1]
  sdata = mat["sdata", 1, 1][[1]]
  if (is.character(sdata) || length(sdata) == 0) {
    # we do a radical search for sdata
    for (i in 1:length(var)) {
      dn = unlist(dimnames(var[[i]]))[1]
      if (!is.null(dn) && dn == "click.pt") {
        sdata = var[[i]]
        #warning(filename, " used radical search for structure <sdata>")
        break
      }
    }
  }
  
  # Cleanup unused: we only use sdata, img and var
  rm(mat)
  # Process image
  assert_that(dim(img)[1] >= 10) # Expected size, might change, just for sure
  assert_that(dim(img)[2] == 1)
  nSlices = dim(img)[3]
  xy = list()
  area = list()
  
  # Process Slices
  for (i in 1:nSlices) {
    slice = img[, 1, i]
    nsegments = slice$segment[1, 1]
    if (nsegments == 0)
      next # skip empty
    for (segment in 1:nsegments) {
      idx = 1:slice$scon.count[segment]
      if (length(idx) < 3)
        next
      group = slice$scon.group[1, segment]
      x = slice$scon.x[segment, idx]
      y = slice$scon.y[segment, idx]
      xy0 = try(simplifyContour(x, y, dmin), silent = TRUE)
      if (inherits(xy0, "try-error")) {
        #write.csv(data.frame(x=x,y=y),file=paste0("slice",i,"_",segment,".csv"),
        #          row.names=FALSE)
        stop(
          paste(
            "simplifyContour failed in slice", i, ", segment",
            segment, ": ", attr(xy0, "condition")$message
          )
        )
      }
      xy0$contour$segment = segment
      xy0$contour$slice = i
      xy0$contour$group = group
      xy = c(xy, list(xy0$contour)) # casting to list is required
      area = c(area,
               list(
                 data.frame(
                   area = xy0$area, segment = segment,
                   slice = i, group = group
                 )
               ))
    }
  }
  xy = do.call(rbind, xy)[,c(4, 3, 5, 1, 2)] # Slice vertices
  area = do.call(rbind, area)[,c(3, 2, 4, 1)] # Slice area
  
  # Extract parameters
  vars = lapply(var, function(x) {
    if (is.matrix(x) && dim(x)[1] == 1) {
      as.vector(x)
    } else
      x
  })
  vars$group.names = unlist(vars$group.names)
  vars$examination.date = strptime(vars$examination.date, "%Y.%m.%d / %H:%M:%S")
  vars$file.date = file.info(filename)$mtime
  vars$exeVersion = exeVersion
  
  indices = NULL
  vertices = NULL
  d3 = NULL
  # we do a radical search for d3
  for (i in 1:length(var)) {
    dn = unlist(dimnames(var[[i]]))[1]
    if (!is.null(dn) && dn == "flag") {
      d3 = var[[i]]
      break
    }
  }
  if (is.null(d3)) {
    warning("File <", basename(filename),
            "> has no 3D data.\n", dimnames(d3[[1]]))
  } else {
    #d3names = dimnames(d3)[[1]]
    assert_that(dim(d3)[1] >= 16) # Check if size of structure changed
    assert_that(dim(d3)[2] == 1)
    #assert_that(dim(d3)[3] <= length(vars$group.names))
    nVolumes = max(dim(d3)[3], length(vars$group.names))
    vol3d = NULL
    area3d = NULL
    group.names = NULL
    vertices = list()
    indices = list()
    for (vol in 1:nVolumes) {
      area3d0 = NULL
      group.name = vars$group.names[vol]
      if (is.na(group.name))
        next
      # Remove duplicates
      if (group.name %in% group.names) {
        warning(filename, "; Duplicate group name removed: ", group.name)
        next
      }
      vol3d0 = try(unlist(d3["vol", 1, vol]), silent = TRUE)
      if (inherits(vol3d0, "try-error") || is.null(vol3d0)
          || vol3d0 == 0) {
        warning(group.name, " has no volume")
        next
      }
      area3d0 = try(unlist(d3["area", 1, 1]), silent = TRUE)
      if (inherits(vol3d0, "try-error") || is.null(area3d0)
          || area3d0 == 0) {
        warning(group.name, " has no area")
        next
      }
      # for RGL we need the transposed form
      indices0 = d3["f", ,vol]$f
      if (is.null(indices0))
        next
      # do not muck with this strange conversion,
      # the simple solution with x[] = as.integer(x) fails here
      # as required by RGL
      indices0 = t(matrix(as.integer(indices0), nrow = nrow(indices0)))
      vertices0 = d3["v", ,vol]$v
      if (is.null(vertices0))
        next
      vertices0 = t(matrix(as.double(vertices0), nrow = nrow(vertices0)))
      # Append data
      group.names = c(group.names, group.name)
      vol3d = c(vol3d, vol3d0)
      area3d = c(area3d, area3d0)
      vertices[[group.name]] = vertices0
      indices[[group.name]] = indices0
    }
    names(vol3d) = group.names
    vars$vol3d = vol3d
    names(area3d) = group.names
    vars$area3d = area3d
    vars$group.names = group.names
  }
  list(
    vars = vars, slices = xy, slicearea = area,
    sdata = sdata,
    indices = indices, vertices = vertices
  )
}


simplifyContour = function(x,y,dmin = NULL) {
  swap = FALSE
  # Remove duplicates
  xyu = unique(data.frame(x = x, y = y))
  # Not sure about alignment, so try
  ow = try(owin(poly = with(xyu, list(x = x, y = y))), silent = TRUE)
  if (inherits(ow, "try-error")) {
    swap = TRUE
    ow = owin(poly = with(xyu, list(x = y, y = x)))
  }
  if (is.null(dmin)) {
    return(list(area = area.owin(ow), contour = data.frame(x = x, y = y)))
  }
  ows = simplify.owin(ow, dmin)$bdry[[1]]
  if (swap)
    ct = data.frame(x = ows$y, y = ows$x)
  else
    ct = data.frame(x = ows$x, y = ows$y)
  list(area = ows$area, contour = ct)
}