|
1 | | - |
2 | 1 | #= |
3 | 2 | using AbstractNFFTs: AbstractNFFTPlan, convolve!, convolve_transpose! |
| 3 | +using AbstractNFFTs: size_in, size_out |
4 | 4 | using LinearAlgebra: mul! |
5 | 5 | =# |
6 | 6 |
|
7 | | -function sdc(p::AbstractNFFTPlan{T,D,1}; iters=20) where {T,D} |
8 | | - # Weights for sample density compensation. |
9 | | - # Uses method of Pipe & Menon, 1999. Mag Reson Med, 186, 179. |
10 | | - weights = similar(p.tmpVec, Complex{T}, p.J) |
11 | | - weights .= one(Complex{T}) |
12 | | - weights_tmp = similar(weights) |
13 | | - scaling_factor = zero(T) |
| 7 | +""" |
| 8 | + out = _fill_similar(array, v::T, dims) where T |
| 9 | +Return an allocated array similar to `array` filled with value `v`. |
| 10 | +
|
| 11 | +This 2-step initialization helper function |
| 12 | +is needed to accommodate GPU array types. |
| 13 | +The more obvious statement `weights = fill(v, dims)` |
| 14 | +can lead to the wrong array type and cause GPU tests to fail. |
| 15 | +""" |
| 16 | +function _fill_similar(array::AbstractArray, v::T, dims::Union{Integer,Dims}) where T |
| 17 | + weights = similar(array, T, dims) |
| 18 | + fill!(weights, v) |
| 19 | + return weights |
| 20 | +end |
| 21 | + |
| 22 | +#= |
| 23 | +It would be desirable to use the memory pointed to by p.tmpVec |
| 24 | +as working buffer for sdc iterations, |
| 25 | +but this attempt led to intermittent corrupted outputs and errors. |
| 26 | +So it is left commented out for now. |
| 27 | +=# |
| 28 | +#= |
| 29 | +function _reinterpret_real(g::StridedArray{Complex{T}}) where {T <: Real} |
| 30 | + r1 = reinterpret(T, g) |
| 31 | + r2 = @view r1[1:2:end,:] |
| 32 | + return r2 |
| 33 | +end |
| 34 | +=# |
| 35 | + |
| 36 | +""" |
| 37 | + weights = sdc(plan::{Abstract}NFFTPlan; iters=20, ...) |
| 38 | +
|
| 39 | +Compute weights for sample density compensation for given NFFT `plan` |
| 40 | +Uses method of Pipe & Menon, Mag Reson Med, 44(1):179-186, Jan. 1999. |
| 41 | +DOI: 10.1002/(SICI)1522-2594(199901)41:1<179::AID-MRM25>3.0.CO;2-V |
| 42 | +The function applies `iters` iterations of that method. |
| 43 | +
|
| 44 | +The scaling here such that if the plan were 1D with N nodes equispaced |
| 45 | +from -0.5 to 0.5-1/N, then the returned weights are ≈ 1/N. |
| 46 | +
|
| 47 | +The weights are scaled such that ``A' diag(w) A 1_N ≈ 1_N``. |
| 48 | +
|
| 49 | +The returned vector is real, positive values of length `plan.J`. |
| 50 | +
|
| 51 | +This method _almost_ conforms to the `AbstractNFFT` interface |
| 52 | +except that it uses `p.Ñ` and `p.tmpVec` that are not part of that interface. |
| 53 | +
|
| 54 | +There are several named keyword arguments that are work buffers |
| 55 | +that are all mutated: `weights weights_tmp workg workf workv`. |
| 56 | +If the caller provides all of those, |
| 57 | +then this function should make only small allocations. |
| 58 | +""" |
| 59 | +function sdc( |
| 60 | + p::AbstractNFFTPlan{T,D,1}; |
| 61 | + iters::Int = 20, |
| 62 | + # the type of p.tmpVec (if present) is needed for GPU arrays (JLArray): |
| 63 | + _array::AbstractArray = hasfield(typeof(p), :tmpVec) ? |
| 64 | + p.tmpVec : Array{T,D}(undef, zeros(Int, D)...), |
| 65 | + # the following are working buffers that are all mutated: |
| 66 | + weights::AbstractVector{T} = _fill_similar(_array, one(T), only(size_out(p))), |
| 67 | + weights_tmp::AbstractVector = similar(weights), |
| 68 | +# workg::AbstractArray = _reinterpret_real(p.tmpVec), # todo |
| 69 | + workg::AbstractArray = similar(_array, T, p.Ñ), |
| 70 | + workf::AbstractVector = similar(_array, Complex{T}, only(size_out(p))), |
| 71 | + workv::AbstractArray = similar(_array, Complex{T}, size_in(p)), |
| 72 | +) where {T <: Real, D} |
| 73 | + |
| 74 | + return sdc!( |
| 75 | + p, |
| 76 | + iters, |
| 77 | + weights, |
| 78 | + weights_tmp, |
| 79 | + workg, |
| 80 | + workf, |
| 81 | + workv, |
| 82 | + ) |
| 83 | +end |
| 84 | + |
| 85 | + |
| 86 | +""" |
| 87 | + weights = sdc!( p::AbstractNFFTPlan, iters, weights, weights_tmp, workg) |
| 88 | +Compute sampling density compensation `weights` |
| 89 | +without performing any final scaling step. |
| 90 | +
|
| 91 | +Ideally this function should be (nearly) non-allocating. |
| 92 | +""" |
| 93 | +function sdc!( |
| 94 | + p::AbstractNFFTPlan{T,D,1}, |
| 95 | + iters::Int, |
| 96 | + # the following are working buffers that are all mutated: |
| 97 | + weights::AbstractVector{T}, |
| 98 | + weights_tmp::AbstractVector, |
| 99 | + workg::AbstractArray, |
| 100 | +) where {T <: Real, D} |
| 101 | + |
| 102 | + scaling_factor = missing # will be set below |
14 | 103 |
|
15 | 104 | # Pre-weighting to correct non-uniform sample density |
16 | 105 | for i in 1:iters |
17 | | - convolve_transpose!(p, weights, p.tmpVec) |
18 | | - if i==1 |
19 | | - scaling_factor = maximum(abs.(p.tmpVec)) |
| 106 | + convolve_transpose!(p, weights, workg) |
| 107 | + if i == 1 |
| 108 | + scaling_factor = maximum(workg) |
20 | 109 | end |
21 | 110 |
|
22 | | - p.tmpVec ./= scaling_factor |
23 | | - convolve!(p, p.tmpVec, weights_tmp) |
| 111 | + workg ./= scaling_factor |
| 112 | + convolve!(p, workg, weights_tmp) |
24 | 113 | weights_tmp ./= scaling_factor |
25 | | - weights ./= (abs.(weights_tmp) .+ eps(T)) |
| 114 | + any(≤(0), weights_tmp) && throw("non-positive weights") |
| 115 | + weights ./= weights_tmp |
26 | 116 | end |
27 | | - # Post weights to correct image scaling |
28 | | - # This finds c, where ||u - c*v||_2^2 = 0 and then uses |
29 | | - # c to scale all weights by a scalar factor. |
30 | | - u = similar(weights, Complex{T}, p.N) |
31 | | - u .= one(Complex{T}) |
32 | | - |
33 | | - # conversion to Array is a workaround for CuNFFT. Without it we get strange |
34 | | - # results that indicate some synchronization issue |
35 | | - #f = Array( p * u ) |
36 | | - #b = f .* Array(weights) # apply weights from above |
37 | | - #v = Array( adjoint(p) * convert(typeof(weights), b) ) |
38 | | - #c = vec(v) \ vec(Array(u)) # least squares diff |
39 | | - #return abs.(convert(typeof(weights), c * Array(weights))) |
40 | | - |
41 | | - # non converting version |
42 | | - f = similar(p.tmpVec, Complex{T}, p.J) |
43 | | - mul!(f, p, u) |
44 | | - f .*= weights # apply weights from above |
45 | | - v = similar(p.tmpVec, Complex{T}, p.N) |
46 | | - mul!(v, adjoint(p), f) |
47 | | - c = vec(v) \ vec(u) # least squares diff |
48 | | - return abs.(c * weights) |
| 117 | + return weights |
| 118 | +end |
| 119 | + |
| 120 | + |
| 121 | +""" |
| 122 | + sdc!( ... ) |
| 123 | +
|
| 124 | +This version scales the weights such that |
| 125 | +``A' D(w) A 1_N ≈ 1_N``. |
| 126 | +Find ``c ∈ ℝ`` that minimizes ``‖u - c * v‖₂`` |
| 127 | +where ``u ≜ 1_N`` (vector of ones) |
| 128 | +and ``v ≜ A' D(w) A 1_N``, and then scale `w` by that `c`. |
| 129 | +The analytical solution is |
| 130 | +``c = real(u'v) / ‖v‖² = real(sum(v)) / ‖v‖².`` |
| 131 | +
|
| 132 | +Ideally this function should be (nearly) non-allocating. |
| 133 | +""" |
| 134 | +function sdc!( |
| 135 | + p::AbstractNFFTPlan{T,D,1}, |
| 136 | + iters::Int, |
| 137 | + # the following are working buffers that are all mutated: |
| 138 | + weights::AbstractVector{T}, |
| 139 | + weights_tmp::AbstractVector, |
| 140 | + workg::AbstractArray, |
| 141 | + workf::AbstractVector, |
| 142 | + workv::AbstractArray, |
| 143 | +) where {T <: Real, D} |
| 144 | + |
| 145 | + sdc!(p, iters, weights, weights_tmp, workg) |
| 146 | + |
| 147 | + u = workv # trick to save memory |
| 148 | + fill!(u, one(T)) |
| 149 | + mul!(workf, p, u) |
| 150 | + workf .*= weights # apply weights from above |
| 151 | + mul!(workv, adjoint(p), workf) |
| 152 | + c = real(sum(workv)) / sum(abs2, workv) |
| 153 | + weights .*= c |
| 154 | + return weights |
49 | 155 | end |
0 commit comments