MRIRealign.jl

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)

Reference modes

# Robust consensus across all pairwise alignments (default, slowest)
params = realign!(img; ref_mode=:consensus)

# Align to the temporal mean (fast, may be blurred)
params = realign!(img; ref_mode=:mean)

# Align to a specific time frame (fast, quality depends on that frame)
params = realign!(img; ref_mode=1)

Smoothing

For noisy data, applying Gaussian smoothing before estimation can improve robustness. The fwhm keyword accepts a 3-tuple of full-width-at-half-maximum values in voxel units:

params = realign!(img; fwhm=(5.0, 5.0, 5.0))

Note

The default fwhm=nothing (no smoothing) differs from SPM's default of approximately 5 mm. For noisy data, setting fwhm explicitly is recommended.

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.

API Reference

Main interface

MRIRealign.realign! — Function

realign!(img; center, ref_mode, mask, fwhm, realign, voxel_size, radius, x_abstol) -> motion_params

Estimate rigid-body (6-DOF) motion parameters from a 4-D MRI time series and, optionally, reslice the volumes to undo the estimated motion.

The algorithm minimizes the sum of squared intensity differences between each volume and a reference, using a Gauss–Newton trust-region optimizer with exact analytic Jacobians of the rotation matrix.

Arguments

img::AbstractArray{T,4}: image array with dimensions (x, y, z, t). T may be real- or complex-valued. For complex data the motion estimation is performed on the magnitude; the phase is resliced along with the magnitude when realign=true.

Keyword arguments

center=size(img)[1:3] .÷ 2: rotation center in voxel coordinates. Only affects the parameterization of the motion (i.e., the returned translation values); the aligned images are identical regardless of center.
ref_mode=:consensus: reference strategy. One of
- :consensus — align every pair of time frames and compute a robust weighted consensus (IRLS with geodesic rotation distance). Most accurate, but t-fold slower than a single reference.
- :mean — align each frame to the temporal mean. Fast, but the mean may be blurred when motion is large.
- an Integer — align every frame to the given time-frame index. Fast; works well when that frame has good image quality.
mask=trues(size(img)[1:3]): BitArray or Bool array selecting the voxels over which the cost function is evaluated.
fwhm=nothing: optional Gaussian smoothing kernel specified as a 3-tuple of full-width-at-half-maximum values (σx, σy, σz) in voxel units. nothing (default) applies no smoothing. Smoothing can improve robustness for noisy data.
realign=true: if true, img is overwritten in-place with the motion-corrected volumes. If false, only the parameters are estimated.
voxel_size::NTuple{3}=(1.0, 1.0, 1.0): voxel dimensions in millimeters, e.g. (1.5, 1.5, 3.0). Used to define the interpolant on mm-spaced axes so that translations, the affine matrix, and all spatial derivatives are natively in mm.
radius::Real=64.0: characteristic head radius in millimeters. Rotation angles (in radians) are multiplied by radius to obtain arc-length displacements in mm, placing rotations on the same footing as translations inside the optimizer.
x_abstol::Real=1e-3: absolute convergence tolerance for the optimizer, in millimeters. The optimizer terminates when the parameter step falls below this value. The default of 1e-3 mm (1 μm) is far below any practical MRI resolution.

Returns

motion_params::Matrix{T} of size (6, t). Each column holds [rx, ry, rz, tx, ty, tz] — three rotation angles in radians and three translations in millimeters — for the corresponding time frame.

Note that voxel_size, radius, and x_abstol can in any unit of translation, as long as all three are the same. The translations in the return motion_params are in the units of the input parameters.

Examples

# Estimate and apply motion correction (1 mm isotropic voxels)
params = realign!(img)

# Specify voxel size and head radius
params = realign!(img; voxel_size=(1.5, 1.5, 3.0), radius=70.0)

# Estimate only (no reslicing), using frame 1 as reference
params = realign!(img; ref_mode=1, realign=false)

# With a brain mask, smoothing, and a coarser tolerance (10 μm)
params = realign!(img; mask=brain_mask, fwhm=(5.0, 5.0, 5.0), x_abstol=0.01)

realign!(img, motion_params; center, voxel_size) -> motion_params

Reslice img in-place using pre-computed motion_params (e.g., from a previous call to realign!).

Arguments

img::AbstractArray{T,4}: image array with dimensions (x, y, z, t).
motion_params::AbstractMatrix: (6, t) matrix of motion parameters in the format [rx, ry, rz, tx, ty, tz] per column, with rotations in radians and translations in millimeters.

Keyword arguments

center=size(img)[1:3] .÷ 2: rotation center that was used when motion_params was estimated.
voxel_size::NTuple{3}=(1.0,1.0,1.0): voxel dimensions in mm.

Returns

motion_params (the same matrix that was passed in).

Examples

# Two-step workflow: estimate, then apply later
params = realign!(img; realign=false)
# ... inspect params ...
realign!(img, params)

source

Geometric utilities

MRIRealign.create_rotation_matrix — Function

create_rotation_matrix(rx, ry, rz) -> SMatrix{3,3}
create_rotation_matrix(p)          -> SMatrix{3,3}

Build a 3×3 rotation matrix from three Euler angles rx, ry, rz (in radians), using the ZYX intrinsic convention: R = Rz * Ry * Rx.

The single-argument form extracts (p[1], p[2], p[3]), so any indexable collection (vector, tuple, SVector, …) works.

Examples

julia> R = create_rotation_matrix(0.0, 0.0, 0.0)
3×3 StaticArraysCore.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0