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.

Sign up for GitHub

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

Jump to bottom

Custom values for brain_phantom #484

Open
wants to merge 20 commits into
base: master
Choose a base branch
from
+101 −122
Open

Custom values for brain_phantom #484

Changes from all commits
Commits
Show all changes
20 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
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
223 changes: 101 additions & 122 deletions KomaMRIBase/src/datatypes/Phantom.jl
View file
Open in desktop
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" => [2.569, 0.329, 0.058, 1, 0],
"GM" => [1.153, 0.083, 0.069, 0.86, 0],
"WM" => [0.746, 0.070, 0.061, 0.77, 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,74 +251,20 @@ 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])

# 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

# 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 @@ -336,6 +298,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 +309,34 @@ julia> obj = brain_phantom3D(; ss=5)

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

julia> phantom_values =
Dict(
# T1, T2, T2*, ρ, Δw
"CSF" => [2.569, 0.329, 0.058, 1, 0],
"GM" => [1.153, 0.083, 0.069, 0.86, 0],
"WM" => [0.746, 0.070, 0.061, 0.77, 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 +345,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 +355,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 +406,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 +535,61 @@ 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" => [2.569, 0.329, 0.058, 1, 0],
"GM" => [0.833, 0.083, 0.069, 0.86, 0],
"WM" => [0.500, 0.070, 0.061, 0.77, 0],
"FAT1" => [0.350, 0.070, 0.058, 1, -3.84], #-220 Hz
"MUSCLE" => [0.900, 0.047, 0.030, 1, 0],
"SKIN/MUSCLE" => [0.569, 0.329, 0.058, 1, 0],
"SKULL" => [0, 0, 0, 0, 0],
"VESSELS" => [0, 0, 0, 0, 0],
"FAT2" => [0.500, 0.070, 0.061, 0.77, -3.84], #-220 Hz
"DURA" => [2.569, 0.329, 0.058, 1, 0],
"MARROW" => [0.500, 0.070, 0.061, 0.77, 0])
Comment on lines +562 to +574
Copy link
Member

Choose a reason for hiding this comment

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

Improve indentation, and reorder properties (not done here):

Suggested change
default_properties = Dict(
# T1, T2, T2*, ρ, Δw
"CSF" => [2.569, 0.329, 0.058, 1, 0],
"GM" => [0.833, 0.083, 0.069, 0.86, 0],
"WM" => [0.500, 0.070, 0.061, 0.77, 0],
"FAT1" => [0.350, 0.070, 0.058, 1, -3.84], #-220 Hz
"MUSCLE" => [0.900, 0.047, 0.030, 1, 0],
"SKIN/MUSCLE" => [0.569, 0.329, 0.058, 1, 0],
"SKULL" => [0, 0, 0, 0, 0],
"VESSELS" => [0, 0, 0, 0, 0],
"FAT2" => [0.500, 0.070, 0.061, 0.77, -3.84], #-220 Hz
"DURA" => [2.569, 0.329, 0.058, 1, 0],
"MARROW" => [0.500, 0.070, 0.061, 0.77, 0])
default_properties = Dict(
# T1 T2 T2* ρ Δw
"CSF" => [2.569, 0.329, 0.058, 1, 0],
"GM" => [0.833, 0.083, 0.069, 0.86, 0],
"WM" => [0.500, 0.070, 0.061, 0.77, 0],
"FAT1" => [0.350, 0.070, 0.058, 1, -3.84], # -220 Hz
"MUSCLE" => [0.900, 0.047, 0.030, 1, 0],
"SKIN/MUSCLE" => [0.569, 0.329, 0.058, 1, 0],
"SKULL" => [0, 0, 0, 0, 0],
"VESSELS" => [0, 0, 0, 0, 0],
"FAT2" => [0.500, 0.070, 0.061, 0.77, -3.84], # -220 Hz
"DURA" => [2.569, 0.329, 0.058, 1, 0],
"MARROW" => [0.500, 0.070, 0.061, 0.77, 0]
)


tissue_properties = merge(default_properties, tissue_properties)
properties = []
for i=1:5
temp =
(labels .== 23) * tissue_properties["CSF"][i] .+ #CSF
(labels .== 46) * tissue_properties["GM"][i] .+ #GM
(labels .== 70) * tissue_properties["WM"][i] .+ #WM
(labels .== 93) * tissue_properties["FAT1"][i] .+ #FAT1
(labels .== 116) * tissue_properties["MUSCLE"][i] .+ #MUSCLE
(labels .== 139) * tissue_properties["SKIN/MUSCLE"][i] .+ #SKIN/MUSCLE
(labels .== 162) * tissue_properties["SKULL"][i] .+ #SKULL
(labels .== 185) * tissue_properties["VESSELS"][i] .+ #VESSELS
(labels .== 209) * tissue_properties["FAT2"][i] .+ #FAT2
(labels .== 232) * tissue_properties["DURA"][i] .+ #DURA
(labels .== 255) * tissue_properties["MARROW"][i] #MARROW
push!(properties, temp)
end
Comment on lines +578 to +592
Copy link
Member

Choose a reason for hiding this comment

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

Just realized this is iterating for (T1, T2, T2*, ρ, Δw) instead of the tissue properties. This code has many problems:

  • for loop is hard-coded to i = 1:5, what is 5? what if we want more properties? if we change it to 6 it will error.
  • the properties are not preallocated, you push! them.
  • The order of the properties does not make much sense. T1 and T2 before ρ?
  • The Δw = -3.84 doesn't seem correct.

Suggestion:

Nproperties = 5 # ρ, T1, T2, T2*, Δw, explicit ... 5 is the number of properties
# or
Nproperties = size(tissue_properties ... , 2) # Or something like that
properties = zeros(size(...), Nproperties)
tissue_labels = [23, 46, ....]
for (label, tissue) in zip(tissue_labels, tissue_properties)
	for property in 1:Nproperties
		properties[:,:, property] .+= (labels .== label) .* tissue[property]
	end
end


return properties
end
Loading