diff --git a/examples/reco1D.jl b/examples/reco1D.jl
index 7beaccf..416f692 100644
--- a/examples/reco1D.jl
+++ b/examples/reco1D.jl
@@ -1,16 +1,18 @@
-using PyPlot, MPIReco, OpenMPIData
+using CairoMakie, MPIReco, OpenMPIData, MPIReco.RegularizedLeastSquares
 
 include("visualization.jl")
 
 filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","1.mdf")
 #filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","4.mdf")
+calib = MPIFile(filenameCalib)
 
 for (i,phantom) in enumerate(["shapePhantom", "resolutionPhantom", "concentrationPhantom"])
 
   filenameMeas = joinpath(OpenMPIData.basedir(),"data","measurements",phantom,"1.mdf")
+  meas = MPIFile(filenameMeas)
 
   # reconstruct data
-  c = reconstruction(filenameCalib, filenameMeas, iterations=10, lambd=0.001,
+  c = reconstruct("SinglePatch", meas; sf = calib, iterations=10, reg= [L2Regularization(0.001)], solver = Kaczmarz,
                            minFreq=80e3, SNRThresh=1.25,  recChannels=1:3,
                    bgCorrectionInternal=true, spectralLeakageCorrection=false)[:,:,10:28,10:28,:]
 
@@ -20,18 +22,18 @@ for (i,phantom) in enumerate(["shapePhantom", "resolutionPhantom", "concentratio
   # visualization  
   if phantom =="shapePhantom"
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction1D.png")
-    showMIPs(c[1,:,:,:,1],filename=filenameImage, fignum=i)
+    showMIPs(c[1,:,:,:,1],filename=filenameImage)
   elseif phantom =="resolutionPhantom"
     slice=[div(s[1]+1,2),div(s[2]+1,2),div(s[3]+1,2)]
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction1D.png")
-    showSlices(c[1,:,:,:,1],slice,filename=filenameImage, fignum=i)
+    showSlices(c[1,:,:,:,1],slice,filename=filenameImage)
   elseif phantom =="concentrationPhantom"
     slice1=[div(s[1],3)+1,div(s[2],3)+1,div(s[3],3)+1]
     slice2=[2*div(s[1],3)+1,2*div(s[2],3)+1,div(s[3],3)+1]
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction1D_1.png")
-    showSlices(c[1,:,:,:,1],slice1,filename=filenameImage,fignum=i)
+    showSlices(c[1,:,:,:,1],slice1,filename=filenameImage)
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction1D_2.png")
-    showSlices(c[1,:,:,:,1],slice2,filename=filenameImage,fignum=i+1)
+    showSlices(c[1,:,:,:,1],slice2,filename=filenameImage)
   end
 end
 
diff --git a/examples/reco2D.jl b/examples/reco2D.jl
index 35f9f7a..d5c4bfb 100644
--- a/examples/reco2D.jl
+++ b/examples/reco2D.jl
@@ -1,16 +1,18 @@
-using PyPlot, MPIReco, OpenMPIData
+using CairoMakie, MPIReco, OpenMPIData
 
 include("visualization.jl")
 
 filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","2.mdf")
 #filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","5.mdf")
+calib = MPIFile(filenameCalib)
 
 for (i,phantom) in enumerate(["shapePhantom", "resolutionPhantom", "concentrationPhantom"])
 
   filenameMeas = joinpath(OpenMPIData.basedir(),"data","measurements",phantom,"2.mdf")
+  meas = MPIFile(filenameMeas)
 
   # reconstruct data
-  c = reconstruction(filenameCalib, filenameMeas, iterations=10, lambd=0.01,
+  c = reconstruct("SinglePatch", meas; sf = calib, iterations=10, reg= [L2Regularization(0.01)], solver = Kaczmarz,
                     minFreq=80e3, SNRThresh=1.5, recChannels=1:3,
 		    bgCorrectionInternal=true, spectralLeakageCorrection=false)[:,:,:,10:28,:]
 
@@ -20,17 +22,17 @@ for (i,phantom) in enumerate(["shapePhantom", "resolutionPhantom", "concentratio
   # visualization  
   if phantom =="shapePhantom"
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction2D.png")
-    showMIPs(c[1,:,:,:,1],filename=filenameImage, fignum=i)
+    showMIPs(c[1,:,:,:,1],filename=filenameImage)
   elseif phantom =="resolutionPhantom"
     slice=[div(s[1]+1,2),div(s[2]+1,2),div(s[3]+1,2)]
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction2D.png")
-    showSlices(c[1,:,:,:,1],slice,filename=filenameImage, fignum=i)
+    showSlices(c[1,:,:,:,1],slice,filename=filenameImage)
   elseif phantom =="concentrationPhantom"
     slice1=[div(s[1],3)+1,div(s[2],3)+1,div(s[3],3)+1]
     slice2=[2*div(s[1],3)+1,2*div(s[2],3)+1,div(s[3],3)+1]
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction2D_1.png")
-    showSlices(c[1,:,:,:,1],slice1,filename=filenameImage,fignum=i)
+    showSlices(c[1,:,:,:,1],slice1,filename=filenameImage)
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction2D_2.png")
-    showSlices(c[1,:,:,:,1],slice2,filename=filenameImage,fignum=i+1)
+    showSlices(c[1,:,:,:,1],slice2,filename=filenameImage)
   end
 end
diff --git a/examples/reco3D.jl b/examples/reco3D.jl
index 66b986d..531ebfa 100644
--- a/examples/reco3D.jl
+++ b/examples/reco3D.jl
@@ -1,17 +1,20 @@
-using PyPlot, MPIReco, OpenMPIData
+using CairoMakie, MPIReco, OpenMPIData, MPIReco.RegularizedLeastSquares
 
 include("visualization.jl")
 
 filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","3.mdf")
 #filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","6.mdf")#High Resolution
+calib = MPIFile(filenameCalib)
+
 
 for (i,phantom) in enumerate(["shapePhantom", "resolutionPhantom", "concentrationPhantom"])
 
   filenameMeas = joinpath(OpenMPIData.basedir(),"data","measurements",phantom,"3.mdf")
+  meas = MPIFile(filenameMeas)
 
   # reconstruct data
-  c = reconstruction(filenameCalib, filenameMeas, iterations=3, lambd=0.001, bgCorrectionInternal=false,
-                   minFreq=80e3, SNRThresh=2.0, recChannels=1:3, frames=1:1000, numAverages=1000)
+  c = reconstruct("SinglePatch", meas; sf = calib, iterations=3, reg= [L2Regularization(0.001)], bgCorrectionInternal=false,
+                   minFreq=80e3, SNRThresh=2.0, recChannels=1:3, frames=1:1000, numAverages=1000, solver = Kaczmarz)
 
   mkpath( joinpath(OpenMPIData.basedir(),"data/reconstructions/$(phantom)"))
   s = size(c)[2:4]
@@ -19,18 +22,18 @@ for (i,phantom) in enumerate(["shapePhantom", "resolutionPhantom", "concentratio
   # visualization
   if phantom =="shapePhantom"
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D.png")
-    showMIPs(c[1,:,:,:,1],filename=filenameImage,fignum=i)
+    showMIPs(c[1,:,:,:,1],filename=filenameImage)
   elseif phantom =="resolutionPhantom"
     slice=[div(s[1]+1,2),div(s[2]+1,2),div(s[3]+1,2)]
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D.png")
-    showSlices(c[1,:,:,:,1],slice,filename=filenameImage, fignum=i)
+    showSlices(c[1,:,:,:,1],slice,filename=filenameImage)
   elseif phantom =="concentrationPhantom"
     slice1 = [div(s[1],3)+1,div(s[2],3)+1,div(s[3],3)+1]
     slice2 = [2*div(s[1],3)+1,2*div(s[2],3)+1,2*div(s[3],3)+1]
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D_1.png")
-    showSlices(c[1,:,:,:,1],slice1,filename=filenameImage,fignum=i)
+    showSlices(c[1,:,:,:,1],slice1,filename=filenameImage)
     filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D_2.png")
-    showSlices(c[1,:,:,:,1],slice2,filename=filenameImage,fignum=i+1)
+    showSlices(c[1,:,:,:,1],slice2,filename=filenameImage)
   end
 
 end
diff --git a/examples/reco3D_gpu.jl b/examples/reco3D_gpu.jl
new file mode 100644
index 0000000..aa2d6ab
--- /dev/null
+++ b/examples/reco3D_gpu.jl
@@ -0,0 +1,46 @@
+using CairoMakie, MPIReco, OpenMPIData, MPIReco.RegularizedLeastSquares
+
+# Add GPU package
+using CUDA # or AMDGPU, Metal, ...
+arrayType = CuArray # or ROCArray, ....
+
+
+include("visualization.jl")
+
+filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","3.mdf")
+#filenameCalib = joinpath(OpenMPIData.basedir(),"data","calibrations","6.mdf")#High Resolution
+calib = MPIFile(filenameCalib)
+
+
+for (i,phantom) in enumerate(["shapePhantom"]) # , "resolutionPhantom", "concentrationPhantom"])
+
+  filenameMeas = joinpath(OpenMPIData.basedir(),"data","measurements",phantom,"3.mdf")
+  meas = MPIFile(filenameMeas)
+
+  # reconstruct data
+  c = reconstruct("SinglePatch", meas; sf = calib, iterations=300, reg= [L2Regularization(0.1)], rho = 0.95, bgCorrectionInternal=false,
+                   minFreq=80e3, SNRThresh=2.0, recChannels=1:3, frames=1:1000, numAverages=1000,
+                   solver = CGNR, # Kaczmarz is not well suited for GPU recos 
+                   arrayType = arrayType) # Specify arraytype for GPU reconstruction
+
+  mkpath( joinpath(OpenMPIData.basedir(),"data/reconstructions/$(phantom)"))
+  s = size(c)[2:4]
+  
+  # visualization
+  if phantom =="shapePhantom"
+    filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D_$(parentmodule(arrayType)).png")
+    showMIPs(c[1,:,:,:,1],filename=filenameImage)
+  elseif phantom =="resolutionPhantom"
+    slice=[div(s[1]+1,2),div(s[2]+1,2),div(s[3]+1,2)]
+    filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D_$(parentmodule(arrayType)).png")
+    showSlices(c[1,:,:,:,1],slice,filename=filenameImage)
+  elseif phantom =="concentrationPhantom"
+    slice1 = [div(s[1],3)+1,div(s[2],3)+1,div(s[3],3)+1]
+    slice2 = [2*div(s[1],3)+1,2*div(s[2],3)+1,2*div(s[3],3)+1]
+    filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D_1_$(parentmodule(arrayType)).png")
+    showSlices(c[1,:,:,:,1],slice1,filename=filenameImage)
+    filenameImage = joinpath(OpenMPIData.basedir(),"data","reconstructions","$phantom","reconstruction3D_2_$(parentmodule(arrayType)).png")
+    showSlices(c[1,:,:,:,1],slice2,filename=filenameImage)
+  end
+
+end
diff --git a/examples/visualization.jl b/examples/visualization.jl
index f8540e0..dfd77f6 100644
--- a/examples/visualization.jl
+++ b/examples/visualization.jl
@@ -1,33 +1,29 @@
 export showMIPs, showSlices
 
-PyPlot.rc("font", family="serif")
-PyPlot.rc("text", usetex=true)
+Makie.update_theme!(fonts = (regular = Makie.texfont(), bold = Makie.texfont(:bold), italic = Makie.texfont(:italic)))
 
-function showMIPs(c; fignum=1, filename=nothing)
+function showMIPs(c;filename=nothing)
   cAbs = abs.(arraydata(c))
   cxy = maximum(cAbs,dims=3)[:,:,1]
   cxz = maximum(cAbs,dims=2)[:,1,:]
   cyz = maximum(cAbs,dims=1)[1,:,:]
 
-  figure(fignum, figsize=(4,4))
-  
-  subplot(2,2,1)
-  imshow(cxz', interpolation="nearest")
-  title("MIP xz")
-  subplot(2,2,2)
-  imshow(cyz', interpolation="nearest")
-  title("MIP yz")
-  subplot(2,2,4)
-  imshow(cxy, interpolation="nearest")
-  title("MIP xy")
-
-  subplots_adjust(wspace=0.18,hspace=0.3,left=0.06,bottom=0.06,right=1.0,top=0.93)
+  fig = Figure()
+  g = GridLayout(fig[1,1])
+  ax = CairoMakie.Axis(g[1,1], title = "MIP xz", width = 150, height = 150)
+  heatmap!(ax, cxz)
+  ax = CairoMakie.Axis(g[1,2], title = "MIP yz", width = 150, height = 150 )
+  heatmap!(ax, cyz)
+  ax = CairoMakie.Axis(g[2,2], title = "MIP xy", width = 150, height = 150)
+  heatmap!(ax, cxy)
 
+  resize_to_layout!(fig)
   if filename != nothing
-    savefig(filename)
+    save(filename, fig)
   end
 end
-function showSlices(c, slice; fignum=1, filename=nothing)
+
+function showSlices(c, slice;filename=nothing)
   cAbs = abs.(arraydata(c))
   x=slice[1]
   y=slice[2]
@@ -37,21 +33,17 @@ function showSlices(c, slice; fignum=1, filename=nothing)
   cxz = cAbs[:,y,:]
   cyz = cAbs[x,:,:]
 
-  figure(fignum, figsize=(4,4))
-  
-  subplot(2,2,1)
-  imshow(cxz', interpolation="nearest")
-  title("Slice at y=$y xz")
-  subplot(2,2,2)
-  imshow(cyz', interpolation="nearest")
-  title("Slice at x=$x yz")
-  subplot(2,2,4)
-  imshow(cxy, interpolation="nearest")
-  title("Slice at z=$z xy")
-
-  subplots_adjust(wspace=0.18,hspace=0.3,left=0.06,bottom=0.06,right=1.0,top=0.93)
+  fig = Figure()
+  g = GridLayout(fig[1,1])
+  ax = CairoMakie.Axis(g[1,1], title = "Slice at y=$y xz", width = 150, height = 150)
+  heatmap!(ax, cxz)
+  ax = CairoMakie.Axis(g[1, 2], title = "Slice at x=$x yz", width = 150, height = 150)
+  heatmap!(ax, cyz)
+  ax = CairoMakie.Axis(g[2,2], title = "Slice at z=$z xy", width = 150, height = 150)
+  heatmap!(ax, cxy)
 
+  resize_to_layout!(fig)
   if filename != nothing
-    savefig(filename)
+    save(filename, fig)
   end
 end