A small neuroimage interactive plotter

Manipulate Package

The manipulate from RStudio allows you to create simple Tcl/Tk operators for interactive visualization. I will use it for a simple slider to view different slices of an image.

library(manipulate)

fslr package

I'm calling the fslr package because I know that if you have it installed, you will likely have FSL and have a 1mm T1 template from MNI in a specific location. fslr also loads the oro.nifti package so that readNIfTI is accessible after loading fslr. You can download a test NIfTI image here if you don't have access to any and don't have FSL downlaoded.

Here I will read in the template image:

library(fslr)
options(fsl.path='/usr/local/fsl')
template = file.path(fsldir(), "data/standard", 
                     "MNI152_T1_1mm_brain.nii.gz")
img = readNIfTI(template)

The iplot function

The iplot function defined below takes in a nifti object, the specific plane to be plotted and additional options to be passed to oro.nifti::image. The function is located on my GitHub here.

iplot = function(img, plane = c("axial", 
                                "coronal", "sagittal"), ...){
  ## pick the plane
  plane = match.arg(plane, c("axial", 
                             "coronal", "sagittal"))
  # Get the max number of slices in that plane for the slider
  ns=  switch(plane,
              "axial"=dim(img)[3],
              "coronal"=dim(img)[2],
              "sagittal"=dim(img)[1])
  ## run the manipulate command
  manipulate({
    image(img, z = z, plot.type= "single", plane = plane, ...)
    # this will return mouse clicks (future experimental work)
    pos <- manipulatorMouseClick()
    if (!is.null(pos)) {
      print(pos)
    }
  },
  ## make the slider
  z = slider(1, ns, step=1, initial = ceiling(ns/2))
  )
}

Example plots

Here are some examples of how this iplot function would be used:

iplot(img)
iplot(img, plane = "coronal")
iplot(img, plane = "sagittal")

The result will be a plotted image of the slice with a slider. This is most useful if you run it within RStudio.

Below are 2 example outputs of what you see in RStudio:

Slice 91:

Slice 145:

Conclusions

The iplot function allows users to interactively explore neuroimages. The plotting is not as fast as I'd like, I may try to speed up the oro.nifti::image command or implement some subsampling. It does however show a proof of concept how interactive neuroimaging visualization can be done in R.

Note

manipulate must be run in RStudio for manipulation. The fslr function fslview will call FSLView from FSL for interactive visualization. This is an option of interactive neuroimaging “in R”, but not a real or satisfactory implementation for me (even though I use it frequently). If anyone has implemented such a solution in R, I'd love to hear about it.

matlabr: a Package to Calling MATLAB from R with system

In my research, I primarily use R, but I try to use existing code if available. In neuroimaging and other areas, that means calling MATLAB code. There are some existing solutions for the problem of R to MATLAB: namely the R.matlab package and the RMatlab package (which can call R from MATLAB as well). I do not use thse solutions usually though.

Previously, Mandy Mejia wrote “THREE WAYS TO USE MATLAB FROM R”. Option 2 is about how to use R.matlab, and Mandy gives and example with some cod. She also describes in Options 1 and 3 how to use the system command to call MATLAB commands.

I like this strategy options because:

1. I didn’t take the time to learn R.matlab.
2. It worked for me.
3. I wrote a package to wrap the options Mandy described: matlabr.

matlabr: Wrapping together system calls to MATLAB

The matlabr package is located in GitHub and you can install it with the following command:

devtools::install_github("muschellij2/matlabr")

It has a very small set of functions and I will go through each function and describe what they do:

1. get_matlab: Mostly internal command that will return a character string that will be passed to system. If matlab is in your PATH (bash variable), and you are using R based on the terminal, the command would return "matlab". If MATLAB is not in your PATH or using a GUI-based system like RStudio, you must set options(matlab.path='/your/path/to/matlab').
2. have_matlab: Wrapper for get_matlab to return a logical if matlab is found.
3. run_matlab_script: This will pass a .m file to MATLAB. It also wraps the command in a try-catch statement in MATLAB so that if it fails, it will print the error message. Without this try-catch, if MATLAB errors, then running the command will remain in MATLAB and not return to R.
4. run_matlab_code: This takes a character vector of MATLAB code, ends lines with ;, writes it to a temporary .m file, and then runs run_matlab_script on the temporary .m file.
5. rvec_to_matlab: Takes in a numeric R vector and creates a MATLAB column matrix.
6. rvec_to_matlabclist: Takes in a vector from R (usually a character vector) and quotes these strings with single quotes and places them in a MATLAB cell using curly braces: { and }. It then stacks these cells into a “matrix” of cells.

Setting up matlabr

Let’s set up the matlab.path as I’m running in RStudio:

library(matlabr)
options(matlab.path = "/Applications/MATLAB_R2014b.app/bin")
have_matlab()

The result from have_matlab() indicates that the matlab command can be called.

Let’s write some code to test it

Here we will create some code to take a value for x, y, z (scalars) and a matrix named a and then save x, a, z to a text file:

code = c("x = 10", 
         "y=20;", 
         "z=x+y", 
         "a = [1 2 3; 4 5 6; 7 8 10]",
         "save('test.txt', 'x', 'a', 'z', '-ascii')")
res = run_matlab_code(code)

/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpHnOinq/file2f8352c04937.m

Output

First off, we see that test.txt indeed was written to disk.

file.exists("test.txt")

[1] TRUE

We can read in the test.txt from using readLines:

output = readLines(con = "test.txt")
print(output)

[1] "   1.0000000e+01"                                
[2] "   1.0000000e+00   2.0000000e+00   3.0000000e+00"
[3] "   4.0000000e+00   5.0000000e+00   6.0000000e+00"
[4] "   7.0000000e+00   8.0000000e+00   1.0000000e+01"
[5] "   3.0000000e+01"                                

Conclusions

matlabr isn’t fancy and most likely has some drawbacks as using system can have some quirks. However, these functions have been helpful for me to use some SPM routines and other MATLAB commands while remaining “within R“. R.matlab has a better framework, but it may not be as straightforward for batch processing. Also matlabr has some wrappers that will do a try-catch so that you don’t get stuck in MATLAB after calling system.

Let me know if this was helpful or if you have ideas on how to make this better. Or better yet, give a pull request.

White Matter Segmentation in R

Goals and Overall Approach

We will use multiple packages and pieces of software for white matter (and gray matter/cerebro spinal fluid (CSF)) segmentation.

The overall approach will be, with the required packages in parentheses:

1. N4 Inhomogeneity Bias-Field Correction (extrantsr and ANTsR)
2. Brain extraction using BET and additional tools (extrantsr and fslr)
3. FAST for tissue-class segmentation. (fslr)

Installing Packages

Below is a script to install all the current development versions of all packages. The current fslr packages depends on oro.nifti (>= 0.5.0) , which is located at muschellij2/oro.nifti or bjw34032/oro.nifti.

Note, the ITKR and ANTsR packages can take a long time to compile. The extrantsr package builds on ANTsR and makes some convenience wrapper functions.

devtools::install_github("muschellij2/oro.nifti")
devtools::install_github("muschellij2/fslr")
devtools::install_github("stnava/cmaker")
devtools::install_github("stnava/ITKR")
devtools::install_github("stnava/ANTsR")
devtools::install_github("muschellij2/extrantsr")
install.packages("scales")

Load in the packages

Here we will load in the required packages. The scales package is imported just for the alpha function, used below in plotting.

rm(list=ls())
library(fslr)
library(extrantsr)
library(scales)

Specifying FSL path

For fslr to work, FSL must be installed. If run in the Terminal, the FSLDIR environmental variable should be found using R's Sys.getenv("FSLDIR") function.

If run in an IDE (such as RStudio or the R GUI), R must know the path of FSL, as set by the following code:

options(fsl.path="/usr/local/fsl/")

Image Filenames

Here we will set the image name. The nii.stub function will strip off the .nii.gz from img.name.

img.name = "SUBJ0001-01-MPRAGE.nii.gz"
img.stub = nii.stub(img.name)

N4 Bias Field Correction

The first step in most MRI analysis is performing inhomogeneity correction. The extrantsr function bias_correct can perform N3 or N4 bias correction from the ANTsR package.

n4img = bias_correct( img.name, correction = "N4", 
                      outfile = paste0(img.stub, "_N4.nii.gz") )

ortho2(n4img)

Let us note that the image is of the head and a bit of the neck. We wish to perform white matter segmentation only on the brain tissues, so we will do brain extraction.

Brain Extraction

The extrantsr function fslbet_robust performs brain extraction. It relies on the fslr function fslbet which calls bet from FSL. It also performs neck removal (remove.neck = TRUE) and will perform BET once and then estimate a new center of gravity (COG) and then re-run BET. These functions are implemented in fslbet specifically, but these have been re-implemented in fslbet_robust in a slightly different way. fslbet_robust will also perform N4 inhomogeneity correction, but as this has already been performed above, we will set correct = FALSE.

For neck removal, a template brain and mask must be specified. We will use the T1, 1mm resolution, MNI brain included with FSL's installation.

bet = fslbet_robust(img = n4img, 
                    retimg = TRUE,
                    remove.neck = TRUE,
                    robust.mask = FALSE,
                    template.file = file.path( fsldir(), 
                                               "data/standard", 
                                               "MNI152_T1_1mm_brain.nii.gz"),
                    template.mask = file.path( fsldir(), 
                                               "data/standard", 
                                               "MNI152_T1_1mm_brain_mask.nii.gz"), 
                    outfile = "SUBJ0001-01-MPRAGE_N4_BET", 
                    correct = FALSE)

The results look good – the brain tissue is kept (in red) only. Not much brain tissue is discarded nor non-brain-tissue is included.

ortho2(n4img, bet > 0, 
       col.y=alpha("red", 0.5))

FAST Image Segmentation

Now that we have a brain image, we can use FAST for image segmentation. We will use the fslr function fast, which calls fast from FSL. We will pass the -N option so that FAST will not perform inhomogeneity correction (different from N4 and N3), because we had performed this before.

fast = fast(file = bet, 
            outfile = paste0(img.stub, "_BET_FAST"), 
            opts = '-N')

White Matter Results

By default, FAST assumes 3 tissue classes, generally white matter, gray matter, and CSF. These are generally ordered by the mean intensity of the class. For T1-weighted images, white matter is the highest intensity, and assigned class 3. Let's see the results:

ortho2(bet, fast == 3, 
       col.y=alpha("red", 0.5))

Gray Matter / CSF Results

We can also visualize the classes for 1 and 2 for CSF and gray matter, respectively.

ortho2(bet, fast == 1, col.y=alpha("red", 0.5), text="CSF Results")

ortho2(bet, fast == 2, col.y=alpha("red", 0.5), text="Gray Matter\nResults")

The results indicate good segmentation of the T1 image. The fslr function fast result in more than the tissue-class segmentation, see the other files output:

list.files(pattern=paste0(img.stub, "_BET_FAST"))

[1] "SUBJ0001-01-MPRAGE_BET_FAST_mixeltype.nii.gz"
[2] "SUBJ0001-01-MPRAGE_BET_FAST_pve_0.nii.gz"    
[3] "SUBJ0001-01-MPRAGE_BET_FAST_pve_1.nii.gz"    
[4] "SUBJ0001-01-MPRAGE_BET_FAST_pve_2.nii.gz"    
[5] "SUBJ0001-01-MPRAGE_BET_FAST_pveseg.nii.gz"   
[6] "SUBJ0001-01-MPRAGE_BET_FAST_seg.nii.gz"      

Conclusions

It's a exciting time to be working in neuroimaging in R. The fslr and ANTsR packages provide functionality to perform operations for neuroimaging processing. I will be doing a series on some of the options for analysis in the coming weeks. The code for this analysis (and the data) is located at https://github.com/muschellij2/HopStat/blob/gh-pages/White_Matter_Segmentation_in_R/

Notes

The fslr function ortho2 is a rewrite of the oro.nifti::orthographic function, but with different defaults and will set values of 0 in the second image (y argument) to NA.