updating example with cairomakie

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 ImageUtils = "8ad4436d-4835-5a14-8bce-3ae014d2950b"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

docs/src/API.md

@@ -9,27 +9,20 @@ Operators are implemented as subtypes of `AbstractLinearOperator`, which is defi
 
 
 ```@docs
-MRIReco.encodingOps2d_simple
-MRIReco.encodingOps3d_simple
-MRIReco.encodingOps2d_parallel
-MRIReco.encodingOps3d_parallel
-MRIReco.encodingOp2d_multiEcho
-MRIReco.encodingOp3d_multiEcho
-MRIReco.encodingOp2d_multiEcho_parallel
-MRIReco.encodingOp3d_multiEcho_parallel
-MRIReco.fourierEncodingOp2d
-MRIReco.fourierEncodingOp3d
-MRIReco.ExplicitOp(shape::NTuple{D,Int64}, tr::Trajectory, correctionmap::Array{ComplexF64,D}; MRIReco.echoImage::Bool=false, kargs...) where D
-RegularizedLeastSquares.FFTOp(T::Type, shape::Tuple, shift=true)
-MRIReco.NFFTOp(shape::Tuple, tr::Trajectory; nodes=nothing, kargs...)
-MRIReco.FieldmapNFFTOp(shape::NTuple{D,Int64}, tr::Trajectory,
-                        correctionmap::Array{ComplexF64,D};
-                        method::String="nfft",
-                        echoImage::Bool=true,
-                        alpha::Float64=1.75,
-                        m::Float64=3.0,
-                        K=20,
-                        kargs...) where D
+MRIReco.EncodingOp
+MRIReco.lrEncodingOp
+MRIReco.fourierEncodingOp
+MRIReco.encodingOps_simple
+MRIReco.encodingOps_parallel
+MRIReco.encodingOp_multiEcho
+MRIReco.encodingOp_multiEcho_parallel
+MRIReco.fourierEncodingOp
+MRIReco.ExplicitOp(shape::NTuple{D,Int64}, tr::Trajectory{T}, correctionmap::AbstractArray{Tc,D}
+                        ; echoImage::Bool=false, S = storage_type(correctionmap), kargs...) where {T, Tc <: Union{Complex{T}, T}, D}
+MRIReco.SubspaceOp
+MRIReco.FFTop                   
+MRIReco.NFFTOp
+MRIReco.FieldmapNFFTOp
 MRIReco.SamplingOp
 MRIReco.SensitivityOp
 MRIReco.SparseOp

docs/src/examples/exampleCS.jl

@@ -1,5 +1,6 @@
-using PyPlot, MRIReco, MRIReco.RegularizedLeastSquares
-
+using CairoMakie, MRIReco, MRIReco.RegularizedLeastSquares
+using MRIFiles, MRISampling, MRICoilSensitivities
+include(joinpath(@__DIR__,"exampleUtils.jl"))
 # load fully sampled data
 f = ISMRMRDFile(@__DIR__()*"/data/knee_3dFSE_slice170.h5")
 acqData = AcquisitionData(f);
@@ -24,7 +25,7 @@ msk[acqDataSub.subsampleIndices[1]] .= 1
 # Estimate the coil sensitivities using ESPIRiT and reconstruct using SENSE
 
 # coil sensitivities
-smaps = espirit(acqData,(6,6),30,eigThresh_1=0.035,eigThresh_2=0.98)
+smaps = espirit(acqData,(6,6),30,eigThresh_2=0)
 
 # SENSE reconstruction
 params = Dict{Symbol, Any}()
@@ -47,32 +48,33 @@ params[:reconSize] = (320,320)
 params[:senseMaps] = smaps
 
 params[:solver] = ADMM
-params[:reg] = TVRegularization(1.e-1, shape = (320, 320))
+params[:reg] = TVRegularization(2e-1, shape = (320, 320))
 params[:iterations] = 50
-params[:ρ] = 0.1
+params[:rho] = 0.1
 params[:absTol] = 1.e-4
 params[:relTol] = 1.e-2
 params[:tolInner] = 1.e-2
 params[:normalizeReg] = MeasurementBasedNormalization()
 
+
 img_tv = reconstruction(acqDataSub, params)
 
 
-# use PyPlot for interactive display
-figure(1)
-clf()
-subplot(2,2,1)
-title("Phantom")
-imshow(abs.(img[:,:,1,1,1]))
-subplot(2,2,2)
-title("Mask")
-imshow(abs.(msk))
-subplot(2,2,3)
-title("CG Reconstruction")
-imshow(abs.(img_cg[:,:,1,1,1]))
-subplot(2,2,4)
-title("TV Reconstruction")
-imshow(abs.(img_tv[:,:,1,1,1]))
+begin
+  colormap=:grays
+  f = Figure(size=(800,800))
+  ax = Axis(f[1,1],title="Phantom")
+  heatmap!(ax,rotr90(abs.(img[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,2],title="Mask")
+  heatmap!(ax,rotr90(abs.(msk));colormap)
+  ax = Axis(f[2,1],title="CG Reconstruction")
+  heatmap!(ax,rotr90(abs.(img_cg[:,:,1,1,1]));colormap)
+  ax = Axis(f[2,2],title="TV Reconstruction")
+  heatmap!(ax,rotr90(abs.(img_tv[:,:,1,1,1]));colormap)
+  [hidedecorations!(f.content[ax]) for ax in eachindex(f.content)]
+  f
+end
+
 
 # export images
 filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/kneeOrig.png")
@@ -82,4 +84,4 @@ exportImage(filename, abs.(msk[:,:,1,1,1]) )
 filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/kneeCG.png")
 exportImage(filename, abs.(img_cg[:,:,1,1,1]) )
 filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/kneeTV.png")
-exportImage(filename, abs.(img_tv[:,:,1,1,1]) )
+exportImage(filename, abs.(img_tv[:,:,1,1,1]) )

docs/src/examples/exampleCustom.jl

@@ -1,6 +1,9 @@
-using Wavelets, DelimitedFiles, LinearAlgebra, PyPlot, MRIReco, MRIReco.RegularizedLeastSquares
+using MRIReco, MRIFiles, MRIOperators, MRISampling
+using MRIReco.MRIOperators.Wavelets, MRIReco.RegularizedLeastSquares
+using DelimitedFiles, LinearAlgebra, CairoMakie
+include(joinpath(@__DIR__,"exampleUtils.jl"))
 
-function analyzeImage(x::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64};t0::Int64=size(D,2),tol=1e-3) where T
+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)
@@ -15,10 +18,10 @@ function analyzeImage(x::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize::NT
       α[:,i,j] .= matchingpursuit(patch, x->D*x, x->transpose(D)*x, tol)
     end
   end
-  return vec(α)
+  return res .= vec(α)
 end
 
-function synthesizeImage(α::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize::NTuple{2,Int64}) where T
+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)
@@ -28,34 +31,52 @@ function synthesizeImage(α::Vector{T},D::Matrix{T},xsize::NTuple{2,Int64},psize
       x[i:i+px-1,j:j+py-1] .+= reshape(D*α[:,i,j],px,py)
     end
   end
-  return vec(x[1:nx,1:ny])/(px*py)
+  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 = x->analyzeImage(x,D,xsize,psize,tol=tol)
-  ctprodu = x->synthesizeImage(x,D,xsize,psize)
-  return LinearOperator(prod(xsize)*size(D,2),prod(xsize),false,false
+        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
-
 ## To test our method, let us load some simulated data and subsample it
 
 # phantom
-img = readdlm("data/mribrain100.tsv")
+img = readdlm(joinpath(@__DIR__,"data/mribrain100.tsv"))
 
-acqData = AcquisitionData(ISMRMRDFile("data/acqDataBrainSim100.h5"))
+acqData = AcquisitionData(ISMRMRDFile(joinpath(@__DIR__,"data/acqDataBrainSim100.h5")))
 nx,ny = acqData.encodingSize
 
 # undersample kspace data
 acqData = sample_kspace(acqData, 2.0, "poisson", calsize=25,profiles=false);
 
 
+# CS reconstruction using Wavelets
+params = Dict{Symbol,Any}()
+params[:reco] = "direct"
+params[:reconSize] = (nx,ny)
+
+@time img_u = reconstruction(acqData,params)
+@info "relative error: $(norm(img-img_u)/norm(img))"
+
 ## Load the Dictionary, build the sparsifying transform and perform the reconstruction
 
 # load the dictionary
-D = ComplexF64.(readdlm("data/brainDict98.tsv"))
+D = ComplexF32.(readdlm(joinpath(@__DIR__,"data/brainDict98.tsv")))
+
+begin
+D2 = reshape(D,6,6,:)# check patch
+f=Figure()
+for i in 1:6, j in 1:6
+  ax=Axis(f[i,j],aspect=1)
+  heatmap!(ax,abs.(D2[:,:,(i-1)*6+j]))
+  hidedecorations!(ax)
+end
+f
+end
 
 # some parameters
 px, py = (6,6)  # patch size
@@ -66,39 +87,62 @@ params = Dict{Symbol,Any}()
 params[:reco] = "standard"
 params[:reconSize] = (nx,ny)
 params[:iterations] = 50
-params[:reg] = L1Regularization(2.e-2)
-params[:sparseTrafo] = dictOp(D,(nx,ny),(px,py),2.e-2)
-params[:ρ] = 0.1
+params[:reg] = L1Regularization(0.0)
+params[:sparseTrafo] = dictOp(D,(nx,ny),(px,py),2e-2)
+params[:rho] = 0.1
 params[:solver] = ADMM
 params[:absTol] = 1.e-4
 params[:relTol] = 1.e-2
+params[:normalizeReg] = MeasurementBasedNormalization()
+params[:kwargWarning] = true
 
-img_d = reconstruction(acqData,params)
+@time img_d = reconstruction(acqData,params)
 @info "relative error: $(norm(img-img_d)/norm(img))"
 
 
+# use CairoMakie for interactive display
+begin
+  colormap=:grays
+  f = Figure(size=(1000,300))
+  ax = Axis(f[1,1],title="Original")
+  heatmap!(ax,rotr90(abs.(img[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,2],title="Undersampled")
+  heatmap!(ax,rotr90(abs.(img_u[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,3],title="Wavelet")
+  heatmap!(ax,rotr90(abs.(img_w[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,4],title="Dictionary")
+  heatmap!(ax,rotr90(abs.(img_d[:,:,1,1,1]));colormap)
+  [hidedecorations!(f.content[ax]) for ax in eachindex(f.content)]
+  f
+end
 # For comparison, let us perform the same reconstruction as above but with a Wavelet transform
 
-delete!(params, :sparseTrafo)
-params[:sparseTrafoName] = "Wavelet"
+params = Dict{Symbol,Any}()
+params[:reco] = "standard"
+params[:reconSize] = (nx,ny)
+params[:iterations] = 50
+params[:reg] = L1Regularization(2.e-2)
+params[:sparseTrafo] = "Wavelet"
 
 img_w = reconstruction(acqData,params)
 @info "relative error: $(norm(img-img_w)/norm(img))"
 
 
-# use PyPlot for interactive display
-figure(34)
-clf()
-subplot(2,2,1)
-title("Original")
-imshow(abs.(img[:,:,1,1,1]))
-subplot(2,2,3)
-title("Wavelet")
-imshow(abs.(img_w[:,:,1,1,1]))
-subplot(2,2,4)
-title("Dictionary")
-imshow(abs.(img_d[:,:,1,1,1]))
-
+# use CairoMakie for interactive display
+begin
+  colormap=:grays
+  f = Figure(size=(1000,300))
+  ax = Axis(f[1,1],title="Original")
+  heatmap!(ax,rotr90(abs.(img[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,2],title="Undersampled")
+  heatmap!(ax,rotr90(abs.(img_u[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,3],title="Wavelet")
+  heatmap!(ax,rotr90(abs.(img_w[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,4],title="Dictionary")
+  heatmap!(ax,rotr90(abs.(img_d[:,:,1,1,1]));colormap)
+  [hidedecorations!(f.content[ax]) for ax in eachindex(f.content)]
+  f
+end
 
 # export images
 filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/brainOrig.png")

docs/src/examples/exampleFieldmap.jl

@@ -1,4 +1,6 @@
-using MRIReco, MRIReco.RegularizedLeastSquares
+using MRIReco, MRIReco.RegularizedLeastSquares, MRISimulation
+using ImageUtils: shepp_logan
+include(joinpath(@__DIR__,"exampleUtils.jl"))
 
 ##### simple example ####
 

docs/src/examples/exampleIO.jl

@@ -1,4 +1,5 @@
-using HTTP, PyPlot, MRIReco
+using HTTP, CairoMakie, MRIReco
+include(joinpath(@__DIR__,"exampleUtils.jl"))
 
 if !isdir("brukerfileCart")
   HTTP.open("GET", "http://media.tuhh.de/ibi/mrireco/brukerfileCart.zip") do http

docs/src/examples/exampleRadial.jl

@@ -1,5 +1,6 @@
-using PyPlot, MRIReco 
-
+using CairoMakie, MRIReco, MRISimulation
+using ImageUtils: shepp_logan
+include(joinpath(@__DIR__,"exampleUtils.jl"))
 ##### simple example ####
 
 N = 256
@@ -21,15 +22,17 @@ params[:reco] = "direct"
 params[:reconSize] = (N,N)
 Ireco = reconstruction(acqData, params)
 
-# use PyPlot for interactive display
-figure(1)
-clf()
-subplot(1,2,1)
-title("Phantom")
-imshow(abs.(I))
-subplot(1,2,2)
-title("Reconstruction")
-imshow(abs.(Ireco[:,:,1,1,1]))
+# use CairoMakie for  display
+begin
+  colormap=:grays
+  f = Figure(size=(500,250))
+  ax = Axis(f[1,1],title="Phantom")
+  heatmap!(ax,rotr90(abs.(I));colormap)
+  ax = Axis(f[1,2],title="Reconstruction")
+  heatmap!(ax,rotr90(abs.(Ireco[:,:,1,1,1])))
+  [hidedecorations!(f.content[ax]) for ax in eachindex(f.content)]
+  f
+end
 
 # export images
 filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/simpleReco.png")

docs/src/examples/exampleSENSE.jl

@@ -1,5 +1,7 @@
-using PyPlot, MRIReco, MRIReco.RegularizedLeastSquares
+using CairoMakie, MRIReco, MRIReco.RegularizedLeastSquares
+using MRICoilSensitivities, MRISimulation
 
+include(joinpath(@__DIR__,"exampleUtils.jl"))
 ##### simple example ####
 
 N = 256
@@ -26,22 +28,31 @@ acqData = simulation(I, params)
 params = Dict{Symbol, Any}()
 params[:reco] = "multiCoil"
 params[:reconSize] = (N,N)
-params[:reg] = L2Regularization(1.e-3)
-params[:iterations] = 40
+params[:reg] = L2Regularization(1.e-1)
+params[:iterations] = 1
 params[:solver] = CGNR
 params[:senseMaps] = coilsens
 params[:toeplitz] = false
+params[:normalizeReg] = MeasurementBasedNormalization()
+
+@time Ireco_u = reconstruction(acqData, params)
+
+params[:iterations] = 40
 @time Ireco = reconstruction(acqData, params)
 
-# use PyPlot for interactive display
-figure(1)
-clf()
-subplot(1,2,1)
-title("Phantom")
-imshow(abs.(I))
-subplot(1,2,2)
-title("Reconstruction")
-imshow(abs.(Ireco[:,:,1,1,1]))
+# use CairoMakie for  display
+begin
+  colormap=:grays
+  f = Figure(size=(750,250))
+  ax = Axis(f[1,1],title="Phantom")
+  heatmap!(ax,rotr90(abs.(I));colormap)
+  ax = Axis(f[1,2],title="Reconstruction iteration 1")
+  heatmap!(ax,rotr90(abs.(Ireco_u[:,:,1,1,1]));colormap)
+  ax = Axis(f[1,3],title="Reconstruction iteration 40")
+  heatmap!(ax,rotr90(abs.(Ireco[:,:,1,1,1]));colormap)
+  [hidedecorations!(f.content[ax]) for ax in eachindex(f.content)]
+  f
+end
 
 # export images
 filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/senseReco.png")