Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

convert3dTO2d #78

Open
aalfi opened this issue Apr 5, 2022 · 10 comments
Open

convert3dTO2d #78

aalfi opened this issue Apr 5, 2022 · 10 comments

Comments

@aalfi
Copy link

aalfi commented Apr 5, 2022

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
Copy link
Contributor

aTrotier commented May 6, 2022

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
Copy link
Contributor

aTrotier commented May 6, 2022

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

  • here : 274x208
  • my data : 220x160

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
Copy link
Contributor

I find the issue in Wavelets.jl. I made a PR : JuliaDSP/Wavelets.jl#79

@aTrotier
Copy link
Contributor

@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.

@JeffFessler
Copy link
Member

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
Copy link
Author

aalfi commented Jun 14, 2022

@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
Copy link
Member

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
Copy link
Contributor

aTrotier commented Jun 15, 2022

@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
Copy link
Member

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
Copy link
Contributor

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 :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants