R/ReadMRIXML.R

In dmenne/dmrixml: Zuerich MRI-XML Support and Database Functions

##' @title Reading MRI files in XML Zuerich Format
##' @description
##' \code{ReadMRIXML} reads MRI files in XML format defined by \code{MRIZH.xsd},
##' returning data seperately for each slice. This is a useful format for
##' storing results in a relational database.
##'
##' \code{ReadMRIXMLxyz} reads MRI files in XML format returning voxel values
##' for each x,y,z tuple, for use in displaying 3D.
##'
##' @name ReadMRIXML
##' @aliases ReadMRIXMLxyz
##' @param mriFile the path of the MRI file in XML format
##' @details
##' \itemize{
##'   \item \code{Scan}: patient initials, number, type...
##'   \item \code{Slices}: slice information (z)
##'   \item \code{Voxels}: voxels (x, y, val)
##' }
#'
##' @export ReadMRIXML
##' @export ReadMRIXMLxyz
##' @import XML stringr
##' @author Dieter Menne <dieter.menne@@menne-biomed.de>
##' @examples
##' mriFile = system.file("extdata", "testMRIZH.xml", package = "Dmrixml")
##' mri = ReadMRIXML(mriFile)
##' print(str(mri))
##' mrixyz = ReadMRIXMLxyz(mriFile)
##' print(str(mrixyz))
# -------------------------- ReadMRIXML  --------------------------------------
ReadMRIXML = function(mriFile) {
  doc = xmlParse(mriFile)
  
  Scan = ExtractScan(doc)
  # Data frame of slices
  slices = getNodeSet(doc, "//slice")
  SliceID = 1
  Slices = NULL
  Voxels = NULL
  for (slice in slices) {
    Slice = data.frame(
      SliceID = SliceID,
      z  =  as.numeric(xmlGetAttr(slice, "z")),
      sliceUnit  =  xmlGetAttr(slice, "unit"),
      orientation = xmlGetAttr(slice, "orientation"),
      what = xmlGetAttr(slice, "what")
    )
    Slices = rbind(Slices, Slice)
    # Voxels for current slice
    Voxels = rbind(
      Voxels, data.frame(
        SliceID = SliceID,
        x = as.numeric(unlist(
          getNodeSet(slice, "voxels/voxel/@x")
        )),
        y = as.numeric(unlist(
          getNodeSet(slice, "voxels/voxel/@y")
        )),
        val = as.numeric(unlist(
          getNodeSet(slice, "voxels/voxel/@val")
        ))
      )
    )
    SliceID = SliceID + 1
  }
  list(Scan = Scan, Slices = Slices, Voxels = Voxels)
}

# -------------------------- ReadMRIXMxyz  --------------------------------------
ReadMRIXMLxyz = function(mriFile) {
  doc = xmlTreeParse(mriFile, useInternalNodes = TRUE)
  Scan = ExtractScan(doc)
  # Data frame of slices
  slices = getNodeSet(doc, "//slice")
  SliceID = 1
  Voxels = NULL
  for (slice in slices) {
    z  =  as.numeric(xmlGetAttr(slice, "z"))
    # Voxels for current slice
    Voxels = rbind(
      Voxels, data.frame(
        SliceID = SliceID,
        x = as.numeric(unlist(
          getNodeSet(slice, "voxels/voxel/@x")
        )),
        y = as.numeric(unlist(
          getNodeSet(slice, "voxels/voxel/@y")
        )),
        z = z,
        val = as.numeric(unlist(
          getNodeSet(slice, "voxels/voxel/@val")
        ))
      )
    )
    SliceID = SliceID + 1
  }
  list(Scan = Scan, Voxels = Voxels)
}
##' @name ValidateMRIXML
##' @title Validate MRI XML with xsd Schema
##' @description Uses schema in MRIZH.xsd to test an XML file.
##'
##' @param mriFile the path of the MRI file in XML format
##' @return \code{valid} on success, or a list of error strings on failure
##' @export
##' @examples
##' badMriFile = system.file("extdata", "testBadMRIZH.xml", package = "Dmrixml")
##' ValidateMRIXML(badMriFile)

# -------------------------- ValidateMRIXML  --------------------------------------
ValidateMRIXML = function(mriFile) {
  xsd = xmlParse(system.file("extdata", "MRIZH.xsd", package = "Dmrixml"),
                 isSchema = TRUE)
  doc = xmlInternalTreeParse(mriFile)
  val = xmlSchemaValidate(xsd, doc)
  if (val$status == 0)
    return("valid")
  return(unlist(lapply(val$errors, "[[", "msg")))
}



# -------------------------- ExtractScan  --------------------------------------
ExtractScan <- function(doc) {
  RandNo =  xpathSApply(doc, "/xhrm/study/patient/inStudyPatientID", xmlValue)
  Initials =  xpathSApply(doc, "/xhrm/study/patient/initials", xmlValue)
  # We assume there is only one record in an xml file, there could be more
  RecordDate =  xpathSApply(doc, "/xhrm/study/records/record/recordDate", xmlValue)
  # Assumes only on scan per file
  ScanTime =  xpathSApply(doc, "//scan/time", xmlValue)
  Time = paste(RecordDate, ScanTime)
  resolutionX = xpathSApply(doc, "//scan/resolutionX", xmlValue)
  resolutionY = xpathSApply(doc, "//scan/resolutionY", xmlValue)
  resolutionZ = xpathSApply(doc, "//scan/resolutionZ", xmlValue)
  patientOrientation = xpathSApply(doc, "//scan/patientOrientation", xmlValue)
  mriType = xpathSApply(doc, "//scan/mriType", function(el)
    xmlGetAttr(el, "val"))
  mriUnit = xpathSApply(doc, "//scan/mriType", function(el)
    xmlGetAttr(el, "unit"))
  list(
    RandNo = as.integer(RandNo),
    Initials = Initials,
    RecordDate = as.POSIXct(RecordDate),
    Time = as.POSIXct(Time),
    resolutionX = as.numeric(resolutionX),
    resolutionY = as.numeric(resolutionY),
    resolutionZ = as.numeric(resolutionZ),
    patientOrientation = patientOrientation,
    mriType = mriType,
    mriUnit = mriUnit
  )
}