MRIRealign.jl

Documentation	Build Status	Test Coverage

MRIRealign.jl performs rigid-body (6-DOF) motion correction for 4-D MRI time-series data. It estimates three rotation angles and three translations per volume by minimizing the sum of squared intensity differences, then reslices (resamples) the volumes to undo the estimated motion.

The algorithm follows the seminal paper by Friston et al. and was heavily inspired by SPM's spm_realign function. Key differences from SPM include:

• Speed — a Gauss–Newton trust-region optimizer with exact analytic Jacobians of the rotation matrix converges in fewer iterations than SPM's re-estimation loop.
• Consensus estimation — all time frames are aligned pairwise and a robust weighted consensus is computed via iteratively reweighted least squares with geodesic rotation distance on SO(3). This is less sensitive to the image quality of any single reference frame and avoids difficulties in mapping to a blurred temporal mean.

Quick Start

On Unix systems, Julia can be installed with

curl -fsSL https://install.julialang.org | sh

and on Windows systems with

winget install --name Julia --id 9NJNWW8PVKMN -e -s msstore

More detailed installation instructions can be found here.

Thereafter, you can start Julia from the command line with

julia

Loading data

This tutorial assumes that you have a folder at the path /path_to_files/ with NIfTI files of the format mask.nii and somename_1.nii, somename_2.nii, … . MRIRealign.jl does not include I/O functions, so you are free to load data from NIfTI, DICOM, MAT, HDF5 files, etc.

Install the packages once:

using Pkg
Pkg.add("MRIRealign")
Pkg.add("NIfTI")

Then load them:

using MRIRealign
using NIfTI

Change to the data directory:

cd("/path_to_files/")

Optionally, load a mask and convert it to a BitArray:

mask = round.(Bool, niread("mask.nii"))

Create a sorted list of volume file names (natural numeric order):

files = filter(f -> isfile(f) && f != "mask.nii", readdir())
files = sort(files, by = file -> parse(Int, match(r"\d+", file).match))

Allocate a 4-D array and read all volumes into it:

img = Array{Float64}(undef, size(mask)..., length(files))

for t in eachindex(files)
    img[:,:,:,t] .= niread(files[t]).raw
end

Estimating and applying motion correction

Call realign!, which overwrites img with the aligned volumes and returns the motion parameters — a (6, t) matrix where each column is [rx, ry, rz, tx, ty, tz] (rotations in radians, translations in voxels):

params = realign!(img; mask=mask)

Write the aligned images back to NIfTI files:

for t in eachindex(files)
    ni = niread(files[t])
    ni.raw .= img[:,:,:,t]
    niwrite(files[t], ni)
end

Estimate-only workflow

To estimate motion parameters without modifying the images:

params = realign!(img; mask=mask, realign=false)

The returned params can later be applied with the two-argument form:

realign!(img, params)

Note: This tutorial assumes that the NIfTI headers of all files are identical and replaces the raw data with interpolated data. To update the NIfTI header instead (preserving the original voxel data), use realign=false and write the parameters into the header. See the NIfTI.jl documentation for details.