Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Custom values for brain_phantom #484

Open
wants to merge 23 commits into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
271d98f
Custom values for brain_phantom
gsahonero Sep 22, 2024
a67c5bd
Update Phantom.jl
gsahonero Sep 23, 2024
afbd05b
Introducing default brain phantom properties
gsahonero Sep 26, 2024
262bf14
Merge branch 'master' of https://github.com/gsahonero/KomaMRI.jl
gsahonero Sep 26, 2024
cf83a9e
Details on documentation
gsahonero Sep 26, 2024
f858ef1
Updating to remove unnecessary things
gsahonero Sep 26, 2024
0294f4b
Returning phantom to original code
gsahonero Sep 26, 2024
165b793
Resolving problems in brain phantom
gsahonero Sep 27, 2024
91d59f5
Merge branch 'master' of https://github.com/gsahonero/KomaMRI.jl
gsahonero Sep 27, 2024
6b990f4
Merge branch 'JuliaHealth:master' into master
gsahonero Sep 27, 2024
5b6e25d
Avoiding flat arrays.
gsahonero Oct 4, 2024
4020c3d
Merge branch 'JuliaHealth:master' into master
gsahonero Oct 4, 2024
4e4896b
Merge branch 'master' into master
cncastillo Nov 1, 2024
3689518
Update Phantom.jl
gsahonero Nov 9, 2024
9bfe36b
Merge branch 'master' into master
cncastillo Nov 13, 2024
9cc9d29
Update KomaMRIBase/src/datatypes/Phantom.jl
gsahonero Nov 15, 2024
1b9aa93
Update KomaMRIBase/src/datatypes/Phantom.jl
gsahonero Nov 15, 2024
bd6b91e
Update KomaMRIBase/src/datatypes/Phantom.jl
gsahonero Nov 15, 2024
f55d387
Update KomaMRIBase/src/datatypes/Phantom.jl
gsahonero Nov 15, 2024
ecc9bbe
Update KomaMRIBase/src/datatypes/Phantom.jl
gsahonero Nov 15, 2024
1079053
Requested changes
gsahonero Dec 11, 2024
b843b60
Removing a println
gsahonero Dec 11, 2024
e61cd85
Reducing lines
gsahonero Dec 12, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
224 changes: 101 additions & 123 deletions KomaMRIBase/src/datatypes/Phantom.jl
Original file line number Diff line number Diff line change
Expand Up @@ -212,7 +212,7 @@ Default ss=4 sample spacing is 2 mm. Original file (ss=1) sample spacing is .5 m
- `axis`: (`::String`, `="axial"`, opts=[`"axial"`, `"coronal"`, `"sagittal"`]) orientation of the phantom
- `ss`: (`::Integer or ::Vector{Integer}`, `=4`) subsampling parameter for all axes if scaler, per axis if 2 element vector [ssx, ssy]
- `us`: (`::Integer or ::Vector{Integer}`, `=1`) upsampling parameter for all axes if scaler, per axis if 2 element vector [usx, usy], if used ss is set to ss=1

- `tissue_properties`: (`::Dict`, `=Dict()`) phantom tissue properties in SI units considering the available tissues

# Returns
- `obj`: (`::Phantom`) Phantom struct
Expand All @@ -223,10 +223,26 @@ julia> obj = brain_phantom2D(; axis="sagittal", ss=1)

julia> obj = brain_phantom2D(; axis="axial", us=[1, 2])

julia> phantom_values =
Dict(
# ρ, T1, T2, T2*, Δw
"CSF" => [1, 2.569, 0.329, 0.058, 0],
"GM" => [0.86, 0.833, 0.083, 0.069, 0],
"WM" => [0.77, 0.500, 0.070, 0.061, 0],
"FAT1" => [0, 0, 0, 0, 0],
"MUSCLE" => [0, 0, 0, 0, 0],
"SKIN/MUSCLE" => [0, 0, 0, 0, 0],
"SKULL" => [0, 0, 0, 0, 0],
"VESSELS" => [0, 0, 0, 0, 0],
"FAT2" => [0, 0, 0, 0, 0],
"DURA" => [0, 0, 0, 0, 0],
"MARROW" => [0, 0, 0, 0, 0])
julia> obj = brain_phantom2D(; tissue_properties=phantom_values)

julia> plot_phantom_map(obj, :ρ)
```
"""
function brain_phantom2D(; axis="axial", ss=4, us=1)
function brain_phantom2D(; axis="axial", ss=4, us=1, tissue_properties = Dict())
# check and filter input
ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(2, ss, us)

Expand All @@ -235,75 +251,25 @@ function brain_phantom2D(; axis="axial", ss=4, us=1)
data = MAT.matread(path * "/phantom/brain2D.mat")

# subsample or upsample the phantom data
class = repeat(data[axis][1:ssx:end, 1:ssy:end]; inner=[usx, usy])
labels = repeat(data[axis][1:ssx:end, 1:ssy:end]; inner=[usx, usy])
# to make it compatible with default_brain_tissue_properties
labels = reshape(labels, (size(labels, 1), size(labels, 2), 1))

# Define spin position vectors
Δx = .5e-3 * ssx / usx
Δy = .5e-3 * ssy / usy
M, N = size(class)
M, N = size(labels)
FOVx = (M - 1) * Δx #[m]
FOVy = (N - 1) * Δy #[m]
x = (-FOVx / 2):Δx:(FOVx / 2) #spin coordinates
y = (-FOVy / 2):Δy:(FOVy / 2) #spin coordinates
x, y = x .+ y' * 0, x * 0 .+ y' #grid points
x = reshape(x, (size(x, 1), size(x, 2), 1))
y = reshape(y, (size(y, 1), size(y, 2), 1))

# Define spin property vectors
T2 =
(class .== 23) * 329 .+ #CSF
(class .== 46) * 83 .+ #GM
(class .== 70) * 70 .+ #WM
(class .== 93) * 70 .+ #FAT1
(class .== 116) * 47 .+ #MUSCLE
(class .== 139) * 329 .+ #SKIN/MUSCLE
(class .== 162) * 0 .+ #SKULL
(class .== 185) * 0 .+ #VESSELS
(class .== 209) * 70 .+ #FAT2
(class .== 232) * 329 .+ #DURA
(class .== 255) * 70 #MARROW
T2s =
(class .== 23) * 58 .+ #CSF
(class .== 46) * 69 .+ #GM
(class .== 70) * 61 .+ #WM
(class .== 93) * 58 .+ #FAT1
(class .== 116) * 30 .+ #MUSCLE
(class .== 139) * 58 .+ #SKIN/MUSCLE
(class .== 162) * 0 .+ #SKULL
(class .== 185) * 0 .+ #VESSELS
(class .== 209) * 61 .+ #FAT2
(class .== 232) * 58 .+ #DURA
(class .== 255) * 61 #MARROW
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
ρ =
(class .== 23) * 1 .+ #CSF
(class .== 46) * 0.86 .+ #GM
(class .== 70) * 0.77 .+ #WM
(class .== 93) * 1 .+ #FAT1
(class .== 116) * 1 .+ #MUSCLE
(class .== 139) * 1 .+ #SKIN/MUSCLE
(class .== 162) * 0 .+ #SKULL
(class .== 185) * 0 .+ #VESSELS
(class .== 209) * 0.77 .+ #FAT2
(class .== 232) * 1 .+ #DURA
(class .== 255) * 0.77 #MARROW
Δw_fat = -220 * 2π
Δw =
(class .== 93) * Δw_fat .+ #FAT1
(class .== 209) * Δw_fat #FAT2
T1 = T1 * 1e-3
T2 = T2 * 1e-3
T2s = T2s * 1e-3

# Get tissue properties
ρ, T1, T2, T2s, Δw = default_brain_tissue_properties(labels, tissue_properties)

# Define and return the Phantom struct
obj = Phantom{Float64}(;
name="brain2D_" * axis,
Expand Down Expand Up @@ -336,6 +302,7 @@ Default ss=4 sample spacing is 2 mm. Original file (ss=1) sample spacing is .5 m
- `ss`: (`::Integer or ::Vector{Integer}`, `=4`) subsampling parameter for all axes if scaler, per axis if 3 element vector [ssx, ssy, ssz]
- `us`: (`::Integer or ::Vector{Integer}`, `=1`) upsampling parameter for all axes if scaler, per axis if 3 element vector [usx, usy, usz]
- `start_end`: (`::Vector{Integer}`, `=[160,200]`) z index range of presampled phantom, 180 is center
- `tissue_properties`: (`::Dict`, `=Dict()`) phantom tissue properties in SI units considering the available tissues

# Returns
- `obj`: (`::Phantom`) 3D Phantom struct
Expand All @@ -346,19 +313,34 @@ julia> obj = brain_phantom3D(; ss=5)

julia> obj = brain_phantom3D(; us=[2, 2, 1])

julia> phantom_values =
Dict(
# ρ, T1, T2, T2*, Δw
"CSF" => [1, 2.569, 0.329, 0.058, 0],
"GM" => [0.86, 0.833, 0.083, 0.069, 0],
"WM" => [0.77, 0.500, 0.070, 0.061, 0],
"FAT1" => [0, 0, 0, 0, 0],
"MUSCLE" => [0, 0, 0, 0, 0],
"SKIN/MUSCLE" => [0, 0, 0, 0, 0],
"SKULL" => [0, 0, 0, 0, 0],
"VESSELS" => [0, 0, 0, 0, 0],
"FAT2" => [0, 0, 0, 0, 0],
"DURA" => [0, 0, 0, 0, 0],
"MARROW" => [0, 0, 0, 0, 0])
julia> obj = brain_phantom3D(; tissue_properties=phantom_values)

julia> plot_phantom_map(obj, :ρ)
```
"""
function brain_phantom3D(; ss=4, us=1, start_end=[160, 200])
function brain_phantom3D(; ss=4, us=1, start_end=[160, 200], tissue_properties=Dict())
# check and filter input
ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(3, ss, us)

# Get data from .mat file
path = @__DIR__
data = MAT.matread(path * "/phantom/brain3D.mat")

# subsample or upsample the phantom data
class = repeat(
labels = repeat(
data["data"][1:ssx:end, 1:ssy:end, start_end[1]:ssz:start_end[2]];
inner=[usx, usy, usz],
)
Expand All @@ -367,7 +349,7 @@ function brain_phantom3D(; ss=4, us=1, start_end=[160, 200])
Δx = .5e-3 * ssx / usx
Δy = .5e-3 * ssy / usy
Δz = .5e-3 * ssz / usz
M, N, Z = size(class)
M, N, Z = size(labels)
FOVx = (M - 1) * Δx #[m]
FOVy = (N - 1) * Δy #[m]
FOVz = (Z - 1) * Δz #[m]
Expand All @@ -377,63 +359,9 @@ function brain_phantom3D(; ss=4, us=1, start_end=[160, 200])
x = 1 * xx .+ 0 * yy .+ 0 * zz
y = 0 * xx .+ 1 * yy .+ 0 * zz
z = 0 * xx .+ 0 * yy .+ 1 * zz

# Define spin property vectors
T2 =
(class .== 23) * 329 .+ #CSF
(class .== 46) * 83 .+ #GM
(class .== 70) * 70 .+ #WM
(class .== 93) * 70 .+ #FAT1
(class .== 116) * 47 .+ #MUSCLE
(class .== 139) * 329 .+ #SKIN/MUSCLE
(class .== 162) * 0 .+ #SKULL
(class .== 185) * 0 .+ #VESSELS
(class .== 209) * 70 .+ #FAT2
(class .== 232) * 329 .+ #DURA
(class .== 255) * 70 #MARROW
T2s =
(class .== 23) * 58 .+ #CSF
(class .== 46) * 69 .+ #GM
(class .== 70) * 61 .+ #WM
(class .== 93) * 58 .+ #FAT1
(class .== 116) * 30 .+ #MUSCLE
(class .== 139) * 58 .+ #SKIN/MUSCLE
(class .== 162) * 0 .+ #SKULL
(class .== 185) * 0 .+ #VESSELS
(class .== 209) * 61 .+ #FAT2
(class .== 232) * 58 .+ #DURA
(class .== 255) * 61 #MARROW
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
ρ =
(class .== 23) * 1 .+ #CSF
(class .== 46) * 0.86 .+ #GM
(class .== 70) * 0.77 .+ #WM
(class .== 93) * 1 .+ #FAT1
(class .== 116) * 1 .+ #MUSCLE
(class .== 139) * 1 .+ #SKIN/MUSCLE
(class .== 162) * 0 .+ #SKULL
(class .== 185) * 0 .+ #VESSELS
(class .== 209) * 0.77 .+ #FAT2
(class .== 232) * 1 .+ #DURA
(class .== 255) * 0.77 #MARROW
Δw_fat = -220 * 2π
Δw =
(class .== 93) * Δw_fat .+ #FAT1
(class .== 209) * Δw_fat #FAT2
T1 = T1 * 1e-3
T2 = T2 * 1e-3
T2s = T2s * 1e-3

# Get tissue properties
ρ, T1, T2, T2s, Δw = default_brain_tissue_properties(labels, tissue_properties)

# Define and return the Phantom struct
obj = Phantom{Float64}(;
Expand Down Expand Up @@ -482,7 +410,6 @@ function pelvis_phantom2D(; ss=4, us=1)

# subsample or upsample the phantom data
class = repeat(data["pelvis3D_slice"][1:ssx:end, 1:ssy:end]; inner=[usx, usy])

# Define spin position vectors
Δx = .5e-3 * ssx / usx
Δy = .5e-3 * ssy / usy
Expand Down Expand Up @@ -612,5 +539,56 @@ function check_phantom_arguments(nd, ss, us)
ssy = ss[2]
end
end

return ssx, ssy, ssz, usx, usy, usz
end

"""
ρ, T1, T2, T2s, Δw = default_brain_tissue_properties(labels, tissue_properties = nothing)

This function returns the default brain tissue properties using a labels identifier Matrix
# Arguments
- `labels` : (`::Matrix`) the labels identifier matrix of the phantom
- `tissue_properties` : (`::Dict`, `=Dict()`) phantom tissue properties in ms and Hz considering the available tissues

# Returns
- `ρ, T1, T2, T2s, Δw`: (`::Matrix`) matrices of the same size of labels with the tissues properties information

# Examples
```julia-repl
julia> ρ, T1, T2, T2s, Δw = default_brain_tissue_properties(labels, tissue_properties)

julia> ρ, T1, T2, T2s, Δw = default_brain_tissue_properties(labels)
```
"""
function default_brain_tissue_properties(labels, tissue_properties = Dict())
# Load default tissue properties
default_properties = Dict(
# ρ, T1, T2, T2*, Δw
"CSF" => [1, 2.569, 0.329, 0.058, 0],
"GM" => [0.86, 0.833, 0.083, 0.069, 0],
"WM" => [0.77, 0.500, 0.070, 0.061, 0],
"FAT1" => [1, 0.350, 0.070, 0.058, -220*2π], #-220 Hz
"MUSCLE" => [1, 0.900, 0.047, 0.030, 0],
"SKIN/MUSCLE" => [1, 0.569, 0.329, 0.058, 0],
"SKULL" => [0, 0, 0, 0, 0],
"VESSELS" => [0, 0, 0, 0, 0],
"FAT2" => [0.77, 0.500, 0.070, 0.061, -220*2π], #-220 Hz
"DURA" => [1, 2.569, 0.329, 0.058, 0],
"MARROW" => [0.77, 0.500, 0.070, 0.061, 0])

tissue_properties = merge(default_properties, tissue_properties)
props = ["ρ", "T1", "T2", "T2s", "Δw"]
Nproperties = length(props)
# Order: CSF, DURA, FAT1, FAT2, GM, MARROW, MUSCLE, SKIN/MUSCLE, SKULL, vESSELS, WM
tissue_labels = [23, 232, 93, 209, 46, 255, 116, 139, 162, 185, 70]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you try one last thing?

The sorting and use of tissue_texts, does not make a lot of sense. It is better if tissue_labels is also a dictionary (Dict("CSF" => 23, ...)), and we index it by the default_properties keys.

You should be able to do:

for i=1:Nproperties
    for (tissue, property) in tissue_properties
        data_properties[i, :, :, :] += (labels .== tissue_labels[tissue]) * property[i]
    end
end

Also, don't forget the spaces between the *!

tissue_texts = sort(collect(keys(default_properties)))#
data_properties = zeros(Nproperties, size(labels)...)
for i=1:Nproperties
for (label, tissue) in zip(tissue_labels, tissue_texts)
data_properties[i, :, :, :] += (labels .== label)*tissue_properties[tissue][i]
end
end

return (data_properties[i,:,:,:] for i in 1:Nproperties)
end
Loading