MriResearchTools

NumART2star(image::AbstractArray{T,4}, TEs) where T

Performs T2* calculation on 4D-multi-echo magnitude data. https://doi.org/10.1002/mrm.10283

RSS(mag; dims=ndims(mag))

Performs root-sum-of-squares combination along the last dimension of mag. The dimension can be specificed via the dims keyword argument.

brain_mask(mask)

Tries to extract the brain from a mask with skull and a gap between brain and skull.

Examples

julia> mask = robustmask(mag)
julia> brain = brain_mask(mask)

See also robustmask, ROMEO.romeovoxelquality, phase_based_mask

calculateB0_unwrapped(unwrapped_phase, mag, TEs)

Calculates B0 in [Hz] from unwrapped phase. TEs in [ms]. The phase offsets have to be removed prior.

See also mcpc3ds and romeo

estimatenoise(image::AbstractArray)

Estimates the noise from the corners of the image. It assumes that at least one corner is without signal and only contains noise.

estimatequantile(array, p)

Quickly estimates the quantile p of a possibly large array by using a subset of the data.

gaussiansmooth3d(image)

gaussiansmooth3d(image, sigma=[5,5,5];
    mask=nothing,
    nbox=ifelse(isnothing(mask), 3, 6), 
    weight=nothing, dims=1:min(ndims(image),3), 
    boxsizes=getboxsizes.(sigma, nbox)
    )

Performs Gaussian smoothing on image with sigma as standard deviation of the Gaussian. By application of nbox times running average filters in each dimension. The length of sigma and the length of the dims that are smoothed have to match. (Default 3)

Optional arguments:

• mask: Smoothing can be performed using a mask to inter-/extrapolate missing values.
• nbox: Number of box applications. Default is 3 for normal smoothing and 6 for masked smoothing.
• weight: Apply weighted smoothing. Either weighted or masked smoothing can be porformed.
• dims: Specify which dims should be smoothed. Corresponds to manually looping of the other dimensions.
• boxizes: Manually specify the boxsizes, not using the provided sigma. length(boxsizes)==length(dims) && length(boxsizes[1])==nbox

gaussiansmooth3d(image)

gaussiansmooth3d(image, sigma=[5,5,5];
    mask=nothing,
    nbox=ifelse(isnothing(mask), 3, 6), 
    weight=nothing, dims=1:min(ndims(image),3), 
    boxsizes=getboxsizes.(sigma, nbox)
    )

Performs Gaussian smoothing on image with sigma as standard deviation of the Gaussian. By application of nbox times running average filters in each dimension. The length of sigma and the length of the dims that are smoothed have to match. (Default 3)

Optional arguments:

• mask: Smoothing can be performed using a mask to inter-/extrapolate missing values.
• nbox: Number of box applications. Default is 3 for normal smoothing and 6 for masked smoothing.
• weight: Apply weighted smoothing. Either weighted or masked smoothing can be porformed.
• dims: Specify which dims should be smoothed. Corresponds to manually looping of the other dimensions.
• boxizes: Manually specify the boxsizes, not using the provided sigma. length(boxsizes)==length(dims) && length(boxsizes[1])==nbox

gaussiansmooth3d_phase(phase, sigma=[5,5,5]; weight=1, kwargs...)

Smoothes the phase via complex smoothing. A weighting image can be given. The same keyword arguments are supported as in gaussiansmooth3d:

gaussiansmooth3d(image)

gaussiansmooth3d(image, sigma=[5,5,5];
    mask=nothing,
    nbox=ifelse(isnothing(mask), 3, 6), 
    weight=nothing, dims=1:min(ndims(image),3), 
    boxsizes=getboxsizes.(sigma, nbox)
    )

Performs Gaussian smoothing on image with sigma as standard deviation of the Gaussian. By application of nbox times running average filters in each dimension. The length of sigma and the length of the dims that are smoothed have to match. (Default 3)

Optional arguments:

• mask: Smoothing can be performed using a mask to inter-/extrapolate missing values.
• nbox: Number of box applications. Default is 3 for normal smoothing and 6 for masked smoothing.
• weight: Apply weighted smoothing. Either weighted or masked smoothing can be porformed.
• dims: Specify which dims should be smoothed. Corresponds to manually looping of the other dimensions.
• boxizes: Manually specify the boxsizes, not using the provided sigma. length(boxsizes)==length(dims) && length(boxsizes[1])==nbox

getHIP(mag, phase; echoes=[1,2])

getHIP(compl; echoes=[1,2])

Calculates the Hermitian Inner Product between the specified echoes.

getVSM(B0, rbw, dim, threshold=5.0)

Calculates a voxel-shift-map. B0 is given in [Hz]. rbw is the receiverbandwidth or PixelBandwidth in [Hz]. dim is the dimension of the readout (in which the distortion occurs)

getsensitivity(mag; sigma, nbox=15)

getsensitivity(mag, pixdim; sigma_mm=7, nbox=15)

getsensitivity(mag::NIVolume, datatype=eltype(mag); sigma_mm=7, nbox=15)

Calculates the bias field using the boxsegment approach. It assumes that there is a "main tissue" that is present in most areas of the object. Published in CLEAR-SWI.

See also makehomogeneous

header(v::NIVolume)

Returns a copy of the header with the orientation information.

Examples

julia> vol = readmag("image.nii")
julia> hdr = header(vol)
julia> savenii(vol .+ 10, "vol10.nii"; header=hdr)

homodyne(mag, phase)

homodyne(mag, phase; dims, sigma)

homodyne(I; kwargs...)

Performs homodyne filtering via division of complex complex smoothing.

homodyne(mag, phase)

homodyne(mag, phase; dims, sigma)

homodyne(I; kwargs...)

Performs homodyne filtering via division of complex complex smoothing.

laplacianunwrap(ϕ::AbstractArray)

Requires using ImageFiltering. Performs laplacian unwrapping on the input phase. (1D - 4D) The phase has to be scaled to radians. The implementation is close to the original publication: Schofield and Zhu 2003, https://doi.org/10.1364/OL.28.001194. It is not the fastest implementation of laplacian unwrapping (doesn't use discrete laplacian).

laplacianunwrap(ϕ::AbstractArray)

Requires using ImageFiltering. Performs laplacian unwrapping on the input phase. (1D - 4D) The phase has to be scaled to radians. The implementation is close to the original publication: Schofield and Zhu 2003, https://doi.org/10.1364/OL.28.001194. It is not the fastest implementation of laplacian unwrapping (doesn't use discrete laplacian).

makehomogeneous(mag::NIVolume; sigma_mm=7, nbox=15)

Homogeneity correction for NIVolume from NIfTI files.

Keyword arguments:

• sigma_mm: sigma size for smoothing to obtain bias field. Takes NIfTI voxel size into account
• nbox: Number of boxes in each dimension for the box-segmentation step.

makehomogeneous(mag; sigma, nbox=15)

Homogeneity correction of 3D arrays. 4D volumes are corrected using the first 3D volume to obtain the bias field.

Keyword arguments:

• sigma: sigma size in voxel for smoothing to obtain bias field. (mandatory)
• nbox: Number of boxes in each dimension for the box-segmentation step.

Larger sigma-values make the bias field smoother, but might not be able to catch the inhomogeneity. Smaller values can catch fast varying inhomogeneities but new inhomogeneities might be created. The stronger the bias field, the more boxes are required for segmentation. With too many boxes, it can happen that big darker structures are captured and appear overbrightened.

Calculates the bias field using the boxsegment approach. It assumes that there is a "main tissue" that is present in most areas of the object. Published in CLEAR-SWI.

See also getsensitivity

makehomogeneous(mag; sigma, nbox=15)

Homogeneity correction of 3D arrays. 4D volumes are corrected using the first 3D volume to obtain the bias field.

Keyword arguments:

• sigma: sigma size in voxel for smoothing to obtain bias field. (mandatory)
• nbox: Number of boxes in each dimension for the box-segmentation step.

Larger sigma-values make the bias field smoother, but might not be able to catch the inhomogeneity. Smaller values can catch fast varying inhomogeneities but new inhomogeneities might be created. The stronger the bias field, the more boxes are required for segmentation. With too many boxes, it can happen that big darker structures are captured and appear overbrightened.

Calculates the bias field using the boxsegment approach. It assumes that there is a "main tissue" that is present in most areas of the object. Published in CLEAR-SWI.

See also getsensitivity

mask_from_voxelquality(qmap::AbstractArray, threshold=:auto)

Creates a mask from a quality map. Another option is to use robustmask(qmap)

Examples

julia> qmap = romeovoxelquality(phase_3echo; TEs=[1,2,3]);
julia> mask = mask_from_voxelquality(qmap);

See also ROMEO.romeovoxelquality, ROMEO.romeo, robustmask, brain_mask

mcpc3ds(phase, mag; TEs, keyargs...)

mcpc3ds(compl; TEs, keyargs...)

mcpc3ds(phase; TEs, keyargs...)

Perform MCPC-3D-S coil combination and phase offset removal on 4D (multi-echo) and 5D (multi-echo, uncombined) input.

Optional Keyword Arguments

• echoes: only use the defined echoes. default: echoes=[1,2]
• sigma: smoothing parameter for phase offsets. default: sigma=[10,10,5]
• bipolar_correction: removes linear phase artefact. default: bipolar_correction=false
• po: phase offsets are stored in this array. Can be used to retrieve phase offsets or work with memory mapping.

Examples

julia> phase = readphase("phase5D.nii")
julia> mag = readmag("mag5D.nii")
julia> combined = mcpc3ds(phase, mag; TEs=[4,8,12])

For very large files that don't fit into memory, the uncombined data can be processed with memory mapped to disk:

julia> phase = readphase("phase5D.nii"; mmap=true)
julia> mag = readmag("mag5D.nii"; mmap=true)
julia> po_size = (size(phase)[1:3]..., size(phase,5))
julia> po = write_emptynii(po_size, "po.nii")
julia> combined = mcpc3ds(phase, mag; TEs=[4,8,12], po)

mcpc3ds_meepi(phase, mag; TEs, keyargs...)

mcpc3ds_meepi(compl; TEs, keyargs...)

mcpc3ds_meepi(phase; TEs, keyargs...)

Perform MCPC-3D-S phase offset removal on 5D MEEPI (multi-echo, multi-timepoint) input. The phase offsets are calculated for the template timepoint and removed from all volumes.

Optional Keyword Arguments

• echoes: only use the defined echoes. default: echoes=[1,2]
• sigma: smoothing parameter for phase offsets. default: sigma=[10,10,5]
• po: phase offsets are stored in this array. Can be used to retrieve phase offsets or work with memory mapping.
• template_tp: timepoint for the template calculation. default: template_tp=1

r2s_from_t2s(t2s) = 1000 ./ t2s

Converts from T2* [ms] to R2* [1/s] values.

readmag(filename; rescale=false, keyargs...)

Reads the NIfTI magnitude with sanity checking and optional rescaling to [0;1].

Examples

julia> mag = readmag("mag.nii")

Optional keyargs are forwarded to niread:

niread(file; mmap=false, mode="r")

Read a NIfTI file to a NIVolume. Set mmap=true to memory map the volume.

readphase(filename; rescale=true, keyargs...)

Reads the NIfTI phase with sanity checking and optional rescaling to [-π;π].

Examples

julia> phase = readphase("phase.nii")

Optional keyargs are forwarded to niread:

niread(file; mmap=false, mode="r")

Read a NIfTI file to a NIVolume. Set mmap=true to memory map the volume.

robustmask(weight::AbstractArray; factor=1, threshold=nothing)

Creates a mask from an intensity/weight images by estimating a threshold and hole filling. It assumes that at least one corner is without signal and only contains noise. The automatic threshold is multiplied with factor.

Examples

julia> mask1 = robustmask(mag); # Using magnitude
julia> mask2 = phase_based_mask(phase); # Using phase
julia> mask3 = robustmask(romeovoxelquality(phase; mag)); # Using magnitude and phase
julia> brain = brain_mask(robustmask(romeovoxelquality(phase; mag); threshold=0.9));

See also ROMEO.romeovoxelquality, phase_based_mask, brain_mask

robustmask(weight::AbstractArray; factor=1, threshold=nothing)

Creates a mask from an intensity/weight images by estimating a threshold and hole filling. It assumes that at least one corner is without signal and only contains noise. The automatic threshold is multiplied with factor.

Examples

julia> mask1 = robustmask(mag); # Using magnitude
julia> mask2 = phase_based_mask(phase); # Using phase
julia> mask3 = robustmask(romeovoxelquality(phase; mag)); # Using magnitude and phase
julia> brain = brain_mask(robustmask(romeovoxelquality(phase; mag); threshold=0.9));

See also ROMEO.romeovoxelquality, phase_based_mask, brain_mask

robustrescale(array, newmin, newmax; threshold=false, mask=trues(size(array)), datatype=Float64)

Rescales the image to the the new range, disregarding outliers. Only values inside mask are used for estimating the rescaling option

savenii(image::AbstractArray, filepath; header=nothing)

savenii(image::AbstractArray, name, writedir, header=nothing)

Warning: MRIcro can only open images with types Int32, Int64, Float32, Float64

Examples

julia> savenii(ones(64,64,5), "image.nii")

julia> savenii(ones(64,64,5), "image2", "folder")

to_dim(V::AbstractVector, dim::Int)

to_dim(a::Real, dim::Int)

Converts a vector or number to a higher dimension.

Examples

julia> to_dim(5, 3)
1×1×1 Array{Int64, 3}:
[:, :, 1] =
 5
julia> to_dim([1,2], 2)
 1×2 Matrix{Int64}:
  1  2

write_emptynii(size, path; datatype=Float32, header=NIVolume(zeros(datatype, 1)).header)

Writes an empty NIfTI image to disk that can be used for memory-mapped access.

Examples

julia> vol = write_emptynii((64,64,64), "empty.nii")
julia> vol.raw[:,:,1] .= ones(64,64) # synchronizes mmapped file on disk

Warning: MRIcro can only open images with types Int32, Int64, Float32, Float64