Getting Started

To begin, we first need to gather some MPI data. To do this we could use data from the OpenMPIData initiative, however for these examples we will use the data used in testing MPIReco.

To access this test data, first, enter the Pkg mode in Julia (]) and execute the unit tests of MPIReco::

test MPIReco

This will download and unpack some MPI measurements and calibration MDF files and perform tests with the data. The download location will be printed at the start of the test execution. You can cancel the test execution once the data is downloaded by pressing Ctrl + C or closing the Julia process.

After the download, several MPI files will be present in the specified directory. All subsequent examples assume that you have assigned the path to this directory to a variable named datadir:

const datadir = joinpath("...","artifacts", "...") # enter path here

First Reconstruction

We will start looking at a small reconstruction script. First, we load MPIReco:

using MPIReco

Next, we open handles for the system matrix and measurement data. Both are created via MPIFile function, which can be, for instance, an MDFFile or a BrukerFile.

bSF = MPIFile(joinpath(datadir, "calibrations", "12.mdf"))
b = MPIFile(joinpath(datadir, "measurements", "20211226_203916_MultiPatch", "1.mdf"))

MDFFile
	Study: "FocusField", Dates.DateTime("2014-11-13T16:20:01.217")
	Experiment: "Ph5DotsDF10Overlap0LL (E99)", Dates.DateTime("2014-11-27T14:40:01.280")

To interact with the files, you can use the functionality of MPIFiles, which is also exported by MPIReco:

acqNumFrames(b)

We refer to the documentation for MPIFiles.jl and MDF for more information on the available data and functions.

Now we will perform a system-matrix based reconstruction using the reconstruct function. The function takes as arguments the algorithm to use, the measurement data, and various keyword arguments. The available parameters and their interpretation depend on the chosen algorithm. In this case, the algorithm also expects the system matrix to be set. We set the SNR threshold to 5, meaning that only matrix rows with an SNR above 5 will be used during reconstruction. The parameter frame decides which frame of the measured data should be reconstructed.

c = reconstruct("SinglePatch", b;
                   SNRThresh=5,
                   sf = bSF,
                   frames=1:acqNumFrames(b),
                   minFreq=80e3,
                   recChannels=1:rxNumChannels(b),
                   iterations=1,
                   spectralLeakageCorrection=true)

Float32 ImageMeta with:
  data: 5-dimensional AxisArray{Float32,5,...} with axes:
    :color, 1:1
    :x, (-27.375:1.25:11.375) mm
    :y, (-11.375:1.25:27.375) mm
    :z, (0.0:1.0:0.0) mm
    :time, (0.0:6.5280000000000005:0.0) s
And data, a 1×32×32×1×1 Array{Float32, 5}
  properties:
    tracerName: ["Resovist"]
    version: 2.0.1
    experimentDescription: Ph5DotsDF10Overlap0LL (E99)
    scannerOperator: nmrsu
    rxNumSamplingPoints: 1632
    acqOffsetField: [-0.01; 0.01; -0.0;;;]
    dfStrength: [0.01 0.01 0.0;;;]
    scannerTopology: FFP
    rxDataConversionFactor: [0.0001 0.0001 0.0001; 0.0 0.0 0.0]
    experimentNumber: 99
    dfDivider: [102; 96; 99;;]
    tracerBatch: ["0"]
    acqStartTime: 2014-11-27T14:40:01.280
    tracerInjectionTime: [DateTime("2014-11-27T14:40:01.280")]
    dfWaveform: ["sine";;]
    dim: 2
    uuid: 8eadd786-d630-4f53-ae91-c8f90721bbdc
    dfBaseFrequency: 2.5e6
    dfCycle: 0.0006528
    experimentName: Ph5DotsDF10Overlap0LL (E99)
    studyUuid: 12246d79-646a-41b8-a278-5274fc9ae7c5
    tracerVolume: [0.0]
    studyDescription: n.a.
    experimentSubject: FocusField
    experimentUuid: 787b3b7b-734a-4872-8bad-a38a64478cca
    acqNumFrames: 1
    tracerConcentration: [0.5]
    dfNumChannels: 3
    dfPhase: [1.5708 1.5708 1.5708;;;]
    scannerManufacturer: Bruker/Philips
    studyTime: 2014-11-13T16:20:01.217
    scannerFacility: Universitätsklinikum Hamburg Eppendorf
    studyName: FocusField
    scannerName: Preclinical MPI System
    tracerSolute: ["Fe"]
    size: [32, 32, 1]
    time: 2021-12-26T19:52:45.666
    experimentIsSimulation: false
    datatype: MPI
    acqNumPeriodsPerFrame: 1
    experimentIsCalibration: false
    studyNumber: 1
    rxBandwidth: 1.25e6
    tracerVendor: ["n.a."]
    rxUnit: a.u.
    acqNumAverages: 10000
    rxNumChannels: 3
    acqGradient: [-1.25 0.0 0.0; 0.0 -1.25 0.0; 0.0 0.0 2.5;;;;]

Notice that the result is not just an array of particle concentration, but contains a variety of metadata as well. This is because c is a:

typeof(c)

ImageMeta{Float32, 5, AxisArray{Float32, 5, Array{Float32, 5}, Tuple{Axis{:color, UnitRange{Int64}}, Axis{:x, StepRangeLen{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Int64}}, Axis{:y, StepRangeLen{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Int64}}, Axis{:z, StepRangeLen{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Int64}}, Axis{:time, StepRangeLen{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}, Int64}}}}, Dict{Symbol, Any}}

You can access metadata as properties, such as:

c.tracerName

1-element Vector{String}:
 "Resovist"

The result can also be treated like a normal 5-dimensional array. The first dimension is for different contrasts, followed by the three spatial dimensions and lastly the temporal dimension of the different frames.

size(c)

(1, 32, 32, 1, 1)

You can access the underlying data just like any other array in Julia. For example, you can access an XY slice as follows:

slice = c[1, :, :, 1, 1]
typeof(slice)

ImageMeta{Float32, 2, AxisArrays.AxisMatrix{Float32, Matrix{Float32}, Tuple{Axis{:x, StepRangeLen{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Int64}}, Axis{:y, StepRangeLen{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Base.TwicePrecision{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(mm,), 𝐋, nothing}}}, Int64}}}}, Dict{Symbol, Any}}

Note that the slice is still an ImageMeta array. Inside the array is an AxisArray, another array type with metadata. Not every Julia package has methods for such metadata arrays. To visualize our results with CairoMakie, we need to extract the underlying data.

sliceRaw = slice.data.data
typeof(sliceRaw)

Matrix{Float32} (alias for Array{Float32, 2})

We can now visualize our first MPI reconstruction using MPIReco:

using CairoMakie
fig = heatmap(sliceRaw)
hidedecorations!(fig.axis)
fig

This page was generated using Literate.jl.