updating example with cairomakie

[aTrotier]

I am working on the example #211

I have some differences between the reconstruction but I guess their is a lot of changement in the backend and the reg parameters might be tuned ?

[nHackel]

We could also generate those example with Literate.jl, I recently updated the documentation of two packages with that and I found it very helpful. For example, we wouldn't need to keep track of PNGs with git.

I don't know when/if I'll have time to help migrate the examples to Literate.jl though, so it might be easier to stick with the current structure. Or I could show you how I used it with links or maybe a meeting.

I haven't done any MRI reconstruction where I played around with the least squares part. I tried to keep things similar via the tests, but it might not have worked everywhere.

Are the differences something you would expect from just different lambdas? The shape looks fairly different at the upper and lower borders of the image

[aTrotier]

fieldmapReco gives similar results

[aTrotier]

The example I have implemented before are generated with literate : https://github.com/MagneticResonanceImaging/MRIReco.jl/tree/master/docs/lit/examples

Maybe some of the example might be too long to be generated in the CI (the custom example ?)

[aTrotier]

I haven't done any MRI reconstruction where I played around with the least squares part. I tried to keep things similar via the tests, but it might not have worked everywhere.

Are the differences something you would expect from just different lambdas? The shape looks fairly different at the upper and lower borders of the image

I am not able to tune the parameters to reach something close to the initial reconstruction maybe I am missing something

[nHackel]

Ah is the issue with the ADMM + TV Term? I think ρ was changed to be passed in like params[:rho] = 0.1 instead of params[:ρ] = 0.1

[nHackel]

If you set params[:kwargWarning] = true you should see the reconstruction complain about keyword arguments that the solver doesn't understand/filters out.

Since MRIReco does no filtering of its, I made the default hide these messages. When you enable it and see a keyword argument you would expect to have an effect, like rho here, then something changed there

[aTrotier]

Actually it was mostly an issue with the coil sensitivity maps (which explains the artifacts at the bottom of the images)

 Warning: The following arguments were passed but filtered out: reco, ρ, S. Please watch closely if this introduces unexpexted behaviour in your code.

Nice we should add that in the doc

[aTrotier]

For now, I'll first update like that all the example -> then moved to literate

[aTrotier]

senseReco is failing :

[aTrotier]

not sure what happens but it works now

[aTrotier]

I have an issue with the dictionary sparsifying transform

ERROR: MethodError: no method matching (::var"#48#50"{Matrix{…}, Tuple{…}, Tuple{…}, Float64})(::Vector{ComplexF32}, ::Vector{ComplexF32})

Closest candidates are:
  (::var"#48#50")(::Any)
   @ Main ~/Documents/JULIA_PROJECT/MRIReco_fix_subspace/MRIReco.jl/docs/src/examples/exampleCustom.jl:38

Stacktrace:
  [1] prod3!(res::Vector{…}, prod!::var"#48#50"{…}, v::Vector{…}, α::ComplexF32, β::ComplexF32, Mv5::Vector{…})
    @ LinearOperators ~/.julia/packages/LinearOperators/ixrUQ/src/operations.jl:12
  [2] mul!(res::Vector{…}, op::LinearOperator{…}, v::Vector{…}, α::ComplexF32, β::ComplexF32)
    @ LinearOperators ~/.julia/packages/LinearOperators/ixrUQ/src/operations.jl:31
  [3] mul!(res::Vector{…}, op::LinearOperator{…}, v::Vector{…})
    @ LinearOperators ~/.julia/packages/LinearOperators/ixrUQ/src/operations.jl:40
  [4] *(op::LinearOperator{ComplexF32, Int64, var"#48#50"{…}, Nothing, var"#49#51"{…}, Vector{…}}, v::Vector{ComplexF32})
    @ LinearOperators ~/.julia/packages/LinearOperators/ixrUQ/src/operations.jl:47
  [5] prox!(reg::TransformedRegularization{…}, x::Vector{…}, args::Float64)
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/Regularization/TransformedRegularization.jl:29
  [6] iterate(solver::ADMM{…}, state::RegularizedLeastSquares.ADMMState{…})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/ADMM.jl:237
  [7] iterate(solver::ADMM{…})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/RegularizedLeastSquares.jl:190
  [8] iterate
    @ ./iterators.jl:206 [inlined]
  [9] iterate
    @ ./iterators.jl:205 [inlined]
 [10] solve!(solver::ADMM{…}, b::Vector{…}; callbacks::Function, kwargs::@Kwargs{…})
    @ RegularizedLeastSquares ~/.julia/packages/RegularizedLeastSquares/LvIWc/src/RegularizedLeastSquares.jl:111
 [11] reconstruction_simple(acqData::AcquisitionData{…}; reconSize::Tuple{…}, reg::Vector{…}, sparseTrafo::Vector{…}, weights::Vector{…}, solver::Type{…}, encodingOps::Nothing, arrayType::Type{…}, params::@Kwargs{…})
    @ MRIReco ~/Documents/JULIA_PROJECT/MRIReco_fix_subspace/MRIReco.jl/src/Reconstruction/IterativeReconstruction.jl:65
 [12] reconstruction(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
    @ MRIReco ~/Documents/JULIA_PROJECT/MRIReco_fix_subspace/MRIReco.jl/src/Reconstruction/Reconstruction.jl:45
 [13] top-level scope
    @ ~/Documents/JULIA_PROJECT/MRIReco_fix_subspace/MRIReco.jl/docs/src/examples/exampleCustom.jl:79
Some type information was truncated. Use `show(err)` to see complete types.

But that's the only remaining issue for the examples

[nHackel]

The issue is that the LinearOperator expects a two argument function.

To fix it you can do:

function analyzeImage(res, x::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64};t0::Int64=size(D,2),tol=1e-3) where T
  nx,ny = xsize
  px,py = psize
  x = reshape(x,nx,ny)
  x_pad = repeat(x,2,2)[1:nx+px-1,1:ny+py-1] # pad image using periodic boundary conditions
  α = zeros(T,size(D,2),nx,ny)
  patch = zeros(T,px*py)
  for j=1:ny
    for i=1:nx
      patch[:] .= vec(x_pad[i:i+px-1,j:j+py-1])
      norm(patch)==0 && continue
      # matchingpursuit is contained in Wavelets.jl
      α[:,i,j] .= matchingpursuit(patch, x->D*x, x->transpose(D)*x, tol)
    end
  end
  return res .= vec(α)
end

function synthesizeImage(res, α::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64}) where T
  nx,ny = xsize
  px,py = psize
  x = zeros(T,nx+px-1,ny+py-1)
  α = reshape(α,:,nx,ny)
  for j=1:ny
    for i=1:nx
      x[i:i+px-1,j:j+py-1] .+= reshape(D*α[:,i,j],px,py)
    end
  end
  return res .= vec(x[1:nx,1:ny])/(px*py)
end

function dictOp(D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64},tol::Float64=1.e-3) where T
        produ = (res,x)->analyzeImage(x,D,xsize,psize,tol=tol)
        ctprodu = (res,x)->synthesizeImage(x,D,xsize,psize)
  return LinearOperator(T,prod(xsize)*size(D,2),prod(xsize),false,false
          , produ
          , nothing
          , ctprodu )
end

[aTrotier]

Thanks I'll fix that.

Of note, the CI failed sometimes on macOS. The error does not always occurs...

[aTrotier]

I am not able to obtain something working by playing with the hyper-parameters :

I don't know if I should regenerate the dictionary or if it is related to the backend in RegularizedLeastSquares

By the way do you have the code for the dictionary @tknopp, did you only apply KSVD on the fully sampled dataset ?

[tknopp]

By the way do you have the code for the dictionary @tknopp, did you only apply KSVD on the fully sampled dataset ?

No that has been done by @migrosser, who is not working on this anymore. I could not find the scripts.

[aTrotier]

No that has been done by @migrosser, who is not working on this anymore. I could not find the scripts.

Ok, I asked the question because I add to parse the dictionary to ComplexF32 rather than F64 but I don't think it is the issue here.

Patch looks good to me :

[aTrotier]

Reminder to myself : check example -> https://github.com/MagneticResonanceImaging/MRIReco.jl/tree/master/examples/mridataorg