diff --git a/dev/examples/1-mirt/index.html b/dev/examples/1-mirt/index.html index a603f184..cf5bf704 100644 --- a/dev/examples/1-mirt/index.html +++ b/dev/examples/1-mirt/index.html @@ -1,2 +1,2 @@ -MIRT overview · MIRT.jl
Edit on GitHub
+MIRT overview · MIRT.jl
Edit on GitHub
diff --git a/dev/index.html b/dev/index.html index 2aa6473a..0995cf7c 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -Home · MIRT.jl
Edit on GitHub

MIRT: The Michigan Image Reconstruction Toolbox

MIRT.jl is a collection of Julia functions for performing image reconstruction and solving related inverse problems.

It is very much still under construction, although there are already enough tools to solve useful problems like compressed sensing MRI reconstruction. Trying the demos is a good way to get started.

The documentation is even more still under construction.

Further examples will likely be shared at JuliaImageRecon; see Examples there.

+Home · MIRT.jl
Edit on GitHub

MIRT: The Michigan Image Reconstruction Toolbox

MIRT.jl is a collection of Julia functions for performing image reconstruction and solving related inverse problems.

It is very much still under construction, although there are already enough tools to solve useful problems like compressed sensing MRI reconstruction. Trying the demos is a good way to get started.

The documentation is even more still under construction.

Further examples will likely be shared at JuliaImageRecon; see Examples there.

diff --git a/dev/reference/index.html b/dev/reference/index.html index 40f827ed..6779d671 100644 --- a/dev/reference/index.html +++ b/dev/reference/index.html @@ -1,25 +1,25 @@ -Function References · MIRT.jl
Edit on GitHub

Function References

MIRT.Afft — Method
A = Afft(samp::AbstractArray{Bool} ; T::DataType = ComplexF32)

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

Option:

  • operator::Bool = true set to false to return a LinearMapAM
  • work::AbstractArray work space for in-place fft operations
source
MIRT.Anufft — Method
Anufft(ω, N ; kwargs ...)

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

source
MIRT.Aodwt — Method
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

source
MIRT.caller_name — Method

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.

source
MIRT.diff_adj — Method
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
source
MIRT.diff_forw — Method
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
source
MIRT.diff_map — Method
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
source
MIRT.diffl! — Method
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.

source
MIRT.diffl! — Method
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)
+Function References · MIRT.jl
Edit on GitHub
Function References
MIRT.Afft — Method
A = Afft(samp::AbstractArray{Bool} ; T::DataType = ComplexF32)
Make a LinearMapAO object for (under-sampled) FFT, of type T, using given sampling pattern samp. Especially for compressed sensing MRI with Cartesian sampling.
Option:
  • operator::Bool = true set to false to return a LinearMapAM
  • work::AbstractArray work space for in-place fft operations
source
MIRT.Anufft — Method
Anufft(ω, N ; kwargs ...)
Make a LinearMapAO object of size length(ω) × prod(N). See nufft_init for options.
source
MIRT.Aodwt — Method
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
source
MIRT.caller_name — Method
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.
source
MIRT.diff_adj — Method
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
source
MIRT.diff_forw — Method
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
source
MIRT.diff_map — Method
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
source
MIRT.diffl! — Method
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.
source
MIRT.diffl! — Method
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
source
MIRT.diffl — Method
g = diffl(x::AbstractArray, dims::AbstractVector{Int} ; ...)
Allocating version of diffl! for dims
source
MIRT.diffl — Method
g = diffl(x::AbstractArray, dim::Int ; ...)
Allocating version of diffl! along dim
Example
julia> x = [1,2] .+ [10 30 70]; g = diffl(x, 2)
+ 1
source
MIRT.diffl — Method
g = diffl(x::AbstractArray, dims::AbstractVector{Int} ; ...)
Allocating version of diffl! for dims
source
MIRT.diffl — Method
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
source
MIRT.diffl — Method
g = diffl(x::AbstractArray ; ...)
Allocating version of diffl! along dim=1
source
MIRT.diffl_adj! — Method
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.
source
MIRT.diffl_adj! — Method
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.
source
MIRT.diffl_adj — Method
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)
+ 0  20  40
source
MIRT.diffl — Method
g = diffl(x::AbstractArray ; ...)
Allocating version of diffl! along dim=1
source
MIRT.diffl_adj! — Method
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.
source
MIRT.diffl_adj! — Method
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.
source
MIRT.diffl_adj — Method
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
source
MIRT.diffl_adj — Method
z = diffl_adj(g::AbstractArray ; ...)
Allocating version of diffl_adj! along dim=1.
Example
julia> g = [0, 2, 5]; z = diffl_adj(g)
+ -1  0  1
source
MIRT.diffl_adj — Method
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
source
MIRT.diffl_adj — Method
z = diffl_adj(g::AbstractArray, dims::AbstractVector{Int} ; ...)
Allocating version of diffl_adj! for dims.
source
MIRT.diffl_map — Method
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!
source
MIRT.downsample1 — Method
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]
source
MIRT.downsample2 — Method
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]
source
MIRT.downsample3 — Method
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]
source
MIRT.downsample_dim1 — Method
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)]
source
MIRT.dtft — Method
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
source
MIRT.dtft — Method
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
source
MIRT.dtft_adj — Method
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
source
MIRT.dtft_adj — Method
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}	signal [N]
source
MIRT.dtft_init — Method
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`
source
MIRT.dtft_init — Method
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`
source
MIRT.embed! — Method
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).
source
MIRT.embed — Method
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]
source
MIRT.embed — Method
array = embed(v, mask ; filler=0)
embed vector v of length sum(mask) into elements of an array where mask is true; see embed!.
source
MIRT.diffl_adj — Method
z = diffl_adj(g::AbstractArray, dims::AbstractVector{Int} ; ...)
Allocating version of diffl_adj! for dims.
source
MIRT.diffl_map — Method
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!
source
MIRT.downsample1 — Method
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]
source
MIRT.downsample2 — Method
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]
source
MIRT.downsample3 — Method
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]
source
MIRT.downsample_dim1 — Method
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)]
source
MIRT.dtft — Method
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
source
MIRT.dtft — Method
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
source
MIRT.dtft_adj — Method
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
source
MIRT.dtft_adj — Method
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}	signal [N]
source
MIRT.dtft_init — Method
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`
source
MIRT.dtft_init — Method
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`
source
MIRT.embed! — Method
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).
source
MIRT.embed — Method
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]
source
MIRT.embed — Method
array = embed(v, mask ; filler=0)
embed vector v of length sum(mask) into elements of an array where mask is true; see embed!.
source
MIRT.eql_root — Method
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.
source
MIRT.exp_mult — Method
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)
source
MIRT.exp_xform — Method
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.
source
MIRT.genkspace — Method
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
source
MIRT.genspi — Method
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.
source
MIRT.getindex! — Method
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
source
MIRT.interp1 — Method
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
source
MIRT.ir_dump — Method
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
source
MIRT.ir_mri_coil_compress — Method
(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
source
MIRT.ir_mri_kspace_ga_radial — Method
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"
2015-07 Mai Le, original Matlab version
2015-07-04 Jeff Fessler, minor changes to Matlab version
source
MIRT.ir_mri_sensemap_sim — Method
(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.
source
MIRT.ir_mri_sensemap_sim — Method
(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
source
MIRT.ir_mri_sensemap_sim — Method
smap = ir_mri_sensemap_sim(...)
Simulate 2D or 3D sensitivity maps for sensitivity-encoded MRI based on grivich:00:tmf http://doi.org/10.1119/1.19461
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)
source
MIRT.ir_mri_smap1 — Method
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
source
MIRT.jinc — Method
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).
source
MIRT.map_many — Method
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])
source
MIRT.mask_or — Method
mask_or(mask)
compress 3D mask to 2D by logical or along z direction
source
MIRT.max_percent_diff — Method
d = 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 use max of both arguments to normalize; default false
  • normalize::Bool normalize each before comparing; default false
source
MIRT.mri_kspace_spiral — Method
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?)
source
MIRT.mri_trajectory — Method
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)
source
MIRT.mri_trajectory_gads — Method
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
source
MIRT.mri_trajectory_radial — Method
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
source
MIRT.mri_trajectory_rosette3 — Method
mri_trajectory_rosette3(N, fov ; ...)
3d rosette, with default parameters from bucholz:08:miw
options: omax: maximum omega nt : time samples (65.536 ms for 4 usec dt) dt : time sample spacing (4 usec) ti : time samples
source
MIRT.ncg — Method
(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		array of $J$ blocks $B_1,...,B_J$
  • gradf	array of $J$ functions return gradients of $f_1,...,f_J$
  • curvf	array 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
source
MIRT.ncg — Method
(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
source
MIRT.ndgrid — Method
(xx,yy,zz) = ndgrid(x::AbstractVector{<:Any}, y::..., z::...) todo - improve?
source
MIRT.ndgrid — Method
(xx,yy) = ndgrid(x::AbstractVector{<:Any}, y::AbstractVector{<:Any}) todo - improve?
source
MIRT.nufft_errors — Method
w, errs = nufft_errors( ; M=?, w=?, N=?, n_shift=?, ...)
Compute NUFFT approximation errors (for signal of length N of unit norm), for given digital frequency values w, i.e., Ω.
source
MIRT.nufft_init — Method
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.
source
MIRT.nufft_init — Method
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!
+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.
source
MIRT.exp_mult — Method
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)
source
MIRT.exp_xform — Method
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.
source
MIRT.genkspace — Method
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
source
MIRT.genspi — Method
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.
source
MIRT.getindex! — Method
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
source
MIRT.interp1 — Method
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
source
MIRT.ir_dump — Method
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
source
MIRT.ir_mri_coil_compress — Method
(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
source
MIRT.ir_mri_kspace_ga_radial — Method
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"
2015-07 Mai Le, original Matlab version
2015-07-04 Jeff Fessler, minor changes to Matlab version
source
MIRT.ir_mri_sensemap_sim — Method
(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.
source
MIRT.ir_mri_sensemap_sim — Method
(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
source
MIRT.ir_mri_sensemap_sim — Method
smap = ir_mri_sensemap_sim(...)
Simulate 2D or 3D sensitivity maps for sensitivity-encoded MRI based on grivich:00:tmf http://doi.org/10.1119/1.19461
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)
source
MIRT.ir_mri_smap1 — Method
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
source
MIRT.jinc — Method
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).
source
MIRT.map_many — Method
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])
source
MIRT.mask_or — Method
mask_or(mask)
compress 3D mask to 2D by logical or along z direction
source
MIRT.max_percent_diff — Method
d = 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 use max of both arguments to normalize; default false
  • normalize::Bool normalize each before comparing; default false
source
MIRT.mri_kspace_spiral — Method
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?)
source
MIRT.mri_trajectory — Method
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)
source
MIRT.mri_trajectory_gads — Method
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
source
MIRT.mri_trajectory_radial — Method
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
source
MIRT.mri_trajectory_rosette3 — Method
mri_trajectory_rosette3(N, fov ; ...)
3d rosette, with default parameters from bucholz:08:miw
options: omax: maximum omega nt : time samples (65.536 ms for 4 usec dt) dt : time sample spacing (4 usec) ti : time samples
source
MIRT.ncg — Method
(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		array of $J$ blocks $B_1,...,B_J$
  • gradf	array of $J$ functions return gradients of $f_1,...,f_J$
  • curvf	array 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
source
MIRT.ncg — Method
(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
source
MIRT.ndgrid — Method
(xx,yy,zz) = ndgrid(x::AbstractVector{<:Any}, y::..., z::...) todo - improve?
source
MIRT.ndgrid — Method
(xx,yy) = ndgrid(x::AbstractVector{<:Any}, y::AbstractVector{<:Any}) todo - improve?
source
MIRT.nufft_errors — Method
w, errs = nufft_errors( ; M=?, w=?, N=?, n_shift=?, ...)
Compute NUFFT approximation errors (for signal of length N of unit norm), for given digital frequency values w, i.e., Ω.
source
MIRT.nufft_init — Method
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.
source
MIRT.nufft_init — Method
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 an 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`.)
source
MIRT.nufft_typer — Method
nufft_typer(T::DataType, x::AbstractArray{<:Real} ; warn::Bool=true)
type conversion wrapper for nfft()
source
MIRT.ogm_ls — Method
(x,out) = ogm_ls(B, gradf, curvf, x0; niter=?, ninner=?, fun=?)
OGM with a line search; Drori&Taylor @arxiv 1803.05676; 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		array of $J$ blocks $B_1,...,B_J$
  • gradf	array of $J$ functions return gradients of $f_1,...,f_J$
  • curvf	array 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 an array of length niter+1
source
MIRT.ogm_ls — Method
(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
source
MIRT.pogm_restart — Method
x, out = pogm_restart(x0, Fcost, f_grad, f_L ;
+(Using `operator=true` allows the `LinearMapAO` to support `do_many`.)
source
MIRT.nufft_typer — Method
nufft_typer(T::DataType, x::AbstractArray{<:Real} ; warn::Bool=true)
type conversion wrapper for nfft()
source
MIRT.ogm_ls — Method
(x,out) = ogm_ls(B, gradf, curvf, x0; niter=?, ninner=?, fun=?)
OGM with a line search; Drori&Taylor @arxiv 1803.05676; 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		array of $J$ blocks $B_1,...,B_J$
  • gradf	array of $J$ functions return gradients of $f_1,...,f_J$
  • curvf	array 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 an array of length niter+1
source
MIRT.ogm_ls — Method
(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
source
MIRT.pogm_restart — Method
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
source
MIRT.poweriter — Method
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
source
MIRT.rect — Method
rect(x::Real)
Unit width rect function. Potential problem? Bring up with fess.
source
MIRT.reverser — Method
y = reverser(x, dims)
reverse array along specified dimensions (or all if unspecified)
source
MIRT.rmsd100 — Method
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 %
source
MIRT.snr2sigma — Method
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.)
source

+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

source
MIRT.poweriter — Method
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
source
MIRT.rect — Method
rect(x::Real)

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

source
MIRT.reverser — Method
y = reverser(x, dims)

reverse array along specified dimensions (or all if unspecified)

source
MIRT.rmsd100 — Method
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 %
source
MIRT.snr2sigma — Method
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.)

source
diff --git a/dev/search/index.html b/dev/search/index.html index 394bcd56..c41e50b4 100644 --- a/dev/search/index.html +++ b/dev/search/index.html @@ -1,2 +1,2 @@ -Search · MIRT.jl

Loading search...

    +Search · MIRT.jl

    Loading search...

      diff --git a/dev/toc/index.html b/dev/toc/index.html index a99133c3..8aff274e 100644 --- a/dev/toc/index.html +++ b/dev/toc/index.html @@ -1,2 +1,2 @@ -Table of Contents · MIRT.jl
      Edit on GitHub
      +Table of Contents · MIRT.jl
      Edit on GitHub