Skip to content

Commit 6db6b5e

Browse files
aTrotieraTrotier
committed
test doc-01
1 parent ebd53d9 commit 6db6b5e

9 files changed

+250
-17
lines changed

‎docs/Artifacts.toml

+7
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
[MP2RAGE_data]
2+
git-tree-sha1 = "04cd4c29bb9e2aeb5384fbc70a9af0e1a37ca369"
3+
lazy = true
4+
5+
[[MP2RAGE_data.download]]
6+
sha256 = "1f1b703c79db66ba6ef620651eca431cb0319d87f1eafa53826cb11a93afe4a8"
7+
url = "https://zenodo.org/records/14051522/files/data.tar.gz"

‎docs/Project.toml

+4
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,7 @@
11
[deps]
2+
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
23
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
4+
HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
5+
LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
6+
Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
37
SEQ_BRUKER_a_MP2RAGE_CS_360 = "4ba9dd0f-9419-4edc-b963-b8b02e27d02f"

‎docs/generate_lit.jl

+22
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
execute = isempty(ARGS) || ARGS[1] == "run"
2+
3+
using Literate
4+
5+
lit = joinpath(@__DIR__, "lit")
6+
src = joinpath(@__DIR__, "src")
7+
gen = joinpath(@__DIR__, "src/generated")
8+
9+
for (root, _, files) in walkdir(lit), file in files
10+
splitext(file)[2] == ".jl" || continue # process .jl files only
11+
ipath = joinpath(root, file)
12+
opath = splitdir(replace(ipath, lit => gen))[1]
13+
Literate.markdown(ipath, opath,
14+
documenter = execute)
15+
#codefence = "```julia" => "```") # codefence force to not execute the code
16+
Literate.notebook(ipath, opath; execute = false)
17+
end
18+
19+
# functions
20+
ismd(f) = splitext(f)[2] == ".md"
21+
pages(folder) =
22+
[joinpath("generated/", folder, f) for f in readdir(joinpath(gen, folder)) if ismd(f)]

‎docs/lit/examples/advanced_reco.jl

+88
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
#---------------------------------------------------------
2+
# # [Compressed-sensing reconstruction](@id 02-CS_reconstruction)
3+
#---------------------------------------------------------
4+
5+
# ## Description
6+
#
7+
# This example describes how to perform a compressed-sensingreconstruction of a CS-2 accelerated acquisition.
8+
9+
# ## Loading Package
10+
using LazyArtifacts # loading data
11+
using SEQ_BRUKER_a_MP2RAGE_CS_360
12+
using CairoMakie # plotting
13+
14+
# In addition we load the package internally used to perform the reconstruction
15+
16+
17+
18+
# ## Loading Package
19+
using Artifacts
20+
using LazyArtifacts # loading data
21+
using SEQ_BRUKER_a_MP2RAGE_CS_360
22+
using CairoMakie # plotting
23+
24+
# ## Download the datasets
25+
artifact_toml = "../../../../Artifacts.toml"
26+
_hash = artifact_hash("MP2RAGE_data", artifact_toml)
27+
if ~isnothing(_hash)
28+
datadir = artifact_path(_hash)
29+
else
30+
datadir = artifact_path(artifact_hash("MP2RAGE_data",find_artifacts_toml(".")))
31+
end
32+
@info "The test data is located at $datadir."
33+
34+
# If you want to perform your own reconstruction, you can change the following line in order to point to another a bruker dataset
35+
path_bruker = joinpath(datadir, "MP2RAGE_CS2")
36+
37+
# ## Compressed-sensing reconstruction
38+
# In order to use an advanced reconstruction we will pass some parameters that will be used by the reconstruction package MRIReco.jl
39+
using SEQ_BRUKER_a_MP2RAGE_CS_360.MRIReco
40+
using SEQ_BRUKER_a_MP2RAGE_CS_360.MRIReco.RegularizedLeastSquares
41+
42+
# We have to create a parameter dictionnary that will be used. If you need more information about it take a look at [MRIReco.jl](https://github.com/MagneticResonanceImaging/MRIReco.jl)
43+
44+
CS = Dict{Symbol,Any}()
45+
CS[:sparseTrafo] = "Wavelet" #sparse trafo
46+
CS[:reg] = L1Regularization(100.) # regularization
47+
CS[:solver] = FISTA # solver
48+
CS[:iterations] = 30
49+
50+
d = reconstruction_MP2RAGE(path_bruker; mean_NR=true,paramsCS = CS)
51+
52+
# for comparison purpose let's perform the undersampled reconstruction (without the paramCS keyword)
53+
d_under = reconstruction_MP2RAGE(path_bruker; mean_NR=true)
54+
55+
56+
# We can check the results
57+
58+
begin
59+
f = Figure(size=(500,400))
60+
ax=Axis(f[1,1],title="TI₁ undersampled")
61+
h=heatmap!(ax,abs.(d_under["im_reco"][:,:,60,1,1,1]),colormap=:grays)
62+
63+
ax=Axis(f[1,2],title="TI₁ CS")
64+
h=heatmap!(ax,abs.(d["im_reco"][:,:,60,1,1,1]),colormap=:grays)
65+
66+
67+
ax=Axis(f[2,1],title="UNIT1 undersampled")
68+
h=heatmap!(ax,d_under["T1map"][:,:,60,1,1],colorrange = (500,2000))
69+
70+
ax=Axis(f[2,2],title="UNIT1 CS")
71+
h=heatmap!(ax,d["T1map"][:,:,60,1,1],colorrange = (500,2000))
72+
73+
for ax in f.content # hide decoration befor adding colorbar
74+
hidedecorations!(ax)
75+
end
76+
77+
Colorbar(f[2,3],h,label = "T₁ [ms]", flip_vertical_label=true)
78+
f
79+
end
80+
81+
82+
83+
# ## Write results in BIDS format
84+
# Results can be written following most of the [qBIDS format recommandation](https://bids-specification.readthedocs.io/en/stable/appendices/qmri.html)
85+
86+
subject_name = "sub_01_cs"
87+
dir_path = "" # directory path where the files will be create
88+
write_bids_MP2RAGE(d,subject_name,dir_path)

‎docs/lit/examples/simple_reco.jl

+106
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,106 @@
1+
#---------------------------------------------------------
2+
# # [Simple reconstruction](@id 01-simple_reconstruction)
3+
#---------------------------------------------------------
4+
5+
# ## Description
6+
#
7+
# This example describes how to perform a reconstruction of a fully acquisition acquired with the a_MP2RAGE_CS_360 sequence.
8+
#
9+
10+
# ## Loading Package
11+
using Artifacts
12+
using LazyArtifacts # loading data
13+
using SEQ_BRUKER_a_MP2RAGE_CS_360
14+
using CairoMakie # plotting
15+
16+
# ## Download the datasets
17+
artifact_toml = "../../../../Artifacts.toml"
18+
_hash = artifact_hash("MP2RAGE_data", artifact_toml)
19+
if ~isnothing(_hash)
20+
datadir = artifact_path(_hash)
21+
else
22+
datadir = artifact_path(artifact_hash("MP2RAGE_data",find_artifacts_toml(".")))
23+
end
24+
@info "The test data is located at $datadir."
25+
26+
# If you want to perform your own reconstruction, you can change the following line in order to point to another a bruker dataset
27+
path_bruker = joinpath(datadir, "MP2RAGE_FULLY")
28+
29+
# ## Perform the reconstruction
30+
# this function will perform a standard reconstruction without compressed-sensing. If your data are subsampled, results will be undersampled reconstruction.
31+
#
32+
# the keyword mean_NR=true will average the image before performing the MP2RAGE/T1 maps estimation.
33+
# Otherwise an image/T₁ map will be generated for each Number Of Repetition (NR)
34+
d = reconstruction_MP2RAGE(path_bruker; mean_NR=true)
35+
36+
37+
# the result is a dictionnary with the following fields :
38+
# - "im_reco" : (x,y,z,Number of Repetition,TI) Complex
39+
# - "MP2RAGE" : (x,y,z,TI) Float
40+
# - "T1map" : (x,y,z,Number of Repetition) Float
41+
# - "params_prot"
42+
# - "params_reco"
43+
# - "params_MP2RAGE"
44+
#
45+
# im_reco corresponds to the TI₁ and \TI₂ images in the complex format with 6 dimensions :
46+
# (x,y,z,Number of Repetition,TI)
47+
48+
49+
# We can check the results
50+
51+
begin
52+
f = Figure(size=(500,400))
53+
ax=Axis(f[1,1],title="TI₁")
54+
h=heatmap!(ax,abs.(d["im_reco"][:,:,60,1,1,1]),colormap=:grays)
55+
56+
ax=Axis(f[1,2],title="TI₂")
57+
h=heatmap!(ax,abs.(d["im_reco"][:,:,60,1,1,2]),colormap=:grays)
58+
59+
ax=Axis(f[2,1],title="UNIT1 / MP2RAGE")
60+
h=heatmap!(ax,d["MP2RAGE"][:,:,60,1,1],colormap=:grays)
61+
62+
ax=Axis(f[2,2],title="UNIT1 / MP2RAGE")
63+
h=heatmap!(ax,d["T1map"][:,:,60,1,1],colorrange = (500,2000))
64+
65+
for ax in f.content # hide decoration befor adding colorbar
66+
hidedecorations!(ax)
67+
end
68+
69+
Colorbar(f[2,3],h,label = "T₁ [ms]", flip_vertical_label=true)
70+
f
71+
end
72+
73+
# The Lookup table used for the reconstruction is stored in the dictionnary (LUT)
74+
# First columns is the range of T1.
75+
f=Figure()
76+
ax = Axis(f[1,1],xlabel="T₁ [ms]")
77+
lines!(ax,d["LUT"])
78+
f
79+
80+
# ## Write results in BIDS format
81+
# Results can be written following most of the [qBIDS format recommandation](https://bids-specification.readthedocs.io/en/stable/appendices/qmri.html)
82+
83+
subject_name = "sub_01"
84+
dir_path = "" # directory path where the files will be create
85+
write_bids_MP2RAGE(d,subject_name,dir_path)
86+
87+
#=
88+
which results in :
89+
```
90+
sub_01/
91+
├─ MP2RAGE.json
92+
└─ anat/
93+
├─ sub_01_T1map.nii.gz
94+
├─ sub_01_UNIT1.nii.gz
95+
├─ sub_01_inv-1-complex_MP2RAGE.nii.gz
96+
├─ sub_01_inv-1-mag_MP2RAGE.nii.gz
97+
├─ sub_01_inv-1-phase_MP2RAGE.nii.gz
98+
├─ sub_01_inv-2-complex_MP2RAGE.nii.gz
99+
├─ sub_01_inv-2-mag_MP2RAGE.nii.gz
100+
└─ sub_01_inv-2-phase_MP2RAGE.nii.gz
101+
```
102+
103+
If you want to generate the T1 map with another tools like qMRLab
104+
the required MP2RAGE parameters are stored in the **MP2RAGE.json** file.
105+
=#
106+

‎docs/make.jl

+5-1
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
11
using SEQ_BRUKER_a_MP2RAGE_CS_360
2-
using Documenter
2+
using Documenter, Literate
3+
4+
include("generate_lit.jl")
35

46
DocMeta.setdocmeta!(SEQ_BRUKER_a_MP2RAGE_CS_360, :DocTestSetup, :(using SEQ_BRUKER_a_MP2RAGE_CS_360); recursive=true)
57

@@ -14,6 +16,8 @@ makedocs(;
1416
),
1517
pages=[
1618
"Home" => "index.md",
19+
"Examples" =>["generated/examples/simple_reco.md",
20+
"generated/examples/advanced_reco.md"],
1721
],
1822
)
1923

‎docs/src/index.md

+1-1
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
CurrentModule = SEQ_BRUKER_a_MP2RAGE_CS_360
33
```
44

5-
# SEQ_BRUKER_a_MP2RAGE_CS_360
5+
# SEQ\_BRUKER\_a\_MP2RAGE\_CS\_360
66

77
Documentation for [SEQ_BRUKER_a_MP2RAGE_CS_360](https://github.com/aTrotier/SEQ_BRUKER_a_MP2RAGE_CS_360.jl).
88

‎src/BIDS.jl

+9-9
Original file line numberDiff line numberDiff line change
@@ -51,14 +51,14 @@ function write_bids_MP2RAGE(d::Dict,subname::AbstractString,folder="")
5151
"_UNIT1",
5252
"_T1map"]
5353

54-
data_ = [ abs.(d["im_reco"][:,:,:,1,1,:]),
55-
angle.(d["im_reco"][:,:,:,1,1,:]),
56-
d["im_reco"][:,:,:,1,1,:],
57-
abs.(d["im_reco"][:,:,:,2,1,:]),
58-
angle.(d["im_reco"][:,:,:,2,1,:]),
59-
d["im_reco"][:,:,:,2,1,:],
54+
data_ = [ abs.(d["im_reco"][:,:,:,1,:,1]),
55+
angle.(d["im_reco"][:,:,:,1,:,1]),
56+
d["im_reco"][:,:,:,1,:,1],
57+
abs.(d["im_reco"][:,:,:,1,:,2]),
58+
angle.(d["im_reco"][:,:,:,1,:,2]),
59+
d["im_reco"][:,:,:,1,:,2],
6060
d["MP2RAGE"],
61-
d["T1maps"]]
61+
d["T1map"]]
6262

6363
voxel_size = tuple(parse.(Float64,d["params_prot"]["PVM_SpatResol"])...) #mm
6464
for (name,data) in zip(path_type, data_)
@@ -68,8 +68,8 @@ function write_bids_MP2RAGE(d::Dict,subname::AbstractString,folder="")
6868

6969

7070
# pass parameters
71-
b["ACQ_operator"] # required to read ACQ
72-
MagneticField = parse(Float64,b["BF1"]) / 42.576
71+
d["params_prot"]["ACQ_operator"] # required to read ACQ
72+
MagneticField = parse(Float64, d["params_prot"]["BF1"]) / 42.576
7373
p_MP2 = d["params_MP2RAGE"]
7474

7575
# define JSON dict

‎src/reconstruction.jl

+8-6
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
export reconstruction_MP2RAGE
22

3-
function reconstruction_MP2RAGE(path_bruker;mean_NR::Bool = false)
3+
function reconstruction_MP2RAGE(path_bruker;mean_NR::Bool = false,paramsCS=Dict())
44
b = BrukerFile(path_bruker)
55
raw = RawAcquisitionData_MP2RAGE(b)
66
acq = AcquisitionData(raw,OffsetBruker=true)
@@ -18,27 +18,29 @@ function reconstruction_MP2RAGE(path_bruker;mean_NR::Bool = false)
1818
params[:senseMaps] = sens_fully
1919
params[:iterations] = 1
2020

21+
params = merge(params,paramsCS)
22+
2123
x_approx = reconstruction(acq, params).data
22-
x_approx = permutedims(x_approx,(1,2,3,5,6,4))
2324
if mean_NR
2425
x_approx = mean(x_approx,dims=6)
2526
end
27+
x_approx = permutedims(x_approx,(1,2,3,5,6,4))
2628

2729
## process data to extract T1 maps
28-
MP2 = mp2rage_comb(x_approx[:,:,:,:,1,:])
30+
MP2 = mp2rage_comb(x_approx[:,:,:,:,:,:])
2931

3032
p = params_from_seq(b)
3133

32-
T1,range_T1 = mp2rage_T1maps(MP2,p)
34+
T1,range_T1,LUT = mp2rage_T1maps(MP2,p)
3335

3436
d = Dict{Any,Any}()
3537
d["im_reco"] = x_approx
3638
d["MP2RAGE"] = MP2
37-
d["T1maps"] = T1
38-
d["range_T1"] = range_T1
39+
d["T1map"] = T1
3940
d["params_reco"] = params
4041
d["params_MP2RAGE"] = p
4142
d["params_prot"] = b
43+
d["LUT"] = [range_T1;;LUT]
4244

4345
return d
4446
end

0 commit comments

Comments
 (0)