#' @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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.