search

MagneticResonanceImaging / MRIReco.jl

Julia Package for MRI Reconstruction
https://magneticresonanceimaging.github.io/MRIReco.jl/latest/
Other
85 stars 22 forks source link

convert3dTO2d #78

Open aalfi opened 2 years ago

aalfi commented 2 years ago

I'm trying to run the example.jl .. I am a newbie learning Julia The error mentioned below acqData2d = convert3dTO2d(acqData)

ERROR: LoadError: MethodError: no method matching CartesianTrajectory(::Int64, ::Int64; TE=0.0f0, AQ=0.001f0) Closest candidates are: CartesianTrajectory(::Type{T}, ::Any, ::Any; TE, AQ, kmin, kmax, kargs...) where T at ~/.julia/packages/MRIBase/it3vS/src/Trajectories/2D/Cartesian2D.jl:22 Stacktrace: [1] convert3dTo2d(acqData::AcquisitionData{Float32}) @ MRIBase ~/.julia/packages/MRIBase/it3vS/src/Datatypes/AcqData.jl:314 [2] top-level scope @ ~/julia/MRIReco.jl-master/example.jl:26 [3] include(fname::String) @ Base.MainInclude ./client.jl:451 [4] top-level scope @ REPL[13]:1 [5] top-level scope @ ~/.julia/packages/CUDA/5jdFl/src/initialization.jl:52

aTrotier commented 2 years ago

We have made some modification to the package which broke the example.jl.

I have a branch to work on the correction : https://github.com/MagneticResonanceImaging/MRIReco.jl/tree/Example

The function convert3dTO2d is now working but I have an other error during the reconstruction which required to change the lambda to params[:λ] = Float32(2.e-1)

julia> reconstruction(acqData2d, params)
ERROR: TaskFailedException

    nested task error: MethodError: no method matching proxL1!(::Vector{ComplexF32}, ::Float64; λ=0.2, reco="multiCoil", absTol=0.0001, normalizeReg=true, sparseTrafo=Linear operator
      nrow: 56992
      ncol: 56992
      eltype: ComplexF64
      symmetric: false
      hermitian: false
      nprod:   0
      ntprod:  0
      nctprod: 0

    , relTol=0.001, tolInner=0.01, ρ=0.1, iterations=30, solver="admm", senseMaps=[0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im;;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;;; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im;;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;;; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;;; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; … ; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im; -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im … -0.0f0 - 0.0f0im -0.0f0 - 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; … ; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im; 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im … 0.0f0 + 0.0f0im 0.0f0 + 0.0f0im;;;; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; … ; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im; 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im … 0.0f0 - 0.0f0im 0.0f0 - 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im;;; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; … ; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im; -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im … -0.0f0 + 0.0f0im -0.0f0 + 0.0f0im], reconSize=(274, 208), regularization="L1", sparseTrafoName="Wavelet", shape=(274, 208))
    Closest candidates are:
      proxL1!(::AbstractArray{Tc}, ::T; sparseTrafo, kargs...) where {T, Tc<:Union{Complex{T}, T}} at ~/.julia/packages/RegularizedLeastSquares/37BYe/src/proximalMaps/ProxL1.jl:13
    Stacktrace:
     [1] iterate(solver::RegularizedLeastSquares.ADMM{Float32, MRIReco.CompositeOp{ComplexF32}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#14#17"{LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, MRIReco.CompositeOp{ComplexF32}, Vector{ComplexF32}}, LinearOperators.var"#15#18"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.var"#16#19"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, Vector{ComplexF32}}, Vector{ComplexF32}, Vector{Float32}, IterativeSolvers.Identity}, iteration::Int64)
       @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/37BYe/src/ADMM.jl:239
     [2] iterate
       @ ~/.julia/packages/RegularizedLeastSquares/37BYe/src/ADMM.jl:220 [inlined]
     [3] iterate
       @ ./iterators.jl:159 [inlined]
     [4] iterate
       @ ./iterators.jl:158 [inlined]
     [5] solve(solver::RegularizedLeastSquares.ADMM{Float32, MRIReco.CompositeOp{ComplexF32}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#14#17"{LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, MRIReco.CompositeOp{ComplexF32}, Vector{ComplexF32}}, LinearOperators.var"#15#18"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.var"#16#19"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, Vector{ComplexF32}}, Vector{ComplexF32}, Vector{Float32}, IterativeSolvers.Identity}, b::Vector{ComplexF32}; A::MRIReco.CompositeOp{ComplexF32}, AᴴA::LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#14#17"{LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, MRIReco.CompositeOp{ComplexF32}, Vector{ComplexF32}}, LinearOperators.var"#15#18"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.var"#16#19"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, Vector{ComplexF32}}, startVector::Vector{ComplexF32}, solverInfo::Nothing, kargs::Base.Pairs{Symbol, Any, NTuple{14, Symbol}, NamedTuple{(:absTol, :normalizeReg, :sparseTrafo, :relTol, :tolInner, :solver, :senseMaps, :regularization, :sparseTrafoName, :λ, :reco, :ρ, :iterations, :reconSize), Tuple{Float64, Bool, String, Float64, Float64, String, Array{ComplexF32, 4}, String, String, Float64, String, Float64, Int64, Tuple{Int64, Int64}}}})
       @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/37BYe/src/ADMM.jl:207
     [6] macro expansion
       @ ~/.julia/dev/MRIReco/src/Reconstruction/IterativeReconstruction.jl:192 [inlined]
     [7] (::MRIReco.var"#201#202"{AcquisitionData{Float32}, Tuple{Int64, Int64}, Vector{Regularization}, Vector{Vector{ComplexF32}}, String, Array{ComplexF32, 4}, Nothing, Dict{Symbol, Any}, Dict{Any, Any}, Int64, Int64, Int64})()
       @ MRIReco ./threadingconstructs.jl:178

...and 3 more exceptions.

Stacktrace:
 [1] sync_end(c::Channel{Any})
   @ Base ./task.jl:381
 [2] macro expansion
   @ ./task.jl:400 [inlined]
 [3] reconstruction_multiCoil(acqData::AcquisitionData{Float32}, reconSize::Tuple{Int64, Int64}, reg::Vector{Regularization}, sparseTrafo::LinearOperators.LinearOperator{ComplexF64, Int64, SparsityOperators.var"#1#2"{SparsityOperators.var"#28#30"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Nothing, SparsityOperators.var"#1#2"{SparsityOperators.var"#29#31"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Vector{ComplexF64}}, weights::Vector{Vector{ComplexF32}}, solvername::String, senseMaps::Array{ComplexF32, 4}, normalize::Bool, encodingOps::Nothing, params::Dict{Symbol, Any})
   @ MRIReco ~/.julia/dev/MRIReco/src/Reconstruction/IterativeReconstruction.jl:178
 [4] reconstruction(acqData::AcquisitionData{Float32}, recoParams::Dict{Symbol, Any})
   @ MRIReco ~/.julia/dev/MRIReco/src/Reconstruction/Reconstruction.jl:39
 [5] top-level scope
   @ REPL[20]:1

I have more error after that again

@tknopp Maybe we should correct that and then include it as a test ?

aTrotier commented 2 years ago

Next error is related to the Wavelet transform. I also experiment that on one of my datasets for specific rectangular matrix size :

Should we increase the image size in order to be squared in the linearOperator (and cut that part during the adjoint operation ?)

julia> reconstruction(acqData2d, params)
ERROR: TaskFailedException

    nested task error: ArgumentError: size must have a sufficient power of 2 factor
    Stacktrace:
      [1] _dwt!(y::Matrix{ComplexF32}, x::Matrix{ComplexF32}, filter::Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}, L::Int64, fw::Bool, dcfilter::Vector{ComplexF32}, scfilter::Vector{ComplexF32}, si::Vector{ComplexF32}, tmpbuffer::Vector{ComplexF32})
        @ Wavelets.Transforms ~/.julia/packages/Wavelets/RrFaf/src/mod/transforms_filter.jl:133
      [2] _dwt!(y::Matrix{ComplexF32}, x::Matrix{ComplexF32}, filter::Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}, L::Int64, fw::Bool)
        @ Wavelets.Transforms ~/.julia/packages/Wavelets/RrFaf/src/mod/transforms_filter.jl:121
      [3] dwt
        @ ~/.julia/packages/Wavelets/RrFaf/src/mod/Transforms.jl:117 [inlined]
      [4] dwt
        @ ~/.julia/packages/Wavelets/RrFaf/src/mod/Transforms.jl:116 [inlined]
      [5] (::SparsityOperators.var"#28#30"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}})(x::Vector{ComplexF32})
        @ SparsityOperators ~/.julia/packages/SparsityOperators/eHu9U/src/WaveletOp.jl:16
      [6] (::SparsityOperators.var"#1#2"{SparsityOperators.var"#28#30"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}})(res::Vector{ComplexF64}, x::Vector{ComplexF32}, α::ComplexF32, β::ComplexF32)
        @ SparsityOperators ~/.julia/packages/SparsityOperators/eHu9U/src/SparsityOperators.jl:0
      [7] mul!(res::Vector{ComplexF64}, op::LinearOperators.LinearOperator{ComplexF64, Int64, SparsityOperators.var"#1#2"{SparsityOperators.var"#28#30"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Nothing, SparsityOperators.var"#1#2"{SparsityOperators.var"#29#31"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Vector{ComplexF64}}, v::Vector{ComplexF32}, α::ComplexF32, β::ComplexF32)
        @ LinearOperators ~/.julia/packages/LinearOperators/rejMT/src/operations.jl:29
      [8] mul!
        @ ~/.julia/packages/LinearOperators/rejMT/src/operations.jl:36 [inlined]
      [9] *
        @ ~/.julia/packages/LinearOperators/rejMT/src/operations.jl:43 [inlined]
     [10] proxL1!(x::Vector{ComplexF32}, λ::Float32; sparseTrafo::LinearOperators.LinearOperator{ComplexF64, Int64, SparsityOperators.var"#1#2"{SparsityOperators.var"#28#30"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Nothing, SparsityOperators.var"#1#2"{SparsityOperators.var"#29#31"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Vector{ComplexF64}}, kargs::Base.Pairs{Symbol, Any, NTuple{14, Symbol}, NamedTuple{(:λ, :reco, :absTol, :normalizeReg, :relTol, :tolInner, :ρ, :iterations, :solver, :senseMaps, :reconSize, :regularization, :sparseTrafoName, :shape), Tuple{Float32, String, Float64, Bool, Float64, Float64, Float64, Int64, String, Array{ComplexF32, 4}, Tuple{Int64, Int64}, String, String, Tuple{Int64, Int64}}}})
        @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/37BYe/src/proximalMaps/ProxL1.jl:17
     [11] iterate(solver::RegularizedLeastSquares.ADMM{Float32, MRIReco.CompositeOp{ComplexF32}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#14#17"{LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, MRIReco.CompositeOp{ComplexF32}, Vector{ComplexF32}}, LinearOperators.var"#15#18"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.var"#16#19"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, Vector{ComplexF32}}, Vector{ComplexF32}, Vector{Float32}, IterativeSolvers.Identity}, iteration::Int64)
        @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/37BYe/src/ADMM.jl:239
     [12] iterate
        @ ~/.julia/packages/RegularizedLeastSquares/37BYe/src/ADMM.jl:220 [inlined]
     [13] iterate
        @ ./iterators.jl:159 [inlined]
     [14] iterate
        @ ./iterators.jl:158 [inlined]
     [15] solve(solver::RegularizedLeastSquares.ADMM{Float32, MRIReco.CompositeOp{ComplexF32}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#14#17"{LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, MRIReco.CompositeOp{ComplexF32}, Vector{ComplexF32}}, LinearOperators.var"#15#18"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.var"#16#19"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, LinearOperators.var"#112#113"{Int64}, Vector{ComplexF32}}, Vector{ComplexF32}, Vector{Float32}, IterativeSolvers.Identity}, b::Vector{ComplexF32}; A::MRIReco.CompositeOp{ComplexF32}, AᴴA::LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#14#17"{LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, MRIReco.CompositeOp{ComplexF32}, Vector{ComplexF32}}, LinearOperators.var"#15#18"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, LinearOperators.var"#16#19"{MRIReco.CompositeOp{ComplexF32}, LinearOperators.AdjointLinearOperator{ComplexF32, MRIReco.CompositeOp{ComplexF32}}, Vector{ComplexF32}}, Vector{ComplexF32}}, startVector::Vector{ComplexF32}, solverInfo::Nothing, kargs::Base.Pairs{Symbol, Any, NTuple{14, Symbol}, NamedTuple{(:absTol, :normalizeReg, :sparseTrafo, :relTol, :tolInner, :solver, :senseMaps, :regularization, :sparseTrafoName, :λ, :reco, :ρ, :iterations, :reconSize), Tuple{Float64, Bool, String, Float64, Float64, String, Array{ComplexF32, 4}, String, String, Float32, String, Float64, Int64, Tuple{Int64, Int64}}}})
        @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/37BYe/src/ADMM.jl:207
     [16] macro expansion
        @ ~/.julia/dev/MRIReco/src/Reconstruction/IterativeReconstruction.jl:192 [inlined]
     [17] (::MRIReco.var"#201#202"{AcquisitionData{Float32}, Tuple{Int64, Int64}, Vector{Regularization}, Vector{Vector{ComplexF32}}, String, Array{ComplexF32, 4}, Nothing, Dict{Symbol, Any}, Dict{Any, Any}, Int64, Int64, Int64})()
        @ MRIReco ./threadingconstructs.jl:178

...and 3 more exceptions.

Stacktrace:
 [1] sync_end(c::Channel{Any})
   @ Base ./task.jl:381
 [2] macro expansion
   @ ./task.jl:400 [inlined]
 [3] reconstruction_multiCoil(acqData::AcquisitionData{Float32}, reconSize::Tuple{Int64, Int64}, reg::Vector{Regularization}, sparseTrafo::LinearOperators.LinearOperator{ComplexF64, Int64, SparsityOperators.var"#1#2"{SparsityOperators.var"#28#30"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Nothing, SparsityOperators.var"#1#2"{SparsityOperators.var"#29#31"{Tuple{Int64, Int64}, Wavelets.WT.OrthoFilter{Wavelets.WT.PerBoundary}}}, Vector{ComplexF64}}, weights::Vector{Vector{ComplexF32}}, solvername::String, senseMaps::Array{ComplexF32, 4}, normalize::Bool, encodingOps::Nothing, params::Dict{Symbol, Any})
   @ MRIReco ~/.julia/dev/MRIReco/src/Reconstruction/IterativeReconstruction.jl:178
 [4] reconstruction(acqData::AcquisitionData{Float32}, recoParams::Dict{Symbol, Any})
   @ MRIReco ~/.julia/dev/MRIReco/src/Reconstruction/Reconstruction.jl:39
 [5] top-level scope
   @ REPL[20]:1
aTrotier commented 2 years ago

I find the issue in Wavelets.jl. I made a PR : https://github.com/JuliaDSP/Wavelets.jl/pull/79

aTrotier commented 2 years ago

@aalfi the PR has been merged. I think you should be able to reconstruct the "non-squared" image now :)

However in that case (220x160) we only have 2 levels of wavelets transform, I don't know if it will be enough to get a good reconstruction with CS.

It is possible to zero filled the image to the next 2^n size (see https://github.com/tknopp/SparsityOperators.jl/issues/12) . I am waiting for @tknopp, maybe @JeffFessler also have some thoughts about that.

JeffFessler commented 2 years ago

I don't have any experience with zero-padding prior to a wavelet transform. The DWT is a unitary transform that simplifies proximal operators, whereas the product of the DTW with a zero-padding matrix is not unitary so that seems likely to complicate fast proximal optimization methods, but for splitting methods like ADMM it would be no problem. Zero-padding of axial slices that have air all around seems benign. But for coronal or sag. slices (or in the z dimension for 3D imaging) I would not recommend zero padding because of possible discontinuities at the top/bottom slices, unless a slab excitation is used that makes those end slices close to zero anyway.

aalfi commented 2 years ago

@aalfi the PR has been merged. I think you should be able to reconstruct the "non-squared" image now :)

However in that case (220x160) we only have 2 levels of wavelets transform, I don't know if it will be enough to get a good reconstruction with CS.

It is possible to zero filled the image to the next 2^n size (see tknopp/SparsityOperators.jl#12) . I am waiting for @tknopp, maybe @JeffFessler also have some thoughts about that.

Thank you so much!

JeffFessler commented 2 years ago

Addendum. The comment I wrote above assumed that you were thinking of zero padding only for the regularizer, i.e., a regularizer of the form || DTW ZeroPad x ||. Perhaps you meant instead to reconstruct a larger sized image that has more factors of 2 in its dimensions, so the "zero padding" takes place just once at the start of the iterative recon. That would be a very reasonable approach for non-Cartesian k-space sampling because the NUFFT interpolator can handle any spacing. If you have Cartesian k-space sampling then you would want to think a bit more about the data-fit term. I think it could be done efficiently with some small modifications to existing algorithms. Looks like the example uses ADMM and it should be do-able to modify ADMM to accommodate such "padding" for Cartesian sampling.

aTrotier commented 2 years ago

@aalfi let us now if it works :)

@JeffFessler thanks for the answer I was thinking about Zero padding only the regularizer. Something like that :

"""
    WaveletOp_custom(T::Type, shape, wt=wavelet(WT.db2))

    Compute a wavelet operator that zerofill the image to the next 2^n squared matrix

    Return the Wavelet Operator and the size of the zero filled matrix
"""
function WaveletOp_custom(T::Type, shape; wt=wavelet(WT.db2),L=nothing)
    shape_z = maximum(shape)

    n = 0
    while (2^n < shape_z)
        n = n + 1
    end

    sizeN = 2^n
    shape_z = (sizeN,sizeN)

    Lmax = minimum(Wavelets.Util.maxtransformlevels.(shape_z))

    if L === nothing 
        L = Lmax
    else
        (L > Lmax) ? error("L=$L > maximum possible wavelet transform $Lmax") : L
    end

    @info "Number of wavelet transform = $L and matrix = $sizeN"

    function zeroS(ima)
        ima = reshape(ima,shape)
        ima_z = zeros(typeof(ima[1]),shape_z)
        ima_z[1:shape[1],1:shape[2]] .= ima
        return ima_z
    end

    function cropS(ima_z)
        ima = ima_z[1:shape[1],1:shape[2]]
        return ima
    end

    function EH(W_ima)
        ima = vec(cropS(idwt(reshape(W_ima,shape_z),wt,L)))
        return ima
    end

    function E(ima)
        W_ima = vec(dwt(zeroS(reshape(ima,shape)),wt,L))
        return W_ima
    end

    return LinearOperator{T}(prod(shape_z), prod(shape), false, false
              , wrapProd( x->E(x))
              , nothing
              , wrapProd( y->EH(y))) , shape_z
end

My guess was that the level of DWT has an impact on the CS reconstruction (with only 2 levels of DWT we are less sparse)

My idea was to not restrict the level of DWT to 2 due to the "weird" matrix size along one dimension. Hopefully, most of the 2D matrix size let us reach more level. But for 3D acquisitions, it might be more common to use a matrix which leads to a level = 2/3.

JeffFessler commented 2 years ago

That operator is potentially useful but it is not unitary so you probably would need to modify the optimization algorithm. BTW, nextpow is useful here.

aTrotier commented 2 years ago

BTW, I asked the BART community about their implementation of the DWT Martin answer :

it will do more levels for the larger dimensions.

And in the sigpy example : https://github.com/mikgroup/sigpy-mri-tutorial/blob/master/03-building-an-l1-wavelet-recon-app.ipynb

It seems the DWT is using some signal extension (I guess it is zero-padding but the wavelet module have multiple options image

image

@JeffFessler Thanks, I will leave the code like that for now and maybe discuss that during one of the virtual meeting :)