diff --git a/docs/Project.toml b/docs/Project.toml
index 82dd8046..5685e8b7 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/docs/src/assets/bruker.png b/docs/src/assets/bruker.png
index 788c7f9f..16257c5f 100644
Binary files a/docs/src/assets/bruker.png and b/docs/src/assets/bruker.png differ
diff --git a/docs/src/assets/fieldmap.png b/docs/src/assets/fieldmap.png
index 90109555..0f020dad 100644
Binary files a/docs/src/assets/fieldmap.png and b/docs/src/assets/fieldmap.png differ
diff --git a/docs/src/assets/fieldmapReco1.png b/docs/src/assets/fieldmapReco1.png
index 8b4e1a3c..31ea59a0 100644
Binary files a/docs/src/assets/fieldmapReco1.png and b/docs/src/assets/fieldmapReco1.png differ
diff --git a/docs/src/assets/fieldmapReco2.png b/docs/src/assets/fieldmapReco2.png
index 29d1de91..16c49950 100644
Binary files a/docs/src/assets/fieldmapReco2.png and b/docs/src/assets/fieldmapReco2.png differ
diff --git a/docs/src/assets/kneeCG.png b/docs/src/assets/kneeCG.png
index ab4e7cac..812b1bbd 100644
Binary files a/docs/src/assets/kneeCG.png and b/docs/src/assets/kneeCG.png differ
diff --git a/docs/src/assets/kneeOrig.png b/docs/src/assets/kneeOrig.png
index c5fdcaf6..3da5a1fc 100644
Binary files a/docs/src/assets/kneeOrig.png and b/docs/src/assets/kneeOrig.png differ
diff --git a/docs/src/assets/kneeTV.png b/docs/src/assets/kneeTV.png
index c0258572..029a9cf9 100644
Binary files a/docs/src/assets/kneeTV.png and b/docs/src/assets/kneeTV.png differ
diff --git a/docs/src/assets/mask.png b/docs/src/assets/mask.png
index db8bc058..201e1324 100644
Binary files a/docs/src/assets/mask.png and b/docs/src/assets/mask.png differ
diff --git a/docs/src/assets/phantom.png b/docs/src/assets/phantom.png
index 7beccaa6..43ce58de 100644
Binary files a/docs/src/assets/phantom.png and b/docs/src/assets/phantom.png differ
diff --git a/docs/src/assets/senseReco.png b/docs/src/assets/senseReco.png
index c8b1b237..70ad77ef 100644
Binary files a/docs/src/assets/senseReco.png and b/docs/src/assets/senseReco.png differ
diff --git a/docs/src/assets/simpleReco.png b/docs/src/assets/simpleReco.png
index f941045e..7204945e 100644
Binary files a/docs/src/assets/simpleReco.png and b/docs/src/assets/simpleReco.png differ
diff --git a/docs/src/assets/trajectories.png b/docs/src/assets/trajectories.png
index 7c3f4170..aacc8ef7 100644
Binary files a/docs/src/assets/trajectories.png and b/docs/src/assets/trajectories.png differ
diff --git a/docs/src/examples/exampleCS.jl b/docs/src/examples/exampleCS.jl
index 1cf81411..471f63b7 100644
--- a/docs/src/examples/exampleCS.jl
+++ b/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]) )
\ No newline at end of file
diff --git a/docs/src/examples/exampleCustom.jl b/docs/src/examples/exampleCustom.jl
index 5978f15e..1bc57b48 100644
--- a/docs/src/examples/exampleCustom.jl
+++ b/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,41 @@ 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")))
 
 # some parameters
 px, py = (6,6)  # patch size
@@ -67,38 +77,61 @@ 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[: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")
diff --git a/docs/src/examples/exampleFieldmap.jl b/docs/src/examples/exampleFieldmap.jl
index 44f3e1c3..6681b10a 100644
--- a/docs/src/examples/exampleFieldmap.jl
+++ b/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 ####
 
diff --git a/docs/src/examples/exampleIO.jl b/docs/src/examples/exampleIO.jl
index 00a045eb..e8531b25 100644
--- a/docs/src/examples/exampleIO.jl
+++ b/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
diff --git a/docs/src/examples/exampleRadial.jl b/docs/src/examples/exampleRadial.jl
index 7de57fa4..4def58bf 100644
--- a/docs/src/examples/exampleRadial.jl
+++ b/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")
diff --git a/docs/src/examples/exampleSENSE.jl b/docs/src/examples/exampleSENSE.jl
index 9f52fec7..e93a70e0 100644
--- a/docs/src/examples/exampleSENSE.jl
+++ b/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")
diff --git a/docs/src/examples/exampleTrajectories.jl b/docs/src/examples/exampleTrajectories.jl
index a3a5de3f..a10b655a 100644
--- a/docs/src/examples/exampleTrajectories.jl
+++ b/docs/src/examples/exampleTrajectories.jl
@@ -1,45 +1,41 @@
-using PyPlot
+using CairoMakie
 using MRIReco
 
 #### trajectories ####
 
 trajectories = ["Spiral", "Radial", "Cartesian", "EPI", "SpiralVarDens"]
 
-figure(1)
-clf()
+f = Figure()
 
-tr = trajectory("Spiral", 1, 600, windings=10)
-subplot(2,2,1)
+T=Float32
+
+tr = trajectory(Float32,"Spiral", 1, 600, windings=10)
 nodes = kspaceNodes(tr)
-plot(nodes[1,:], nodes[2,:],"b.",lw=2)
-title("Spiral")
+ax=Axis(f[1,1],title="Spiral")
+plot!(ax,nodes[1,:], nodes[2,:])
+
 
-tr = trajectory("Cartesian", 13, 50, EPI_factor=1)
-subplot(2,2,2)
+tr = trajectory(T,"Cartesian", 13, 50, EPI_factor=1)
 nodes = kspaceNodes(tr)
-plot(nodes[1,:], nodes[2,:],"b.",lw=2)
-title("Cartesian")
+ax=Axis(f[1,2],title="Cartesian")
+plot!(ax,nodes[1,:], nodes[2,:])
+
 
 
-tr = trajectory("Radial", 13, 50)
-subplot(2,2,3)
+tr = trajectory(T,"Radial", 13, 50)
 nodes = kspaceNodes(tr)
-plot(nodes[1,:], nodes[2,:],"b.",lw=2)
-title("Radial")
+ax=Axis(f[2,1],title="Radial")
+plot!(ax,nodes[1,:], nodes[2,:])
+
+
 
 
-tr = trajectory("SpiralVarDens", 3, 200)
-subplot(2,2,4)
+tr = trajectory(T,"SpiralVarDens", 3, 200)
 nodes = kspaceNodes(tr)
-plot(nodes[1,:], nodes[2,:],"b.",lw=2)
-title("SpiralVarDens")
+ax=Axis(f[2,2],title="SpiralVarDens")
+plot!(ax,nodes[1,:], nodes[2,:])
+f
 
-subplots_adjust(top=0.95,
-       bottom=0.05,
-       left=0.09,
-       right=0.95,
-       hspace=0.3,
-       wspace=0.2)
 
 filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/trajectories.png")
-savefig(filename, dpi=200)
+save(filename,f)
diff --git a/docs/src/examples/exampleUtils.jl b/docs/src/examples/exampleUtils.jl
new file mode 100644
index 00000000..17fad3e8
--- /dev/null
+++ b/docs/src/examples/exampleUtils.jl
@@ -0,0 +1,11 @@
+using CairoMakie
+function exportImage(filepath::AbstractString,img)
+  f = Figure(figure_padding = 0)
+  ax=Axis(f[1,1],aspect=1)
+  heatmap!(ax,rotr90(img ),colormap=:grays)
+  hidedecorations!(ax)
+  colsize!(f.layout, 1, Aspect(1,1))
+  resize_to_layout!(f)
+  f
+  save(filepath,f)
+end
\ No newline at end of file
diff --git a/docs/src/gettingStarted.md b/docs/src/gettingStarted.md
index 300ada89..2815a779 100644
--- a/docs/src/gettingStarted.md
+++ b/docs/src/gettingStarted.md
@@ -1,12 +1,10 @@
 # Getting Started
 
 Throughout the entire documentation we assume that you have loaded `MRIReco`
-as well as `PyPlot` via
+as well as `CairoMakie` via
 ```julia
-using PyPlot, MRIReco
+using CairoMakie, MRIReco
 ```
-Its important to load these packages in that order, since otherwise `PyPlot`
-will not work correctly on some systems.
 
 All examples discussed in the documentation can be found in the package source
 code in the folder
@@ -71,7 +69,7 @@ Ireco = reconstruction(acqData, params)
 The resulting image is of type `AxisArray` and has 5 dimensions. One can
 display the image object by calling
 ```julia
-imshow(abs.(Ireco[:,:,1,1,1]))
+heatmap(abs.(Ireco[:,:,1,1,1]))
 ```
 Alternatively one can store the image into a file, which will be discussed in
 the section on [Images](@ref).