Add Asense for sensitivity encoded MRI (#121)

* Add Asense

* Test Asense

* Fix test

* v0.17

* Allow smaps to be vector or array

Project.toml

@@ -1,7 +1,7 @@
 name = "MIRT"
 uuid = "7035ae7a-3787-11e9-139a-5545ed3dc201"
 authors = ["fessler <[email protected]>"]
-version = "0.16.0"
+version = "0.17.0"
 
 [deps]
 AVSfldIO = "b6189060-daf9-4c28-845a-cc0984b81781"

src/system/Asense.jl

@@ -0,0 +1,76 @@
+#=
+Asense.jl
+=#
+
+export Asense
+using FFTW: fftshift!, ifftshift!, fft!, bfft!
+
+
+"""
+    Asense(samp, smaps; T)
+
+Construct a MRI encoding matrix 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)`.
+
+Returns a `LinearMapAO` object.
+"""
+function Asense(
+    samp::AbstractArray{<:Bool},
+    smaps::Vector{<:AbstractArray{<:Number}},
+    ;
+    T::DataType = ComplexF32,
+)
+
+    Base.require_one_based_indexing(samp, smaps...)
+    dims = size(samp)
+    ncoil = length(smaps)
+    for ic in 1:ncoil
+        size(smaps[ic]) == dims || throw("size mismatch")
+    end
+
+    N = prod(dims)
+    work1 = Array{T}(undef, dims)
+    work2 = Array{T}(undef, dims)
+    function forw!(y, x)
+        for ic in 1:ncoil
+            @. work1 = x * smaps[ic]
+            fftshift!(work1, fft!(ifftshift!(work2, work1)))
+            y[:,ic] .= work1[samp]
+        end
+    end
+    function back!(x, y)
+        for ic in 1:ncoil
+            embed!(work1, (@view y[:,ic]), samp)
+            fftshift!(work1, bfft!(ifftshift!(work2, work1)))
+            copyto!(work2, smaps[ic])
+            conj!(work2)
+            if ic == 1
+                @. x = work1 * work2
+            else
+                @. x += work1 * work2
+            end
+        end
+    end
+    A = LinearMapAA(forw!, back!, (ncoil*count(samp), N);
+        prop = (name = "Asense", samp, smaps),
+        odim = (count(samp),ncoil), idim=dims, T,
+    )
+    return A
+end
+
+
+function Asense(
+    samp::AbstractArray{<:Bool, D},
+    smaps::AbstractArray{<:Number},
+    ;
+    kwargs...
+) where D
+    ndims(smaps) == D+1 || error("dimension mismatch")
+    smapv = [eachslice(smaps, dims = D+1)...]
+    return Asense(samp, smapv; kwargs...)
+end

src/system/z-list.jl

@@ -1,3 +1,4 @@
 # system/z-list,jl
 
 include("Afft.jl")
+include("Asense.jl")

test/system/Afft.jl

@@ -1,9 +1,9 @@
-# Afft.jl
+# test/Afft.jl
 
 using MIRT: Afft
 
 using LinearMapsAA: LinearMapAM, LinearMapAO
-using Test: @test
+using Test: @test, @testset
 
 
 samp = trues(3,2); samp[2] = false

test/system/Asense.jl

@@ -0,0 +1,29 @@
+# test/Asense.jl
+
+using MIRT: Asense
+using LinearMapsAA: LinearMapAO
+using LinearAlgebra: dot
+using Test: @test, @testset
+
+@testset "Asense" begin
+    dims = (6,5)
+    samp = trues(dims); samp[2] = false
+    T = ComplexF32
+    smaps = [randn(T, dims), randn(T, dims)]
+
+	A = Asense(samp, smaps)
+	@test A isa LinearMapAO
+	@test Matrix(A)' ≈ Matrix(A')
+	@test A.name == "Asense"
+	@test eltype(A) == ComplexF32
+
+    dims = (30,20)
+    samp = rand(dims...) .< 0.5
+    T = ComplexF32
+    ncoil = 2
+    smaps = randn(T, dims..., ncoil)
+	A = Asense(samp, smaps)
+    x = randn(T, dims)
+    y = randn(T, count(samp), ncoil)
+    @test isapprox(dot(y, A * x), dot(A' * y, x))
+end

test/system/z-test.jl

@@ -1,9 +1,10 @@
-# system/z-test.jl
+# test/system/z-test.jl
 
 using Test: @testset
 
 list = [
 "Afft.jl"
+"Asense.jl"
 ]
 
 for file in list