Building SENSE operators

[JAgho]

I'm working on self-navigated multi-shot DW FSE data, and this entails reconstruction of each of the (aliased) shots to compute phase correction maps for each shot.

As far as I can tell this should be possible using the operator / solver interface here, but I'm getting a "shape mismatch" error which I'm struggling to fix. I've precomputed the SENSE maps using BART and they look reasonable. I based my attempt on the encodingOps_parallel present in the code, and as far as I can tell I'm presenting data in the right format

Given a single shot of data with 1440 k-space points across 8 channels, we have

#make some dummy data
data = rand(ComplexF32, 1440, 8)
trajectory = rand(3, 1440).- 0.5
matrix_size = (74, 66, 60)
b0_sense = rand(ComplexF32, matrix_size..., 8)

#build some operators
ft = LinearOperatorCollection.NFFTOp(ComplexF32, shape=matrix_size, nodes=trajectory, S=Vector{ComplexF32})
SOp = SensitivityOp(b0_sense)
forward_op = DiagOp(ft ; copyOpsFn=copy) ∘ SOp

#make a solver
linear_solver = createLinearSolver(CGNR, forward_op)

#(hopefully) solve
I = solve!(linear_solver, data_trans)

This fails with the following stacktrace:

ERROR: LinearOperatorException("shape mismatch")
Stacktrace:
  [1] mul!(res::Vector{…}, op::AdjointLinearOperator{…}, v::Vector{…}, α::ComplexF32, β::ComplexF32)
    @ LinearOperators ~/.julia/packages/LinearOperators/ixrUQ/src/adjtrans.jl:107
  [2] mul!
    @ ~/.julia/packages/LinearOperators/ixrUQ/src/operations.jl:40 [inlined]
  [3] ctprodu!
    @ ~/.julia/packages/LinearOperatorCollection/GRBTA/src/ProdOp.jl:53 [inlined]
  [4] #18
    @ ~/.julia/packages/LinearOperatorCollection/GRBTA/src/ProdOp.jl:60 [inlined]
  [5] prod3!(res::Vector{…}, prod!::LinearOperatorCollection.var"#18#24"{…}, v::Vector{…}, α::ComplexF32, β::ComplexF32, Mv5::Vector{…})
    @ LinearOperators ~/.julia/packages/LinearOperators/ixrUQ/src/operations.jl:12
  [6] mul!(res::Vector{…}, op::AdjointLinearOperator{…}, v::Vector{…}, α::ComplexF32, β::ComplexF32)
    @ LinearOperators ~/.julia/packages/LinearOperators/ixrUQ/src/adjtrans.jl:117
  [7] mul!
    @ ~/.julia/packages/LinearOperators/ixrUQ/src/operations.jl:40 [inlined]
  [8] initCGNR
    @ ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/CGNR.jl:132 [inlined]
  [9] init!(solver::CGNR{…}, state::RegularizedLeastSquares.CGNRState{…}, b::Vector{…}; x0::Int64)
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/CGNR.jl:123
 [10] init!(solver::CGNR{…}, state::RegularizedLeastSquares.CGNRState{…}, b::Vector{…})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/CGNR.jl:107
 [11] init!(solver::CGNR{…}, state::RegularizedLeastSquares.CGNRState{…}, b::Matrix{…}; scheduler::Type{…}, kwargs::@Kwargs{})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/MultiThreading.jl:35
 [12] init!(solver::CGNR{…}, state::RegularizedLeastSquares.CGNRState{…}, b::Matrix{…})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/MultiThreading.jl:30
 [13] init!(solver::CGNR{…}, b::Matrix{…}; kwargs::@Kwargs{})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/RegularizedLeastSquares.jl:189
 [14] init!
    @ ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/RegularizedLeastSquares.jl:189 [inlined]
 [15] solve!(solver::CGNR{…}, b::Matrix{…}; callbacks::Function, kwargs::@Kwargs{})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/RegularizedLeastSquares.jl:108
 [16] solve!(solver::CGNR{ProdOp{…}, LinearOperator{…}, L2Regularization{…}, Vector{…}}, b::Matrix{ComplexF32})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/RegularizedLeastSquares.jl:102
 [17] top-level scope
    @ ~/.julia/dev/RefocusingRecon/src/load_data.jl:191
Some type information was truncated. Use `show(err)` to see complete types.

As far as I can tell it's just something about the operator interface that I'm not understanding, but if anyone knows what mistake I'm making, it'd be a huge help. Thanks!

[alexjaffray]

The k-space data is expected to be a vector of length (samples x channels), not a 2D array in your case. You can call size() on your operators to make sure they have the correct dimensionality to commute with your data and each other.