Compare maps from quantitative sequence to the phantom

[aTrotier]

Thanks to your efforts, I am now able to generate a T1 maps with my MP2RAGE sequence :

What is the best way to compare the results to the ground truth values from the phantom2D ?

[cncastillo]

Wow! Amazing results 😄 !!

If you want a visual comparison, the easiest would be using matread like in 2D brain phantom and replicate what has been done in https://github.com/JuliaHealth/KomaMRI.jl/blob/master/KomaMRIBase/src/datatypes/Phantom.jl#L198

using KomaMRI, MAT
path = dirname(pathof(KomaMRIBase))
data = matread(path*"/datatypes/phantom/brain2D.mat")
us, ss = 1, 1 # Change the resampling accordingly
class = repeat( data["axial"][1:ss:end,1:ss:end], inner=[us, us])
T1 = (class .== 23 ) * 2569 .+ # CSF
     (class .== 46 ) * 833  .+ # GM
     (class .== 70 ) * 500  .+ # WM
     (class .== 93 ) * 350  .+ # FAT1
     (class .== 116) * 900  .+ # MUSCLE
     (class .== 139) * 569  .+ # SKIN/MUSCLE
     (class .== 162) * 0    .+ # SKULL
     (class .== 185) * 0    .+ # VESSELS
     (class .== 209) * 500  .+ # FAT2
     (class .== 232) * 2569 .+ # DURA
     (class .== 255) * 500     # MARROW
p1 = plot_image(T1')

Another way of doing a visual comparison is

plot_phantom_map(brain_phantom2D(), :T1)

If you want a quantitative comparison, you would need to interpolate T1 onto the reconstructed grid:

using Interpolations
# Defining original grid
Δx = .5e-3*ss/us # This could be anisotropic
M, N = size(class)
FOVx = (M-1)*Δx #[m]
FOVy = (N-1)*Δx #[m]
x = -FOVx/2:Δx:FOVx/2 
y = -FOVy/2:Δx:FOVy/2 
# Interpolation object (caches coefficients and such)
itp = linear_interpolation((x,y), T1; extrapolation_bc=0) # extrapolates with T1=0
# Reconstructed grid, based on https://github.com/JuliaHealth/KomaMRI.jl/discussions/311#discussioncomment-8917530
FOVx_recon = FOVy_recon = 256e-3
Nx = Ny = 64
x_recon = range(-FOVx_recon/2, FOVx_recon/2, length=Nx)
y_recon = range(-FOVy_recon/2, FOVy_recon/2, length=Ny)' 
T1_GT = itp.(x_recon, y_recon)
p2 = plot_image(T1_GT')

So right now, comparing reconstructed T1/T2 maps with the GT is not the easiest, but we can:

• Create a function match_map_to_recon(obj::Phantom, recon_params) that does this automatically based on the recon params Create a function match_map_to_recon(obj::Phantom, recon_params) that does interpolate the GT maps to the reconstructed grid #359
• Enable arbitrary color maps in plot_image. Now it is hardcoded to "Greys" in here. This also applies to plot_phantom_map Enable arbitrary color maps in plot_image and plot_phantom_map. #360

If I use Viridis and limits from 0 to 3s:

p2 = plot_image(T1_GT'; zmin=0, zmax=3000) #internally I changed the colormap to Viridis

I hope this is useful!

[aTrotier]

Great !

Results :

Not sure if I have a small voxel offset because I had to rotate and reverse the T1_GT

This is the script I used :

using JLD2,CairoMakie, KomaMRI.MRIReco, QuantitativeMRI, KomaMRI
seq = read_seq("../seq/2D_MP2RAGE_spoilX_RF117.seq")

d=load("raw_simu_3D_us=3.jld2")
raw=d["raw"]

## change the label
ETL = 64
N = (64,64)
profile=0
for echo = 0:1
  for i = 1:ETL
    profile=profile+1
    raw.profiles[profile].head.idx.contrast=echo
    raw.profiles[profile].head.idx.kspace_encode_step_1=i-1
  end
end

acq=AcquisitionData(raw)

Nx, Ny = raw.params["reconSize"][1:2]
recParams = Dict{Symbol,Any}()
recParams[:reconSize] = (Nx, Ny)
recParams[:densityWeighting] = true
rec = reconstruction(acq, recParams)

f=Figure()
ax=Axis(f[1,1],title = "TI1",aspect=1)
heatmap!(ax,abs.(rec[:,:,1,1]))
ax=Axis(f[1,2],title = "TI2",aspect=1)
heatmap!(ax,abs.(rec[:,:,1,2]))
f

## T1 map
MP2 = mp2rage_comb((rec.data[:,:,1,:,1,1]))

p=ParamsMP2RAGE(seq.DEF["TI_1"],
                seq.DEF["TI_2"],
                seq.DEF["small_TR"],
                seq.DEF["TR_MP2RAGE"],
                seq.DEF["ETL"],
                seq.DEF["alpha_1"],
                seq.DEF["alpha_2"])

T1,T1range,LUT = mp2rage_T1maps(MP2,p;T1Range=0.001:0.001:10)

## load corresponding brain_phantom2D
using KomaMRIBase.MAT

data = matread("KomaMRI.jl/KomaMRIBase/src/datatypes/phantom/brain2D.mat")
us, ss = 10, 1 # Change the resampling accordingly
class = repeat( data["axial"][1:ss:end,1:ss:end], inner=[us, us])
T1p = (class .== 23 ) * 2569 .+ # CSF
     (class .== 46 ) * 833  .+ # GM
     (class .== 70 ) * 500  .+ # WM
     (class .== 93 ) * 350  .+ # FAT1
     (class .== 116) * 900  .+ # MUSCLE
     (class .== 139) * 569  .+ # SKIN/MUSCLE
     (class .== 162) * 0    .+ # SKULL
     (class .== 185) * 0    .+ # VESSELS
     (class .== 209) * 500  .+ # FAT2
     (class .== 232) * 2569 .+ # DURA
     (class .== 255) * 500     # MARROW

using KomaMRIBase.Interpolations
# Defining original grid
sx = .5e-3*ss/us # This could be anisotropic
M, N = size(class)
FOVx = (M-1)*sx #[m]
FOVy = (N-1)*sx #[m]
x = -FOVx/2:sx:FOVx/2 
y = -FOVy/2:sx:FOVy/2 
# Interpolation object (caches coefficients and such)
itp = linear_interpolation((x,y), T1p; extrapolation_bc=0) # extrapolates with T1=0
# Reconstructed grid, based on https://github.com/JuliaHealth/KomaMRI.jl/discussions/311#discussioncomment-8917530
FOVx_recon = FOVy_recon = 256e-3
Nx = Ny = 64
x_recon = Base.range(-FOVx_recon/2, FOVx_recon/2, length=Nx)
y_recon = Base.range(-FOVy_recon/2, FOVy_recon/2, length=Ny)' 
T1_GT = itp.(x_recon, y_recon)

# rotate
T1_GT=reverse(rotl90(T1_GT),dims=1)

begin
    f=Figure()

    ax=Axis(f[1,1],title = "T1_GT",aspect=1)
    h=heatmap!(ax,T1_GT,colorrange = (0,3000))
    hidedecorations!(ax)

    ax=Axis(f[1,2],title = "T1",aspect=1)
    h=heatmap!(ax,T1*1000,colorrange = (0,3000))
    hidedecorations!(ax)
    Colorbar(f[1,3],h)

    ax=Axis(f[2,2],title = "(T1 - T1_GT)",aspect=1)
    h=heatmap!(ax,(T1*1000 .- T1_GT),colormap=:bluesreds)
    hidedecorations!(ax)
    Colorbar(f[2,3],h)

    f
end

[cncastillo]

Is T1_GT = reverse(rotl90(T1_GT),dims=1) not the same as T1_GT = T1_GT' (transpose)?

[aTrotier]

Is T1_GT = reverse(rotl90(T1_GT),dims=1) not the same as T1_GT = T1_GT' (transpose)?

Except for complex data ? After checking seems I should use transpose (in that case I don't care but I can rapidly make a mistake with ')

[cncastillo]

True!, better to use transpose, or something that just permutes the dimension for 3D

[cncastillo]

It looks really good! 😮

I could have made a mistake in positioning the center of the recon, and a shift in shiftx = dx_recon/2 is needed to match the reconstructed FOV

dx_recon, dy_recon = FOVx_recon/Nx, FOVy_recon/Ny
shiftx, shifty = dx_recon/2, dy_recon/2
x_recon = Base.range(-FOVx_recon/2 + shiftx, FOVx_recon/2 + shiftx, length=Nx)
y_recon = Base.range(-FOVy_recon/2 + shifty, FOVy_recon/2 + shifty, length=Ny)' 

[aTrotier]

dx_recon, dy_recon = FOVx_recon/Nx, FOVy_recon/Ny
shiftx, shifty = -dx_recon/2, -dy_recon/2
x_recon = Base.range(-FOVx_recon/2 + shiftx, FOVx_recon/2 + shiftx, length=Nx)
y_recon = Base.range(-FOVy_recon/2 + shifty, FOVy_recon/2 + shifty, length=Ny)' 

T1_GT = itp.(x_recon, y_recon)

# rotate
T1_GT=transpose(T1_GT)

gives :

of note difference is now in %

For the example I will try to change the acq in order to perform the reconstruction with the FFTOp

Anyway thanks @cncastillo. I have the proof of concept if I want to build a master course about Quantitative MRI :)

[cncastillo]

Amazing!! Do not hesitate to ask any questions. I would be happy to help with whatever is needed!