In neuroimaginador/wf4ni: Work-Flows for Neuroimage Processing

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = nzchar(Sys.getenv("MY_EVAL_VARIABLE"))
)

Objective

In this vignette, we're going to develop a pipeline (workflow) to build the deterministic tractography from a diffusion image, using the wf4ni package.

The main steps in this pipeline will be:

• Eddy current correction.
• Tensor computation.
• Fiber tracking inside a mask.
• Computation of Fractional Anisotropy (FA) map.

To demonstrate the flexibility of wf4ni, we'll use different R packages to perform several of the above steps.

It must be noted that the philosophy of wf4ni is to separate the pipeline logic from its implementation. So, the same workflow could be used just changing, if needed, the implementation of any (or all) of the steps above.

The organization of this vignette is as follows. First, we'll describe the dependencies or packages needed to run the workflow. Later, we'll describe the workflow, step by step, without presenting the actual implementation of the functions used. Then, we'll use some sample data to run an experiment with the workflow, presenting its results. We'll end the vignette with an annex in which we'll present the actual implementation of the functions.

Requirements

We consider two types of requirements: packages needed to perform the key steps in the workflow and data which may be necessary at any step of the computation.

Packages

First, to build the flow, the wf4ni package is needed:

library(wf4ni)

The following packages are needed in some points of the pipeline:

library(tidyverse)
library(neurobase)
library(oro.nifti)
library(fslr)
library(dti)
library(readr)

fslr and dti will be used for the processing, while neurobase, oro.nifti and readr will be used for reading and writing files. The tidyverse is useful to simplify some code and make it more legible.

Other Resources

In this case, we are not using any other external data for this pipeline.

DTI Processing

The diffusion workflow will be a NIflow object. To create such object, one needs to provide a name, a work_dir where to store temporary results, and the name of the inputs needed to execute the workflow. In this case, our pipeline will need:

• A DWI image in nifti format.
• A matrix of gradient directions with which the DWI image was obtained.
• A vector of b-values for each acquisition in the DWI image.

dti_flow <- NIflow$new(name = "dti",
                       work_dir = "~/flow/dti",
                       inputs = c("dwi", "grads", "bvals"))

Eddy Current Correction

The first step will be to correct the effect of eddy currents in the diffusion acquisition. We'll use the fslr package to this end.

correct_dwi <- function(image) {
  #...
}

This function is added to the flow a follows:

dti_flow %>%
  add(what = correct_dwi,
      inputs = "dwi",
      output = "dwi_corrected")

And a graph representation of this partial flow:

dti_flow$plot()

Tensor Computation

The next step is to incorporate the gradient directions and the b-values in an object suitable to perform the tensor computation by the dti package.

read_data <- function(image, grads, bvals) {
  #...
}

compute_tensors <- function(dwi_data) {
  #...
}

And the way to include these functions in the flow is:

dti_flow %>%
  add(what = read_data,
      inputs = c("dwi_corrected", "grads", "bvals"),
      output = "dwiData") %>%
  add(what = compute_tensors,
      inputs = "dwiData",
      output = "tensors")

The partial flow is:

dti_flow$plot()

Mask Computation

Using tensor information, we can extract the S0 image (the average of the images acquired without gradient) and the mask (using the fslr package):

get_S0 <- function(tensors) {
  #...
}

compute_mask <- function(S0) {
  #...
}

Next, we add these steps to the flow:

dti_flow %>%
  add(what = get_S0,
      inputs = "tensors",
      output = "S0") %>%
  add(what = compute_mask,
      inputs = "S0",
      output = "mask")

dti_flow$plot()

Fiber Tracking

One of the main results of this pipeline is to obtain the fibers associated to the tensors previously computed. We use the following function (a wrapper to a function in the dti package) to perform fiber tracking on the mask already calculated.

fibertracking <- function(tensors, mask) {
  #...
}

And we add this step to the flow:

dti_flow %>%
  add(what = fibertracking,
      inputs = c("tensors", "mask"),
      output = "fibers")

dti_flow$plot()

FA Map

As a side effect of the previous computations, we can compute a FA map from the tensor map (using a mask to remove non-informative voxels outside the brain). This function is also a wrapper to a function in package dti:

compute_FA <- function(tensors, mask) {
  #...
}

And, finally, we add this last step to the flow:

dti_flow %>%
  add(what = compute_FA,
      inputs = c("tensors", "mask"),
      output = "FA_map")

Final Flow

With these final additions, the complete pipeline is as follows:

dti_flow$plot()

We can also obtain a brief summary of the flow by typing:

summary(dti_flow)

Example

To illustrate how this workflow can be used, let us suppose that we have all the functions mentioned before already implemented (an example implementation is given in the Annex), and build the workflow as in the previous steps. We'll apply this flow to compute the fibers and the FA map given a DWI image.

Code

lag <- stats::lag

We'll use sample data for the DWI image (from the Human Connectome Project), and the b-values vector and gradient matrix.

dwi_file <- "./resources/dti/dti_example.nii.gz"
bvals_file <- "./resources/dti/dti_bvals.csv"
grads_file <- "./resources/dti/dti_gradients.csv"

dwi <- readnii(dwi_file)
bvals <- read_csv(bvals_file,
                  col_names = FALSE,
                  col_types = cols()) %>% unlist()
gradients <- read_csv(grads_file,
                      col_names = FALSE,
                      col_types = cols()) %>%
  as.matrix()

res <- dti_flow$execute(inputs = list(dwi = dwi,
                                      grads = gradients,
                                      bvals = bvals),
                        desired_outputs = c("S0", 
                                            "fibers", 
                                            "FA_map"))

Images

Using the function ortho2 in package neurobase, we can inspect both the S0 image and the obtained FA map.

neurobase::ortho2(res$S0)

knitr::include_graphics(path = "./img/dti_S0.png")

neurobase::ortho2(res$FA_map)

knitr::include_graphics(path = "./img/dti_FA_map.png")

Also, we can inspect the computed fiber bundles by using the native function show3d of package dti. To illustrate the results, we present some snapshots of the obtained fibers:

show3d(res$fibers)

knitr::include_graphics(path = "./img/dti_fibers1.png")
knitr::include_graphics(path = "./img/dti_fibers2.png")
knitr::include_graphics(path = "./img/dti_fibers3.png")
knitr::include_graphics(path = "./img/dti_fibers4.png")

Full Flow Log

The NIflow object stores internally a log of all the processes that it goes through and the memory used to store its internal results:

dti_flow$print_log()

lines <- readLines(con = "./resources/dti_log.txt")

cat(lines, sep = "\n")

Annex: Full Code

In this annex, we provide example code for the tractography workflow, which can be loaded and used to replicate the workflow and results given in this vignette.

correct_dwi <- function(image) {

  my_dwi_file <- tempimg(as.nifti(image))

  dwi_corrected <- eddy_correct(infile = my_dwi_file)

  return(dwi_corrected)

}

read_data <- function(image, grads, bvals) {

  my_dwi_corrected_file <- tempimg(image)

  dwi_data <- readDWIdata(gradient = grads,
                          dirlist = my_dwi_corrected_file,
                          bvalue = bvals,
                          format = "NIFTI")

  return(dwi_data)

}

compute_tensors <- function(dwi_data) {

  tensors <- dtiTensor(dwi_data)

  return(tensors)

}

get_S0 <- function(tensors) {

  return(as.nifti(tensors@th0))

}

compute_mask <- function(S0) {

  S0_file <- tempimg(S0)

  brain_mask <- fslbet(infile = S0_file)

  return(brain_mask > 0)

}

fibertracking <- function(tensors, mask) {

  fibers <- tracking(tensors,
                     mask = mask)

  return(fibers)

}

compute_FA <- function(tensors, mask) {

  indices <- dtiIndices(tensors)

  FA_map <- indices@fa * mask

  return(FA_map)

}

neuroimaginador/wf4ni documentation built on Oct. 15, 2019, 9:26 a.m.