GitHub

Multi-Patch Reconstruction

For multi-patch reconstruction the method proposed by Szwargulski et al. is implemented in MPIReco. It is generalized however as described in Boberg et al..

We first discuss the measurements for the multi-patch case. On modern MPI scanners the BrukerFile or MDFFile can be used as is. However, the data that we use in our unit tests consists of several single-patch measurements. We therefore combine them manually into an MultiMPIFile:

files = ["1.mdf", "2.mdf", "3.mdf", "4.mdf"]
b = MultiMPIFile(joinpath.(datadir, "measurements", "20211226_203916_MultiPatch", files))
Multi MPI File:
1. MDFFile
	Study: "FocusField", Dates.DateTime("2014-11-13T16:20:01.217")
	Experiment: "Ph5DotsDF10Overlap0LL (E99)", Dates.DateTime("2014-11-27T14:40:01.280")
2. MDFFile
	Study: "FocusField", Dates.DateTime("2014-11-13T16:20:01.217")
	Experiment: "Ph5DotsDF10Overlap0LR  (E101)", Dates.DateTime("2014-11-27T14:41:47.871")
3. MDFFile
	Study: "FocusField", Dates.DateTime("2014-11-13T16:20:01.217")
	Experiment: "Ph5DotsDF10Overlap0RL (E100)", Dates.DateTime("2014-11-27T14:40:55.227")
4. MDFFile
	Study: "FocusField", Dates.DateTime("2014-11-13T16:20:01.217")
	Experiment: "Ph5DotsDF10Overlap0RR (E98)", Dates.DateTime("2014-11-27T14:39:13.121")

And now b can be used as if it was a multi-patch file.

Next we take a look at the system matrix. The most simple approach is to use a single system matrix that was measured at the center. This can be done using

bSF = MultiMPIFile([joinpath(datadir, "calibrations", "12.mdf")])
Multi MPI File:
1. MDFFile
	Study: "FocusField", Dates.DateTime("2014-11-13T16:20:01.217")
	Experiment: "SF_DF10_Middle (E39)", Dates.DateTime("2014-11-19T19:31:31.527")

Afterwards we can perform a reconstruction with

c1 = reconstruct("MultiPatch", b; sf = bSF, SNRThresh=5, frames=1:acqNumFrames(b), minFreq=80e3,
                   recChannels=1:rxNumChannels(b), iterations=1, spectralLeakageCorrection=false, tfCorrection = false);
[ Info: Loading Multi Patch operator

The parameters are essentially the same as in the previous reconstructions, we only change the algorithm and the MPIFiles.

It is also possible to use multiple system matrices, which is currently the best way to take field imperfection into account. Our test data has four patches and we therefore can use

sf_files =  ["8.mdf", "9.mdf", "10.mdf", "11.mdf"]
bSFs = MultiMPIFile(joinpath.(datadir, "calibrations", sf_files))

c2 = reconstruct("MultiPatch", b; sf = bSFs, SNRThresh=5, frames=1:acqNumFrames(b), minFreq=80e3,
                   recChannels=1:rxNumChannels(b), iterations=1, spectralLeakageCorrection=false, tfCorrection = false);
[ Info: Loading Multi Patch operator

Now we want somewhat more flexibility and

  • define a mapping between the system matrix and the patches, here we allow to use the same system matrix for multiple patches
  • make it possible to change the FFP position. Usually the value stored in the file is not 100% correct due to field imperfections.

To achieve this, we swap out parts of the MultiPatch blueprint with a parameter which allows us to explicitly set all these parameters

We perform our own frequency filtering.

freq = filterFrequencies(bSFs, SNRThresh=5, minFreq=80e3);

And load four different system matrices for each patch.

S = [getSF(SF,freq,nothing,"Kaczmarz", bgcorrection=false, tfCorrection = false)[1] for SF in bSFs]
4-element Vector{LinearAlgebra.Transpose{ComplexF32, Matrix{ComplexF32}}}:
 [-17.215738f0 + 3.6095254f0im -2.1378102f0 - 11.237434f0im … 6.6960473f0 + 1.5645299f0im 11.383284f0 - 6.8314924f0im; -1.6692096f0 + 15.541314f0im -8.507132f0 + 0.33559847f0im … 6.9464974f0 - 9.465756f0im 4.9152975f0 - 1.861226f0im; … ; -0.12820308f0 + 0.44176456f0im -1.528524f0 - 1.4499999f0im … -0.6129609f0 + 4.7460117f0im -3.0676806f0 + 0.080153696f0im; 0.99223995f0 + 2.846132f0im -0.515808f0 + 3.7751706f0im … -0.029549912f0 + 0.12914856f0im 1.7864641f0 - 1.4239423f0im]
 [17.633753f0 + 2.541937f0im 13.33199f0 + 17.208555f0im … 11.224207f0 - 6.787769f0im 18.090244f0 - 17.947556f0im; 121.56271f0 - 192.6314f0im 121.15709f0 - 178.92175f0im … 128.29367f0 - 216.8749f0im 127.61627f0 - 213.36526f0im; … ; 0.93170273f0 - 0.6052504f0im 0.8655644f0 + 1.2170328f0im … 0.20082603f0 - 0.053652916f0im -4.9453936f0 + 0.2074228f0im; 2.3107705f0 - 1.3597172f0im -0.11082725f0 + 0.57934487f0im … 2.7711606f0 - 1.3732481f0im 2.9816546f0 - 1.540642f0im]
 [15.31848f0 + 9.865552f0im 18.600489f0 - 7.4320703f0im … 9.798755f0 - 4.6202745f0im 6.7259164f0 + 8.711154f0im; 120.42122f0 - 211.66068f0im 116.374855f0 - 209.65962f0im … 102.539505f0 - 239.4076f0im 100.766525f0 - 232.4441f0im; … ; -2.5822783f0 - 1.2518568f0im 1.2427015f0 - 0.7835897f0im … -0.40484527f0 + 2.4156482f0im 1.8506767f0 - 0.96858203f0im; 1.9769307f0 - 1.6546817f0im -1.1609071f0 - 1.4179816f0im … 2.0310106f0 - 2.438818f0im 1.2885772f0 - 0.4624802f0im]
 [-3.1459377f0 + 10.193573f0im -9.436847f0 - 4.1981077f0im … 12.40924f0 - 9.423947f0im 6.1208096f0 - 10.054392f0im; 161.08578f0 - 201.87471f0im 160.50711f0 - 196.88f0im … 149.59499f0 - 213.31972f0im 146.94023f0 - 213.37024f0im; … ; 2.7833283f0 - 1.8007818f0im -0.81187344f0 - 2.141356f0im … -0.8005118f0 - 2.475382f0im 1.080762f0 - 2.2345705f0im; -0.11865397f0 - 1.9258574f0im -2.1296504f0 - 0.80380076f0im … -1.0891958f0 + 0.36377603f0im -0.34626928f0 - 1.5952204f0im]

Afterwards we can describe our patch mapping and positions:

mapping = [1,2,3,4]
SFGridCenter = zeros(3,4)
FFPos = zeros(3,4)
FFPos[:,1] = [-0.008, 0.008, 0.0]
FFPos[:,2] = [-0.008, -0.008, 0.0]
FFPos[:,3] = [0.008, 0.008, 0.0]
FFPos[:,4] = [0.008, -0.008, 0.0];

Lastly, we wrap everything in a parameter structs:

opParams = ExplicitMultiPatchParameter(;tfCorrection = false, systemMatrices = S, SFGridCenter = SFGridCenter, mapping = mapping)
ffPos = CustomFocusFieldPositions(FFPos);

and perform our reconstruction:

c3 = reconstruct("MultiPatch", b; sf = bSFs, opParams = opParams, ffPos = ffPos, ffPosSF = ffPos,
                  SNRThresh=5, frames=1:acqNumFrames(b), minFreq=80e3,
                  recChannels=1:rxNumChannels(b), iterations=1, spectralLeakageCorrection=false, tfCorrection = false);
[ Info: Loading Multi Patch operator

We can again visualize our different multi-patch reconstructions:

fig = Figure();
hidedecorations!(heatmap(fig[1, 1], c1[1, :, :, 1, 1].data.data, axis = (title = "Single SM",)).axis)
hidedecorations!(heatmap(fig[1, 2], c2[1, :, :, 1, 1].data.data, axis = (title = "Multiple SM",)).axis)
hidedecorations!(heatmap(fig[1, 3], c3[1, :, :, 1, 1].data.data, axis = (title = "Explicit SM",)).axis)
fig
Example block output

This page was generated using Literate.jl.