QSM.jl

julia> ]add QSM

using QSM

# constants
γ = 267.52      # gyromagnetic ratio
B0 = 3          # main magnetic field strength

# load 3D single-, or multi-echo data using your favourite
# package, e.g. MAT.jl, NIfTI.jl, ParXRec.jl, ...
mag, phas = ...

bdir = (...,)   # direction of B-field
vsz  = (...,)   # voxel size
TEs  = [...]    # echo times

# extract brain mask from last echo using FSL's bet
mask0 = bet(@view(mag[:,:,:,end]), vsz, "-m -n -f 0.5")

# erode mask
mask1 = erode_mask(mask0, 5)

# unwrap phase + harmonic background field correction
uphas = unwrap_laplacian(phas, mask1, vsz)

# convert units
@views for t in axes(uphas, 4)
    uphas[:,:,:,t] .*= inv(B0 * γ * TEs[t])
end

# remove non-harmonic background fields
fl, mask2 = vsharp(uphas, mask1, vsz)

# dipole inversion
x = rts(fl, mask2, vsz, bdir=bdir)

julia> QSM.reset_threading()