From a8b71388e1c789a9d4c6f7f3f7d5db19b5e7be10 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 12 Jul 2024 18:17:02 +0200 Subject: [PATCH 1/2] Update MPIReco API and move to CairoMakie for plotting --- examples/reco1D.jl | 14 +++++----- examples/reco2D.jl | 14 +++++----- examples/reco3D.jl | 17 +++++++----- examples/visualization.jl | 56 +++++++++++++++++---------------------- 4 files changed, 50 insertions(+), 51 deletions(-) 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/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 From 84c4c258662d6c354ec8521c7c799255abdd386a Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 12 Jul 2024 18:19:36 +0200 Subject: [PATCH 2/2] Init example for GPU reco (no good pictures yet) --- examples/reco3D_gpu.jl | 46 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 examples/reco3D_gpu.jl 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