MIRT is the Michigan Image Reconstruction Toolbox

LineSearchMM{...}

Mutable struct for MM-based line searches.

LineSearchMM(gradf, curvf, u, v; α0 ...)
LineSearchMM(u, v, dot_gradf, dot_curvf; α0 ...)

Construct iterator for line-search based on majorize-minimize (MM) approach for a general family of 1D cost functions of the form $h(α) = \sum_{j=1}^J f_j(u_j + α v_j)$ where each function $f_j(t)$ has a quadratic majorizer of the form

\[q_j(t;s) = f_j(t) + \nabla f_j(s) (t - s) + 1/2 \|t - s\|^2_{C_j(s)}\]

where $C_j(⋅)$ is diagonal matrix of curvatures. (It suffices for each $f_j$ to have a Lipschitz smooth gradient.)

Each function $f_j : \mathcal{X}_j ↦ \mathbb{R}$ where conceptually $\mathcal{X}_j ⊆ \mathbb{R}^{M_j}$, but we allow more general domains.

There are two outer constructors (based on the positional arguments):

• The simple way (not type stable) provides
  • gradf vector of $J$ functions return gradients of $f_1,…,f_J$
  • curvf vector of $J$ functions z -> curv(z) that return a scalar or a vector of curvature values for each element of $z$
• The fancier way (type stable) provides
  • dot_gradf::AbstractVector{<:Function} = make_dot_gradf.(gradf) See make_dot_gradf.
  • dot_curvf::AbstractVector{<:Function} = make_dot_curvf.(curvf) See make_dot_curvf.

in

• u vector of $J$ arrays $u_1,…,u_J$ (typically vectors)
• v vector of $J$ arrays $v_1,…,v_J$ (typically vectors)

We require axes(u_j) == axes(v_j) for all $j=1,…,J$.

option

• α0::Real = 0f0 initial guess for step size
• ninner::Int = 5 # max number of inner iterations of MM line search
• work = LineSearchMMWork(u, v, α) pre-allocated work space for $u_j+α v_j$

LineSearchMMWork{Tz <: AbstractVector{<:AbstractArray}}

Workspace for storing $z_j = u_j + α v_j$ in MM-based line search.

If all of those $z_j$ arrays had the same eltype, then we could save memory by allocating just the longest vector needed. But for Unitful data they could have different eltypes and sizes, which would require a lot of reinterpret and reshape to handle. So we just allocate separate work arrays for each $j$.

A = Afft(samp::AbstractArray{Bool} ; T, dims, operator, work, ...)

Make a LinearMapAO operator object for (under-sampled) FFT, of type T, using given sampling pattern samp. Especially for compressed sensing MRI with Cartesian sampling.

Option:

• T::DataType = ComplexF32
• dims = 1:D apply fft/bfft only along these dimensions
• fft_forward::Bool = true Use false to have bfft! in forward model.
• operator::Bool = true set to false to return a LinearMapAM
• work::AbstractArray work space for in-place fft operations
• remaining arguments are passed to plan_fft

Output

Returns a LinearMapsAA.LinearMapA[M|O] object.

A = Afft(xdim::Dims, fdim ; T, operator, unitary, work, ...)

Make a <:LinearMapAX operator object for the FFT of an array of size xdim, along dimensions fdim, of type T.

In:

• xdim::Dims{D} array size
• fdim = 1:D apply fft/bfft only along these dimensions

Option:

• T::DataType = ComplexF32
• operator::Bool = true return a LinearMapAO set to false to return a LinearMapAM
• unitary::Bool = false set to true for unitary DFT
• work::AbstractArray work space for in-place fft operations
• remaining arguments are passed to plan_fft

Output

Returns a LinearMapsAA.LinearMapA[M|O] object.

Anufft(ω, N ; kwargs ...)

Make a LinearMapAO object of size length(ω) × prod(N). See nufft_init for options.

A, levels, mfun = Aodwt(dims ; level::Int=3, wt=wavelet(WT.haar))

Create orthogonal discrete wavelet transform (ODWT) LinearMapAA

in

• dims::Dims tuple of dimensions

option

• level::Int # of levels; default 3
• wt wavelet transform type (see Wavelets package); default Haar
• operator::Bool=true default to LinearMapAO
• T::DataType : Float32 by default; use ComplexF32 if needed

out

• A a LinearMapAX object
• scales array of size dims showing the scale of each coefficient

which is useful when imposing scale-dependent regularization

• mfun convenience function for A*X when X is a Matrix or Array (not vector)

2019-02-23 Jeff Fessler, University of Michigan

Asense(samp, smaps; T)

Construct a MRI encoding operator model for D-dimensional Cartesian sampling pattern samp and coil sensitivity maps smaps.

The input smaps can either be a D+1 dimensional array of size (size(samp)..., ncoil), or a Vector of ncoil arrays of size size(samp).

Input

• samp::AbstractArray{<:Bool} D-dimensional sampling pattern.
• smaps::Vector{<:AbstractArray{<:Number}} or AbstractArray{<:Number}

Option

• T::DataType = ComplexF32
• dims = 1:D apply fft/bfft only along these dimensions
• fft_forward::Bool = true Use false to have bfft! in forward model.
• unitary::Bool = false set to true for unitary DFT
• remaining arguments are passed to plan_fft

Output

Returns a LinearMapsAA.LinearMapAO object.

_show_struct(io::IO, ::MIME"text/plain", st::Any)

Informative way to show fields of a struct (composite type).

caller_name() or caller_name(;level=4)

Return "filename line fun():" as a string to describe where this function was called.

Stack levels:

• 1: #caller_name
• 2: caller_name()
• 3: function that invoked caller()
• 4: the calling function we want to return

Hence the default level is 4, but we increment it by one in case user says @show caller_name() in which case stack[3] is a macro expansion.

Z = diff_adj(dx, N::Dims{D} ; dims = 1:D)

Adjoint of finite differences of arrays along one or more dimensions. By default performs the same operations as $vec(Z) = [(I_{N_d} \otimes \cdots \otimes D_{N_1}); \dots; (D_{N_d} \otimes \cdots \otimes I_{N_1})]' * d$ where D_N denotes the N-1 × N 1D finite difference matrix and ⊗ denotes the Kronecker product, but does it efficiently without using spdiagm (or any SparseArrays function).

in

• dx vector of typical length N_d*...*(N_1-1) + ... + (N_d-1)*...*N_1
• N::Dims desired output size

option

• dims dimension(s) for performing adjoint finite differences; default 1:ndims(X)

out

• Z N_1 × ... × N_d array by default

d = diff_forw(X ; dims = 1:ndims(X))

Finite differences along one or more dimensions of an array, e.g., for anisotropic TV regularization.

By default performs the same operations as $d = [(I_{N_d} \otimes \cdots \otimes D_{N_1}); \dots; (D_{N_d} \otimes \cdots \otimes I_{N_1})] vec(X)$ where $D_N$ denotes the N-1 × N 1D finite difference matrix and ⊗ denotes the Kronecker product, but does it efficiently without using spdiagm (or any SparseArrays function).

Input dimension N must exceed 1 for each dimension specified by dims.

in

• X N_1 × ... × N_d array (typically an N-D image).

option

• dims dimension(s) for performing finite differences; default 1:ndims(X)

must have unique elements and be a nonempty subset of 1:ndims(X)

out

• d vector of default length N_d*...*(N_1-1) + ... + (N_d-1)*...*N_1

T = diff_map(N::Dims{D} ; dims = 1:D)

in

• N::Dims image size

out

• T LinearMapAA object for computing finite differences via T*x

diffl!(g::AbstractArray, x::AbstractArray, dims::AbstractVector{Int} ; ...)

When x is a N-dimensional array, the ith slice of the g array (along its last dimension) is the diffl! of x along dims[i]. This is useful for total variation (TV) and other regularizers that need finite differences along multiple dimensions.

diffl!(g::AbstractArray, x::AbstractArray, dim::Int ; ...)

Apply left finite difference operator to input array x, storing the result "in-place" in pre-allocated output array g. (The letter g is mnemonic for "gradient".)

Arrays g and x must have the same size, and cannot alias. By default, the "first" elements of g are zero for dimension dim. The default is dim=1.

Option:

• add::Bool = false use x[i] + x[i-1] instead of x[i] - x[i-1] (useful for certain diagonal majorizers).
• edge::Symbol = :zero set the first elements of dimension dim to 0
  • Choose edge=:circ to use circulant (aka periodic) boundary conditions.
  • Choose edge=:none to leave the first elements untouched.

Example

julia> x = [2, 6, 7]; g = similar(x); diffl!(g, x)
3-element Vector{Int64}:
 0
 4
 1

g = diffl(x::AbstractArray, dims::AbstractVector{Int} ; ...)

Allocating version of diffl! for dims

g = diffl(x::AbstractArray, dim::Int ; ...)

Allocating version of diffl! along dim

Example

julia> x = [1,2] .+ [10 30 70]; g = diffl(x, 2)
2×3 Matrix{Int64}:
 0  20  40
 0  20  40

g = diffl(x::AbstractArray ; ...)

Allocating version of diffl! along dim=1

diffl_adj!(z::AbstractArray, g::AbstractArray, dims::AbstractVector{Int} ; ...)

Adjoint of diffl! for multiple dimensions dims. Here g must have one more dimension than z.

diffl_adj!(z, g, dim::Int ; ...)

Adjoint of left finite difference diffl!, in-place. Arrays z and g must be same size. See diffl_adj for details.

z = diffl_adj(g::AbstractArray, dim::Int ; ...)

Allocating version of diffl_adj! along dim.

Example

julia> g = ones(Int,2,3); z = diffl_adj(g, 2)
2×3 Matrix{Int64}:
 -1  0  1
 -1  0  1

z = diffl_adj(g::AbstractArray ; ...)

Allocating version of diffl_adj! along dim=1.

Example

julia> g = [0, 2, 5]; z = diffl_adj(g)
3-element Vector{Int64}:
 -2
 -3
  5

z = diffl_adj(g::AbstractArray, dims::AbstractVector{Int} ; ...)

Allocating version of diffl_adj! for dims.

T = diffl_map(N::Dims{D}, dims::AbstractVector{Int} ; kwargs...)
T = diffl_map(N::Dims{D}, dim::Int ; kwargs...)

in

• N::Dims image size

options: see diffl!

• T::Type for LinearMapAA, default Float32
• operator::Bool = true use false for LinearMapAM

out

• T LinearMapAA object for computing finite differences via T*x

using diffl! and diffl_adj!

y = downsample1(x, down ; warn=true)

Downsample 1D vector by factor down.

in

• x [n1]
• down::Int downsampling factor

option

• warn::Bool warn if noninteger multiple; default isinteractive()

out

• y [n1/down]

y = downsample2(x, down ; warn=true, T)

Downsample by averaging by integer factors.

in

• x [nx ny]
• down can be a scalar (same factor for both dimensions) or a NTuple{2,Int}

option

• warn::Bool warn if noninteger multiple; default isinteractive()
• T::DataType specify output eltype; default eltype(x[1] / down[1])

out

• y [nx/down ny/down]

y = downsample3(x, down ; warn=true, T)

Downsample by averaging by integer factors.

in

• x (nx,ny,nz)
• down can be a scalar (same factor for all dimensions) or a NTuple{3,Int}

option

• warn::Bool warn if noninteger multiple; default true
• T::DataType specify output eltype; default eltype(x[1] / down[1])

out

• y (nx/down,ny/down,nz/down)

y = downsample_dim1(x, down ; warn::Bool)

Down-sample x by factor down along first dimension by averaging.

in

• x [n1 (Nd)]
• down::Int downsampling factor

option

• warn::Bool warn if non-integer multiple; default isinteractive()

out

• y [n1÷down (Nd)]

X = dtft(w, x ; n_shift=?)

multi-dimensional DTFT (DSFT)

$X[m] = \sum_{n=0}^{N-1} x[n] \exp(-i w[m,:] (n - n_{shift})), m=1,…,M$ where here n is a CartesianIndex

in

• w::AbstractMatrix{<:Real} (M,D) frequency locations ("units" radians/sample)
• x::AbstractArray{<:Number} [(Nd)] D-dimensional signal

option

• n_shift::AbstractVector{<:Real} often is N/2; default zeros(D)

out

• X::AbstractVector{ComplexF64} (M) DTFT

X = dtft(w, x ; n_shift=?)

1D DTFT

$X[m] = \sum_{n=0}^{N-1} x[n] \exp(-i w[m] (n - n_{shift})), m=1,…,M$

in

• w::AbstractVector{<:Real} (M) frequency locations ("units" radians/sample)
• x::AbstractVector{<:Number} (N) 1D signal

option

• n_shift::Real often is N/2; default 0

out

• X::AbstractVector{ComplexF64} [M] DTFT

x = dtft_adj(w, X, N ; n_shift=?)

adjoint for multi-dimensional DTFT (DSFT)

$x[n] = \sum_{m=1}^M X[m] \exp(i w[m,:] (n - n_{shift})), n=0,…,N-1$ where here n is a CartesianIndex

in

• X::AbstractVector{ComplexF64} (M) DTFT
• w::AbstractMatrix{<:Real} (M,D) frequency locations ("units" radians/sample)
• N::Dims (D) dimensions of signal x

option

• n_shift::AbstractVector{<:Real} often is N/2; default zeros(D)
• T::DataType default (eltype(w) == Float64) ? ComplexF64 : ComplexF32

out

• x::AbstractArray{<:Number} [(N)] D-dimensional signal

x = dtft_adj(w, X, N ; n_shift=?)

adjoint for 1D DTFT

$x[n] = \sum_{m=1}^M X[m] \exp(i w[m] (n - n_{shift})), n=0,…,N-1$

This is the adjoint (transpose) of dtft, not an inverse DTFT.

in

• w::AbstractVector{<:Real} (M) frequency locations ("units" radians/sample)
• X::AbstractVector{ComplexF64} (M) spectrum values
• N::Int size of signal x

option

• n_shift::Real often is N/2; default 0
• T::DataType output data type; default ComplexF64

out

• x::AbstractVector{<:Number} (N) signal

d = dtft_init(w, N ; n_shift=?)

for multi-dimensional DTFT (DSFT)

in

• w::AbstractMatrix{<:Real} (M,D) frequency locations ("units" radians/sample)
• N::Dims [D] dimensions of signal x

option

• n_shift::AbstractVector{<:Real} often is N/2; default zeros(D)
• T::DataType default ComplexF64 for testing NUFFT accuracy

out

• d::NamedTuple with fields dtft = x -> dtft(x), adjoint = y -> dtft_adj(y), A=LinearMapAO

d = dtft_init(w, N ; n_shift=?)

For 1D DTFT

in

• w::AbstractVector{<:Real} (M) frequency locations ("units" radians/sample)
• N::Int size of signal x

option

• n_shift::Real often is N/2; default 0
• T::DataType default ComplexF64 for testing NUFFT accuracy

out

• d::NamedTuple with fields dtft = x -> dtft(x), adjoint = y -> dtft_adj(y), A=LinearMapAO

embed!(array, v, mask ; filler=0)

embed vector v of length sum(mask) into elements of array where mask is true, setting the remaining elements to filler (default 0).

array = embed(matrix::AbstractMatrix{<:Number}, mask::AbstractArray{Bool})

Embed each column of matrix into mask then cat along next dimension In:

• matrix [sum(mask) L]
• mask [(N)]

Out:

• array [(N) L]

array = embed(v, mask ; filler=0)

embed vector v of length sum(mask) into elements of an array where mask is true; see embed!.

x = eql_root(a,b,c)

Numerically stable method for computing the positive root of the quadratic polynomial -ax^2 - 2bx + c, a >= 0. Assumes solvable equations; will throw otherwise.

in

• a The negative of the x^2 term. Must be positive.
• b Half the negative of the x term.
• c The constant term.

out

• x The positive root that satisfies 0 = -ax^2 - 2bx + c.

D = exp_mult(A, u, v ; warnboth)

Memory efficient and fast implementation of D = A' * exp(-u * v^T) that is useful for B0-field-corrected MRI image reconstruction.

in

• A [N L] matrix
• u [N] vector
• v [M] vector
• warnboth warn if both u and v are complex; default: true

out

• D [L M] complex vector: D = A' * exp(-u * v^T)

D_lm = sum_n A_nl^* exp(-u_n v_m)

exp_xform(x, u, v ; mode::Symbol = :matrix)

in

• x [N L] possibly complex vector(s)
• u [D N] possibly complex vectors
• v [D M] possibly complex vectors
• mode::Symbol :matrix (default) | :element | :row | :column

out

• y [M L] typically complex vector

y[m,l] = sum_n x[n,l] exp(-sum(u[:,n] .* v[:,m]))

Iterates through subsets of the ML matrix designated by :mode (i.e. row, column, element, or just computing the entire matrix) This is the 'slow' 'exact' transform model for MRI.

Output type will depend on input types.

genkspace

Generate the proper length of k-space trajectory.

It linearly interpolates the output of genspiral to the correct length() & takes care of the rotations for the interleaves.

• ld is the length of the data
• nint is the number of interleaves

Brad Sutton; University of Michigan

Gx, Gy, kx, ky, sx, sy, gts = genspi(...)

This is translation of C code from scanner: exactly what is played out to gradients at 4us.

Multi-shot spiral design uses Duyn's approximate slewrate limited design augmented with archimedian gmax limit

in [args]

• D = FOV; cm

• N = matrix size()

• Tmax = longest acquisition allowed; s

• dts = output sample spacing; s

• gtype = trajectory type()

option [CVs]

• nl = number of interleaves
• gamp = design grad max; G/cm
• gslew = design slew rate; mT/m/ms
• nramp = number of rampdown points; default 0

out

• Gx; Gy

time is in sec()

• rev 0 12/26/98 original
• rev 1 4/15/99 little better calc of ts

Borrowed from Doug Noll; Univ. of Michigan. Modified to take more input cv's.

getindex!(y::AbstractVector, x::AbstractArray{T,D}, mask::AbstractArray{Bool,D})

Equivalent to the in-place y .= x[mask] but is non-allocating.

For non-Boolean indexing, just use @views y .= A[index], per https://discourse.julialang.org/t/efficient-non-allocating-in-place-getindex-for-bitarray/42268

yi = interp1(x, y, xi ; how=Gridded(Linear()), extrap=0)

1D interpolation of y = f(x) at points xi

In:

• x::AbstractVector{<:Real}
• y::AbstractVector{<:Number}

Option:

• how::Interpolations.InterpolationType default Gridded(Linear())
• extrap::Any how to extrapolate, e.g., Flat(); default 0

other options from Interpolations.jl are Line() Periodic() Reflect() Throw()

Output is same size as input xi

ir_dump(x::Any ; io::IO = stdout)
ir_dump(io::IO, x::Any)

Show all the fields of a structure or NamedTuple more nicely than dump() does

data = ir_load_brainweb_t1_256()

Load brainweb T1-weighted MRI slice of size 256 × 256.

(odata, σ, Vr) = ir_mri_coil_compress(idata ; ncoil)

MRI coil compression via PCA. Given multiple MRI surface coil images (idata), use SVD/PCA to find a smaller number of virtual coil images (odata).

In:

• idata [(N) n_in]: noisy complex images (2D or 3D) for each coil

Option:

• ncoil Desired # of virtual coils (default: 1)

Out:

• odata [(N) ncoil]: virtual coil images
• σ [n_in]: singular values.
• Vr [n_in, ncoil]: compression matrix for reducing other data.

todo: currently ignores noise correlations

kspace = ir_mri_kspace_ga_radial(; Nro=?, Nspoke=?, ...)

Generate k-space sampling pattern for "golden angle" radial sampling.

option

• Nro:Int number of samples in each readout/spoke, default 256
• Nspoke::Int number of spokes, default 1
• start::Real first angle in series [radians], default π/2
• angle::Real angular spacing [radians], default GA
• delta_ro::Real readout spacing, default 1/Nro
• shift::Real shift due to gradient delays, default 0
  • radial sample locations are ir * delta_ro
  • where ir = [-(Nro/2 - 1):1:Nro/2] + shift

out

• kspace [Nro Nspoke 2] (Float32)

kx and ky k-space locations for Nspoke*Nro samples in interval (-0.5 0.5] for default shift, delta_ro so default units are "cycles / sample"

(smap,info) = ir_mri_sensemap_sim( :all ; kwargs)

Like ir_mri_sensemap_sim but also returns info with data for all coils, mainly for testing and plotting.

(smap,info) = ir_mri_sensemap_sim( ir_ic_pair ; kwargs)

Like ir_mri_sensemap_sim but also returns info with data for specific coils where ir_ic_pair::Vector{Tuple{Int,Int}}. (Usually used internally only.)

• info::NamedTuple geometry information for plots

smap = ir_mri_sensemap_sim(...)

Simulate 2D or 3D sensitivity maps for sensitivity-encoded MRI based on grivich:00:tmf.

This code makes maps for multiple coils, but does not model coupling between coils, so most likely it is an approximation at best.

option

• dims::Dims image size; default (64, 64)
• dx::Real pixel/voxel dimension; default: 3
• dy::Real pixel/voxel dimension; default: dx
• dz::Real ""
• ncoil::Int # of coils total; default 4
• nring::Int # of rings of coils; default 1
• rcoil::Real coil radius; default dx * nx / 2 * 0.50
• dz_coil ring spacing in z; default nz*dz/nring
  • (3D geometry is a cylinder)
• coil_distance::Real distance of coil center from isocenter
  • for central ring of coils as a multiple of FOVx,
  • where FOVx=nx*dx; default 1.2
• orbit::Real default 360 [degrees]
• orbit_start::AbstractVector{<:Real} = fill(0, nring) [degrees]
• scale::Symbol
  • :none (default)
  • ssos_center make SSoS of center = 1

out

• smap [dims ncoil] simulated sensitivity maps (complex!)

All length parameters must have same units (e.g., mm or cm)

(smap, info) = ir_mri_sensemap_sim_do()

ir_mri_smap1()

Based on grivich:00:tmf for a circular coil in "x-y plane" of radius "a"

Note that coil x-y plane is not same as object x-y plane!

Returns (i,j,k) components of $B$ vector for each (x,y,z) location.

ir_mri_smap_r(r, z)

Function for testing near 0.

jinc(x)

Return jinc(x) = J1(pi*x)/(2x), where J1 is a Bessel function of the first kind.

Units of x are typically cycles/m.

Return type is promote_type(typeof(x), Float32).

α = line_search_mm(args...; opt, fun, kwargs...)

Line-search based on majorize-minimize (MM) approach. This is a wrapper around the iterator LineSearchMM. See its constructors for args and other kwargs.

option

• fun(state) User-defined function to be evaluated with the state initially and then after each iteration.
• out::Union{Nothing,Vector{Any}} = nothing optional place to store result of fun for iterates 0,…,ninner: (All missing by default.) This is aVector{Any}of lengthninner+1`.

output

• α final iterate

This function mutates the optional arguments out and work.

make_dot_curvf(curv::Function, [x, Tf::DataType; w = similar(x)])

Make a function with arguments (v, x) that computes the dot product between abs2.(v) and curv(x) where curv(x) is a curvature function associated with a quadratic majorizer for some real-valued cost function f(x), evaluated at array x, typically for use in a line-search method.

• For a single-argument gradient function curv(x), this returns the version equivalent to (v, x) -> dot(abs2.(v), curv(x)). This version will be allocating, unless curv has its own internal workspace due to a closure. This version expects only the argument curv.

• For a two-argument mutating gradient function curv!(w, x), this returns a function equivalent to (v, x) -> dot(abs2.(v), curv!(w, x)) using the keyword argument w as the work array.

• If curv is simply a real number (a Lipschitz constant), possibly with units, then this returns the function (v, x) -> curv * sum(abs2, v).

If f(x) maps an Array x of elements with units Ux into real numbers with units Uf, then its curvature has units Uf/Ux^2. Those units are relevant to defining the work array w.

in

• curv::Function see above
• x an array whose size and eltype is used to allocate w

option

• Tf::DataType = typeof(one(eltype(x))) Specify eltype of function f(x), defaulting to unitless.
• w = similar(x, typeof(oneunit(Tf) / oneunit(eltype(x))^2)) work space for gradient calculation, with appropriate units (if needed).

out

• (v, x) -> dot(abs2.(v), curv([w,] x)) or equivalent.

make_dot_gradf(grad::Function, [x, Tf::DataType; w = similar(x)])

Make a function with arguments (v, x) that computes the dot product between an array v and the gradient grad(x) of some real-valued cost function f(x), evaluated at array x, typically for use in a line-search method.

• For a single-argument gradient function grad(x), this returns the version (v, x) -> dot(v, grad(x)). This version will be allocating, unless grad has its own internal workspace due to a closure. This version expects only the argument grad.

• For a two-argument mutating gradient function grad!(w, x), this returns a function (v, x) -> dot(v, grad!(w, x)) using the keyword argument w as the work array.

If f(x) maps an Array x of elements with units Ux into real numbers with units Uf, then the gradient ∇f has units Uf/Ux. Those units are relevant to defining the work array w.

in

• grad::Function see above
• x an array whose size and eltype is used to allocate w
• Tf::DataType = typeof(one(eltype(x))) Specify eltype of function f(x), defaulting to unitless.

option

• w = similar(x, typeof(oneunit(Tf) / oneunit(eltype(x)))) work space for gradient calculation, with appropriate units (if needed).

out

• (v, x) -> dot(v, grad([w,] x)) or equivalent.

Example. Consider the cost function f(x) = sum(sin.(x)) that has gradient g(x) = cos.(x). The simple way here is make_dot_gradf(g), which will return the function (v, x) = dot(v, g(x)) === dot(v, cos.(x)). That will work fine, but the cos.(x) step will be allocating. A better way (for repeated use) is

function g!(work, x)
   @. work = cos(x) # mutating
   return work
end

Then make_dot_gradf(g!, x) will return a function (v,x) that does not allocate (except when first generated via a closure). For this particular example, an even better approach is to directly define dot_gradf(v,x) = sum(vx -> vx[1]*cos(vx[2]), zip(v,x)), but make_dot_gradf is here for more general cases.

y = map_many(fun::Function, x::AbstractArray{<:Any}, idim::Dims)

Apply a function fun to leading slices of input x; cousin of mapslices

in

• fun::Function maps input of size idim to output of some size odim
• x [idim ldim]

out

• y [odim ldim]

Example: if fun maps array of size (1,2) to array of size (3,4,5) and if input x has size (1,2,7,8) then output y will have size (3,4,5,7,8) where y[:,:,:,i,j] = fun(x[:,:,i,j])

mask_or(mask)

compress 3D mask to 2D by logical or along z direction

mask_outline(mask)

return outer boundary of 2D mask (or mask_or for 3D)

maskit(x::AbstractArray{<:Number})

opposite of embed

max_percent_diff(s1, s2, [options])

Compute the "maximum percent difference" between two signals: s1, s2.

Default is to normalize by maximum(abs, s1).

options

• maxboth::Bool = false use max of both arguments to normalize
• normalize::Bool = false normalize each before comparing

kspace, omega, gradxy = mri_kspace_spiral( [options] )

Make k-space spiral trajectory based on GE 3T scanner constraints

Option:

• N dimention of reconstructed image
• Nt # of time points
• fov field of view in cm
• dt time sampling interval out; default 5e-6 sec
• gamp::Real design gradient amplitude max, G/cm; default 2.2
• gslew::Int design slew rate, mT/m/ms; default 180

Out:

• kspace [Nt,2] kspace trajectory [kx ky] in cycles/cm, NO: cycles/FOV
• omega [Nt,2] "" in radians
• gradxy [Nt 2] gradient waveforms in (units?)

kspace, omega, wi = mri_trajectory( ; ktype, N, fov, arg_wi, kwargs...)

Generate kspace trajectory samples and density compensation functions.

option

• ktype::Symbol k-space trajectory type; default :radial
• N::Dims target image size; default (32,30)
• fov field of view in x and y (and z); default (250,250) mm
• arg_wi options to pass to ir_mri_density_comp - not yet done
• kwargs options for the specific trajectory

out

• kspace [Nk 2|3] kspace samples in units 1/fov
• omega [Nk 2|3] trajectory samples over [-π,π)
• wi [Nk 1] (optional) density compensation factors

trajectory types:

• :cartesian
• :radial
• :cart_y_2
• :random
• :half8
• :epi_sin
• :spiral0 :spiral1 :spiral3
• :rosette3
• :epi_under
• :gads (emulate golden-angle data sharing per winkelmann:07:aor)

omega, wi = mri_trajectory_gads(N, fov ; ...)

Emulate 2D golden angle radial sampling with data sharing

option

• Nro # of samples in each readout/spoke
• shift shift along read-out due to gradient delays (stress)
• kmax_frac fractions of maximum krad (0.5) for rings (annuli)
• under under-sampling factor for each annulus

mri_trajectory_radial()

option

• na_nr default ensures proper sampling at edge of k-space
• na angular spokes; default: na_nr * nr
• nr radial samples per spoke
• ir default: 0:nr

todo: generalize to 3D using barger:02:trc

mri_trajectory_rosette3(N, fov ; ...)

3d rosette, with default parameters from bucholz:08:miw

option

• omax maximum omega
• nt time samples (65.536 ms for 4 usec dt)
• dt time sample spacing (4 usec)
• ti time samples

(x,out) = ncg(B, gradf, curvf, x0 ; ...)

Nonlinear preconditioned conjugate gradient algorithm to minimize a general "inverse problem" cost function of the form $\Psi(x) = \sum_{j=1}^J f_j(B_j x)$ where each function $f_j(v)$ has a quadratic majorizer of the form

\[q_j(v;u) = f_j(u) + \nabla f_j(u) (v - u) + 1/2 \|v - u\|^2_{C_j(u)}\]

where $C_j(u)$ is diagonal matrix of curvatures. (It suffices for each $f_j$ to have a Lipschitz smooth gradient.)

This CG method uses a majorize-minimize (MM) line search.

in

• B vector of $J$ blocks $B_1,…,B_J$
• gradf vector of $J$ functions return gradients of $f_1,…,f_J$
• curvf vector of $J$ functions z -> curv(z) that return a scalar or a vector of curvature values for each element of $z$
• x0 initial guess; need length(x) == size(B[j],2) for $j=1,…,J$

Usually x0 is a Vector but it can be an Array if each B_j is a linear operator (e.g., LinearMapAO) of suitable "dimensions".

option

• niter # number of outer iterations; default 50
• ninner # number of inner iterations of MM line search; default 5
• P # preconditioner; default I
• betahow "beta" method for the search direction; default :dai_yuan
• fun User-defined function to be evaluated with two arguments (x,iter).
  • It is evaluated at (x0,0) and then after each iteration.

output

• x final iterate
• out [niter+1] (fun(x0,0), fun(x1,1), ..., fun(x_niter,niter))
  • (all 0 by default). This is an array of length niter+1

(x,out) = ncg(grad, curv, x0, ...)

Special case of ncg (nonlinear CG) for minimizing a cost function whose gradient is grad(x) and that has a quadratic majorizer with diagonal Hessian given by curv(x). Typically curv = (x) -> L where L is the Lipschitz constant of grad.

nufft_eltype(::DataType)

ensure plan_nfft eltype is Float32 or Float64

w, errs = nufft_errors( ; M=401, w=?, N=513, n_shift=0, ...)

Compute NUFFT approximation errors (for signal of length N of unit norm), for given digital frequency values w, i.e., Ω. Default w is range(0, 2π/N, M).

p = nufft_init(w, N ; nfft_m=4, nfft_sigma=2.0, pi_error=true, n_shift=0)

Setup 1D NUFFT, for computing fast $O(N \log N)$ approximation to

$X[m] = \sum_{n=0}^{N-1} x[n] \exp(-i w[m] (n - n_{shift})), m=1,…,M$

in

• w::AbstractArray{<:Real} [M] frequency locations (aka Ω, units radians/sample)
  • eltype(w) determines the plan_nfft type; so to save memory use Float32!
  • size(w) determines odim for A if operator=true
• N::Int signal length

option

• nfft_m::Int see NFFT.jl documentation; default 4
• nfft_sigma::Real "", default 2.0
• n_shift::Real often is N/2; default 0
• pi_error::Bool throw error if $|w| > π$, default true
  • Set to false only if you are very sure of what you are doing!
• do_many::Bool support extended inputs via map_many? default true
• operator::Bool=true set to false to make A an LinearMapAM

out

• p NamedTuple

(nufft = x -> nufft(x), adjoint = y -> nufft_adj(y), A::LinearMapAO)

The default settings are such that for a 1D signal of length N=512, the worst-case error is below 1e-5 which is probably adequate for typical medical imaging applications. To verify this statement, run nufft_plot1() and see plot.

p = nufft_init(w, N ; nfft_m=4, nfft_sigma=2.0, pi_error=true, n_shift=?)

Setup multi-dimensional NUFFT, for computing fast $O(N \log N)$ approximation to

$X[m] = \sum_{n=0}^{N-1} x[n] \exp(-i w[m,:] (n - n_{shift})), m=1,…,M$

in

• w::AbstractArray{<:Real} [M,D] frequency locations (aka Ω, units radians/sample)
  • eltype(w) determines the plan_nfft type; so to save memory use Float32!
  • size(w)[1:(end-1)] determines odim if operator=true
• N::Dims{D} signal dimensions

option

• nfft_m::Int see NFFT.jl documentation; default 4
• nfft_sigma::Real "", default 2.0
• n_shift::AbstractVector{<:Real} (D) often is N/2; default zeros(D)
• pi_error::Bool throw error if $|w| > π$, default true
  • Set to false only if you are very sure of what you are doing!
• do_many::Bool support extended inputs via map_many? default true
• operator::Bool=true set to false to make A a LinearMapAM

The default do_many option is designed for parallel MRI where the k-space sampling pattern applies to every coil. It may also be useful for dynamic MRI with repeated sampling patterns. The coil and/or time dimensions must come after the spatial dimensions.

out

• p NamedTuple with fields nufft = x -> nufft(x), adjoint = y -> nufft_adj(y), A=LinearMapAO (Using operator=true allows the LinearMapAO to support do_many.)

nufft_typer(T::DataType, x::AbstractArray{<:Real} ; warn::Bool=true)

type conversion wrapper for nfft()

(x,out) = ogm_ls(B, gradf, curvf, x0; niter=?, ninner=?, fun=?)

OGM with a line search Drori&Taylor to minimize a general "inverse problem" cost function of the form $\Psi(x) = \sum_{j=1}^J f_j(B_j x)$ where each function $f_j(v)$ has a quadratic majorizer of the form

\[q_j(v;u) = f_j(u) + \nabla f_j(u) (v - u) + 1/2 \|v - u\|^2_{C_j(u)}\]

where $C_j(u)$ is diagonal matrix of curvatures. (It suffices for each $f_j$ to have a Lipschitz smooth gradient.)

This OGM method uses a majorize-minimize (MM) line search.

in

• B vector of $J$ blocks $B_1,…,B_J$
• gradf vector of $J$ functions return gradients of $f_1,…,f_J$
• curvf vector of $J$ functions z -> curv(z) that return a scalar or a vector of curvature values for each element of $z$
• x0 initial guess; need length(x) == size(B[j],2) for $j=1,…,J$

option

• niter # number of outer iterations; default 50
• ninner # number of inner iterations of MM line search; default 5
• fun User-defined function to be evaluated with two arguments (x,iter).
• It is evaluated at (x0,0) and then after each iteration.

output

• x final iterate
• out (niter+1) (fun(x0,0), fun(x1,1), ..., fun(x_niter,niter))
  • (all 0 by default). This is a vector of length niter+1.

(x,out) = ogm_ls(grad, curv, x0, ...)

Special case of ogm_ls (OGM with line search) for minimizing a cost function whose gradient is grad(x) and that has a quadratic majorizer with diagonal Hessian given by curv(x). Typically curv = (x) -> L where L is the Lipschitz constant of grad.

x, out = pogm_restart(x0, Fcost, f_grad, f_L ;
f_mu=0, mom=:pogm, restart=:gr, restart_cutoff=0.,
bsig=1, niter=10, g_prox=(z,c)->z, fun=...)

Iterative proximal algorithms (PGM=ISTA, FPGM=FISTA, POGM) with restart.

in

• x0 initial guess
• Fcost function for computing the cost function value $F(x)$
  • (needed only if restart === :fr)
• f_grad function for computing the gradient of $f(x)$
• f_L Lipschitz constant of the gradient of $f(x)$

option

• f_mu strong convexity parameter of $f(x)$; default 0.
  • if f_mu > 0, $(\alpha, \beta_k, \gamma_k)$ is chosen by Table 1 in [KF18]
• g_prox function g_prox(z,c) for the proximal operator for $g(x)$
  • g_prox(z,c) computes $argmin_x 1/2 \|z-x\|^2 + c \, g(x)$
• mom momentum option
  • :pogm POGM (fastest); default!
  • :fpgm (FISTA), $\gamma_k = 0$
  • :pgm PGM (ISTA), $\beta_k = \gamma_k = 0$
• restart restart option
  • :gr gradient restart; default!
  • :fr function restart
  • :none no restart
• restart_cutoff for :gr restart if cos(angle) < this; default 0.
• bsig gradient "gamma" decrease option (value within [0 1]); default 1
  • see $\bar{\sigma}$ in [KF18]
• niter number of iterations; default 10
• fun function(iter, xk, yk, is_restart) user-defined function evaluated each iter with secondary xk, primary yk, and boolean is_restart indicating whether this iteration was a restart

out

• x final iterate
  • for PGM (ISTA): $x_N = y_N$
  • for FPGM (FISTA): primary iterate $y_N$
  • for POGM: secondary iterate $x_N$, see [KF18]
• out [fun(0, x0, x0, false), fun(1, x1, y1, is_restart), ...] array of length [niter+1]

Optimization Problem: Nonsmooth Composite Convex Minimization

• $argmin_x F(x), F(x) := f(x) + g(x))$
  • $f(x)$ smooth convex function
  • $g(x)$ convex function, possibly nonsmooth and "proximal-friendly" [CP11]

Optimization Algorithms:

Accelerated First-order Algorithms when $g(x) = 0$ [KF18] iterate as below for given coefficients $(\alpha, \beta_k, \gamma_k)$

• For k = 0,1,...
  • $y_{k+1} = x_k - \alpha f'(x_k)$ : gradient update
  • $x_{k+1} = y_{k+1} + \beta_k (y_{k+1} - y_k) + \gamma_k (y_{k+1} - x_k)$ : momentum update

Proximal versions of the above for $g(x) \neq 0$ are in the below references, and use the proximal operater $prox_g(z) = argmin_x {1/2\|z-x\|^2 + g(x)}$.

• Proximal Gradient method (PGM or ISTA) - $\beta_k = \gamma_k = 0$. [BT09]
• Fast Proximal Gradient Method (FPGM or FISTA) - $\gamma_k = 0$. [BT09]
• Proximal Optimized Gradient Method (POGM) - [THG15]
• FPGM(FISTA) with Restart - [OC15]
• POGM with Restart - [KF18]

references

• [CP11] P. L. Combettes, J. C. Pesquet,

"Proximal splitting methods in signal processing," Fixed-Point Algorithms for Inverse Problems in Science and Engineering, Springer, Optimization and Its Applications, 2011.

• [KF18] D. Kim, J.A. Fessler,

"Adaptive restart of the optimized gradient method for convex optimization," 2018 Arxiv:1703.04641, [http://doi.org/10.1007/s10957-018-1287-4]

• [BT09] A. Beck, M. Teboulle:

"A fast iterative shrinkage-thresholding algorithm for linear inverse problems," SIAM J. Imaging Sci., 2009.

• [THG15] A.B. Taylor, J.M. Hendrickx, F. Glineur,

"Exact worst-case performance of first-order algorithms for composite convex optimization," Arxiv:1512.07516, 2015, SIAM J. Opt. 2017 [http://doi.org/10.1137/16m108104x]

Copyright 2017-3-31, Donghwan Kim and Jeff Fessler, University of Michigan 2018-08-13 Julia 0.7.0 2019-02-24 interface redesign

v1,σ1 = poweriter(A; niter=?, ...)

Determine first right singular vector v1 and first singular value σ1 of A by applying power iteration to A'A

in

• A M × N matrix

option

• niter default 200
• x0 initial guess of v1
• tol stopping tolerance for s1, default 1e-6
• chat::Bool verbose? default false

out

• v1 [N] principal right singular vector
• σ1 spectral norm of A

rect(x::Real)

Unit width rect function. Potential problem? Bring up with fess.

y = reverser(x, dims)

reverse array along specified dimensions (or all if unspecified)

rmsd = rmsd100(x, y ; mask)

Compute 100 * RMSD (root mean squared difference) between x and y within domain mask.

in

• x : array
• y : another array of same size

option

• mask::Array{Bool} : domain over which to compute the RMSE; default trues(size(x))

out

• rmsd : rmsd of x vs y within mask in %

(xr,yr) = rotate2d(x, y, theta)

2D rotation

snr2sigma(db, yb)

Convert SNR in dB to noise σ for complex gaussian noise. No sqrt(2) factors is needed here because randn(Complex{Float}) already accounts for that. (See randn documentation.)

@shows expr

Show the type and size of an expression expr (typically a variable).

More concise output than @show and typically this is all that is needed when debugging code.