diff --git a/KomaMRIBase/src/datatypes/Phantom.jl b/KomaMRIBase/src/datatypes/Phantom.jl index 5f155f28e..c804caeda 100644 --- a/KomaMRIBase/src/datatypes/Phantom.jl +++ b/KomaMRIBase/src/datatypes/Phantom.jl @@ -115,7 +115,7 @@ function get_dims(obj::Phantom) push!(dims, any(x -> x != 0, obj.x)) push!(dims, any(x -> x != 0, obj.y)) push!(dims, any(x -> x != 0, obj.z)) - return dims + return sum(dims) > 0 ? dims : Bool[1, 0, 0] end """ diff --git a/KomaMRIBase/src/motion/motionlist/Motion.jl b/KomaMRIBase/src/motion/motionlist/Motion.jl index d035b562b..869d0e06c 100644 --- a/KomaMRIBase/src/motion/motionlist/Motion.jl +++ b/KomaMRIBase/src/motion/motionlist/Motion.jl @@ -32,6 +32,20 @@ julia> motion = Motion( spins ::AbstractSpinSpan = AllSpins() end +# Main constructors +function Motion(action) + T = first(typeof(action).parameters) + return Motion(action, TimeRange(zero(T)), AllSpins()) +end +function Motion(action, time::AbstractTimeSpan) + T = first(typeof(action).parameters) + return Motion(action, time, AllSpins()) +end +function Motion(action, spins::AbstractSpinSpan) + T = first(typeof(action).parameters) + return Motion(action, TimeRange(zero(T)), spins) +end + # Custom constructors """ translate = Translate(dx, dy, dz, time, spins) diff --git a/KomaMRIBase/src/motion/motionlist/MotionList.jl b/KomaMRIBase/src/motion/motionlist/MotionList.jl index 75e6460a2..b75841f7c 100644 --- a/KomaMRIBase/src/motion/motionlist/MotionList.jl +++ b/KomaMRIBase/src/motion/motionlist/MotionList.jl @@ -32,7 +32,7 @@ struct MotionList{T<:Real} <: AbstractMotion{T} end """ Constructors """ -MotionList(motions...) = length([motions]) > 0 ? MotionList([motions...]) : @error "You must provide at least one motion as input argument. If you do not want to define motion, use `NoMotion{T}()`" +MotionList(motions::Motion...) = length([motions]) > 0 ? MotionList([motions...]) : @error "You must provide at least one motion as input argument. If you do not want to define motion, use `NoMotion{T}()`" """ MotionList sub-group """ function Base.getindex(mv::MotionList{T}, p) where {T<:Real} diff --git a/KomaMRIPlots/src/ui/DisplayFunctions.jl b/KomaMRIPlots/src/ui/DisplayFunctions.jl index 81e0bb198..c2aeb29ac 100644 --- a/KomaMRIPlots/src/ui/DisplayFunctions.jl +++ b/KomaMRIPlots/src/ui/DisplayFunctions.jl @@ -1135,9 +1135,21 @@ function plot_phantom_map( l = Layout(;title=obj.name*": "*string(key)) if view_2d # 2D + function get_displayed_dims(v) + if sum(v) == 1 + idx = argmax(v)[1] + return [idx in [1, 3], idx in [1, 2], idx in [2,3]] + else + return v + end + end + # Layout config + dims = get_displayed_dims(KomaMRIBase.get_dims(obj)) + axis = ["x", "y", "z"][dims] + l[:xaxis] = attr( - title="x", + title=axis[1], range=[x0, xf], ticksuffix=" cm", backgroundcolor=plot_bgcolor, @@ -1146,7 +1158,7 @@ function plot_phantom_map( scaleanchor="y" ) l[:yaxis] = attr( - title="y", + title=axis[2], range=[x0, xf], ticksuffix=" cm", backgroundcolor=plot_bgcolor, @@ -1158,8 +1170,8 @@ function plot_phantom_map( # Add traces for i in 1:length(t) push!(traces, scattergl( - x=(x[:,i])*1e2, - y=(y[:,i])*1e2, + x=dims[1] ? (x[:,i])*1e2 : (y[:,i])*1e2, + y=dims[1] & dims[2] ? (y[:,i])*1e2 : (z[:,i])*1e2, mode="markers", marker=attr(color=getproperty(obj,key)*factor, showscale=colorbar, @@ -1244,7 +1256,7 @@ function plot_phantom_map( ) for (i, t0) in enumerate(t) ], - currentvalue_prefix="x = ", + currentvalue_prefix="t = ", currentvalue_suffix="ms", )] l[:margin] = attr(t=50, l=0, r=0) diff --git a/README.md b/README.md index 92bc54cde..8363b666f 100644 --- a/README.md +++ b/README.md @@ -60,7 +60,7 @@ KomaMRI.jl is a Julia package for highly efficient ⚡ MRI simulations. KomaMRI ## News - **(1 Oct 2024)** [KomaMRI v0.9](https://github.com/JuliaHealth/KomaMRI.jl/releases/tag/v0.9.0): device-agnostic simulations, improved performance (**4-5x faster and 80x less memory**), distributed simulations, GPU benchmarking, mix-and-match motion definitions, improved dynamic phantom plotting, and a new phantom file format! -- **(29 Aug 2024)** Our first GSoC student, Ryan Kierulf, presented his fantastic work at the JuliaHealth monthly meeting 🥳! (presentation available [here](https://www.youtube.com/watch?v=R6Z20G0J4bM)) More info in the docs: [GPU Parallelization](https://juliahealth.org/KomaMRI.jl/dev/explanation/4-gpu-explanation/), [Distributed Simulations](https://juliahealth.org/KomaMRI.jl/dev/how-to/4-run-distributed-simulations/) and [Ryan's JuliaHealth blog](https://juliahealth.org/JuliaHealthBlog/posts/ryan-gsoc/Ryan_GSOC.html) +- **(29 Aug 2024)** Our first GSoC student, Ryan Kierulf, presented his fantastic work at the JuliaHealth monthly meeting 🥳! (presentation available [here](https://www.youtube.com/watch?v=R6Z20G0J4bM)) More info in the docs: [GPU Parallelization](https://juliahealth.org/KomaMRI.jl/dev/explanation/7-gpu-explanation/), [Distributed Simulations](https://juliahealth.org/KomaMRI.jl/dev/how-to/4-run-distributed-simulations/) and [Ryan's JuliaHealth blog](https://juliahealth.org/JuliaHealthBlog/posts/ryan-gsoc/Ryan_GSOC.html) - **(7 Dec 2023)** Koma was present in [MRI Together](https://mritogether.esmrmb.org/) 😼. The talk is available [here](https://www.youtube.com/watch?v=9mRQH8um4-A). Also, I uploaded the promised [educational example](https://juliahealth.org/KomaMRI.jl/stable/tutorial-pluto/01-gradient-echo-spin-echo/). - **(17 Nov 2023)** Pretty excited of being part of [ISMRM Pulseq's virtual meeting](https://github.com/pulseq/ISMRM-Virtual-Meeting--November-15-17-2023). The slides available [here](https://github.com/pulseq/ISMRM-Virtual-Meeting--November-15-17-2023/blob/35a8da7eaa0bf42f2127e1338a440ccd4e3ef53c/slides/day3_KomaMRI_simulator_Quantitative_MRI.pdf). - **(27 Jul 2023)** I gave a talk at MIT 😄 for [JuliaCon 2023](https://juliacon.org/2023/)! A video of the presentation can be seen [here](https://www.youtube.com/watch?v=WVT9wJegC6Q). @@ -117,7 +117,7 @@ KomaUI() Press the button that says "Simulate!" to do your first simulation :). Then, a notification will emerge telling you that the simulation was successful. In this notification, you can either select to (1) see the Raw Data or (2) to proceed with the reconstruction. > [!IMPORTANT] -> Starting from **KomaMRI v0.9** we are using [package extensions](https://pkgdocs.julialang.org/v1/creating-packages/#Conditional-loading-of-code-in-packages-(Extensions)) to deal with GPU dependencies, meaning that to run simulations on the GPU, installing (`add CUDA/AMDGPU/Metal/oneAPI`) and loading (`using CUDA/AMDGPU/Metal/oneAPI`) the desired backend will be necessary (see [GPU Parallelization](https://JuliaHealth.github.io/KomaMRI.jl/dev/explanation/4-gpu-explanation) and [Tested compatibility](#tested-compatibility)). +> Starting from **KomaMRI v0.9** we are using [package extensions](https://pkgdocs.julialang.org/v1/creating-packages/#Conditional-loading-of-code-in-packages-(Extensions)) to deal with GPU dependencies, meaning that to run simulations on the GPU, installing (`add CUDA/AMDGPU/Metal/oneAPI`) and loading (`using CUDA/AMDGPU/Metal/oneAPI`) the desired backend will be necessary (see [GPU Parallelization](https://JuliaHealth.github.io/KomaMRI.jl/dev/explanation/7-gpu-explanation) and [Tested compatibility](#tested-compatibility)). ## How to Contribute KomaMRI exists thanks to all our contributors: diff --git a/docs/src/assets/ph-action-types-dark.svg b/docs/src/assets/ph-action-types-dark.svg new file mode 100644 index 000000000..2ef49a7da --- /dev/null +++ b/docs/src/assets/ph-action-types-dark.svg @@ -0,0 +1,824 @@ + + + + diff --git a/docs/src/assets/ph-action-types-light.svg b/docs/src/assets/ph-action-types-light.svg new file mode 100644 index 000000000..9989cb360 --- /dev/null +++ b/docs/src/assets/ph-action-types-light.svg @@ -0,0 +1,817 @@ + + + + diff --git a/docs/src/assets/ph-phantom-file-format-dark.svg b/docs/src/assets/ph-phantom-file-format-dark.svg new file mode 100644 index 000000000..6b2975d03 --- /dev/null +++ b/docs/src/assets/ph-phantom-file-format-dark.svg @@ -0,0 +1,683 @@ + + + + diff --git a/docs/src/assets/ph-phantom-file-format-light.svg b/docs/src/assets/ph-phantom-file-format-light.svg new file mode 100644 index 000000000..6e4e26859 --- /dev/null +++ b/docs/src/assets/ph-phantom-file-format-light.svg @@ -0,0 +1,682 @@ + + + + diff --git a/docs/src/assets/ph-spinspan-types-dark.svg b/docs/src/assets/ph-spinspan-types-dark.svg new file mode 100644 index 000000000..c2f081f4c --- /dev/null +++ b/docs/src/assets/ph-spinspan-types-dark.svg @@ -0,0 +1,200 @@ + + + + diff --git a/docs/src/assets/ph-spinspan-types-light.svg b/docs/src/assets/ph-spinspan-types-light.svg new file mode 100644 index 000000000..7f3839a4c --- /dev/null +++ b/docs/src/assets/ph-spinspan-types-light.svg @@ -0,0 +1,198 @@ + + + + diff --git a/docs/src/assets/ph-timespan-types-dark.svg b/docs/src/assets/ph-timespan-types-dark.svg new file mode 100644 index 000000000..a63c35930 --- /dev/null +++ b/docs/src/assets/ph-timespan-types-dark.svg @@ -0,0 +1,278 @@ + + + + diff --git a/docs/src/assets/ph-timespan-types-light.svg b/docs/src/assets/ph-timespan-types-light.svg new file mode 100644 index 000000000..ba8adbe3f --- /dev/null +++ b/docs/src/assets/ph-timespan-types-light.svg @@ -0,0 +1,278 @@ + + + + diff --git a/docs/src/explanation/3-phantom-format.md b/docs/src/explanation/3-phantom-format.md new file mode 100644 index 000000000..2df137040 --- /dev/null +++ b/docs/src/explanation/3-phantom-format.md @@ -0,0 +1,40 @@ +# Phantom File Format + +## Introduction + +While there is already an open and fairly standardised format for MRI sequences +such as [Pulseq](https://pulseq.github.io/index.html), this is not the case for digital phantoms. +That's why we defined a new ''.phantom'' format, which relies on the [HDF5 standard](https://www.hdfgroup.org/solutions/hdf5/). +HDF5 is specially designed to store large amounts of heterogeneous data and to make it readable +and writable quickly and easily. In addition, it allows the storage of metadata. +For all these reasons, it is the ideal file format for storing phantoms. + +## File Format Specification + +### Phantom File Tree + +```@raw html +
+ +``` + +### Action types + +```@raw html + + +``` + +### TimeSpan types + +```@raw html + + +``` + +### SpinSpan types + +```@raw html + + +``` \ No newline at end of file diff --git a/docs/src/explanation/1-sequence.md b/docs/src/explanation/4-sequence.md similarity index 95% rename from docs/src/explanation/1-sequence.md rename to docs/src/explanation/4-sequence.md index 6f5206e7f..d11c52617 100644 --- a/docs/src/explanation/1-sequence.md +++ b/docs/src/explanation/4-sequence.md @@ -2,7 +2,7 @@ This section delves into some details about how a sequence is constructed. The sequence definition in **KomaMRI** is strongly related to the [Pulseq](https://pulseq.github.io/index.html) definition. After reading this section, you should be able to create your own **Sequence** structs for conducting custom simulations using the **KomaMRI** package. -## Sequence Overview +## KomaMRI Sequence Overview Let's introduce the following simple sequence figure to expand from a visual example to a more general sequence definition: ```@raw html @@ -34,7 +34,7 @@ mutable struct Sequence end ``` -As you can see, a **Sequence** struct contains 5 field names: ''DEF'' contains information for reconstruction steps (so it is not mandatory to fill it), ''DUR'' is a vector that contains the time durations of each block, ''ADC'' is also a vector with the acquisition samples for every block (an vector of **ADC** structs), ''GR'' is a 2D matrix which 3 rows representing the x-y-z gradients and columns having the samples of each block (a matrix of **Grad** structs) and ''RF'' is also a 2D matrix where each row represents a different coil and the columns are for different block samples too (a matrix of **RF** structs). The **RF**, **Grad** and **ADC** are MRI events that will be explained in the section [Events Definitions](2-seq-events.md). +As you can see, a **Sequence** struct contains 5 field names: ''DEF'' contains information for reconstruction steps (so it is not mandatory to fill it), ''DUR'' is a vector that contains the time durations of each block, ''ADC'' is also a vector with the acquisition samples for every block (an vector of **ADC** structs), ''GR'' is a 2D matrix which 3 rows representing the x-y-z gradients and columns having the samples of each block (a matrix of **Grad** structs) and ''RF'' is also a 2D matrix where each row represents a different coil and the columns are for different block samples too (a matrix of **RF** structs). The **RF**, **Grad** and **ADC** are MRI events that will be explained in the section [Events Definitions](5-seq-events.md). !!! warning So far, **KomaMRI** can only manage one coil for RF excitations. However, in future versions, parallel transmit pTX will be managed by adding more ``rows'' to the RF matrix of the Sequence field name. @@ -88,7 +88,7 @@ julia> seq.DUR 0.0004042313086942605 ``` -Additionally, you can access a subset of blocks in a **Sequence** by slicing or indexing. The result will also be a **Sequence** struct, allowing you to perform the same operations as you would with a full Sequence. For example, if you want to analyze the first 11 blocks, you can do the following: +Additionally, you can access a subset of blocks in a **Sequence** by slicing or indexing. The result will also be a **Sequence** struct, allowing you to perform the same operations as you would with a full Sequence (just a heads-up: this is analogous for the [Phantom](1-phantom.md) structure). For example, if you want to analyze the first 11 blocks, you can do the following: ```julia-repl julia> seq[1:11] Sequence[ τ = 3.837 ms | blocks: 11 | ADC: 5 | GR: 11 | RF: 1 | DEF: 5 ] diff --git a/docs/src/explanation/2-seq-events.md b/docs/src/explanation/5-seq-events.md similarity index 100% rename from docs/src/explanation/2-seq-events.md rename to docs/src/explanation/5-seq-events.md diff --git a/docs/src/explanation/3-simulation.md b/docs/src/explanation/6-simulation.md similarity index 99% rename from docs/src/explanation/3-simulation.md rename to docs/src/explanation/6-simulation.md index 5810e8a35..d64f6b4e8 100644 --- a/docs/src/explanation/3-simulation.md +++ b/docs/src/explanation/6-simulation.md @@ -24,7 +24,7 @@ From the programming perspective, it is needed to call the [`simulate`](@ref) fu | `"gpu"` | is a boolean that determines whether to use GPU or CPU hardware resources, as long as they are available on the host computer. | | `"gpu_device"` | sets the index ID of the available GPU in the host computer. | -For instance, if you want to perform a simulation on the CPU with float64 precision using the `BlochDict()` method (assuming you have already defined `obj` and `seq`), you can do so like this: +For instance, if you want to perform a simulation on the CPU with float64 precision using the `BlochDict()` method (assuming you have already defined `obj`, `seq` and `sys`), you can do so like this: ```julia # Set non-default simulation parameters and run simulation sim_params = KomaMRICore.default_sim_params() diff --git a/docs/src/explanation/4-gpu-explanation.md b/docs/src/explanation/7-gpu-explanation.md similarity index 100% rename from docs/src/explanation/4-gpu-explanation.md rename to docs/src/explanation/7-gpu-explanation.md diff --git a/docs/src/explanation/lit-1-phantom.jl b/docs/src/explanation/lit-1-phantom.jl new file mode 100644 index 000000000..540fcc287 --- /dev/null +++ b/docs/src/explanation/lit-1-phantom.jl @@ -0,0 +1,179 @@ +# # Phantom + +using KomaMRI #hide + +# The first input argument that **KomaMRI** needs for simulating is the phantom. + +# This section goes over the concept of digital phantom and shows how it applies to the specific case +# of **KomaMRI**. We'll go into detail about the [`Phantom`](@ref) structure and its supported operations. + +# ## Digital Phantom + +# A digital phantom is basically a computer model of a physical object (like the human body or a body part) +# which is used in simulations to mimic the characteristics and behaviour that would be obtained from +# real MRI. Instead of using a physical object for testing, the digital phantom allows for virtual experiments. + +# This computer model should essentially contain information about the position and/or displacements +# of the tissues, as well as their MRI-related (T1, T2, PD, off-resonance...) values. + +# ## KomaMRI Phantom Overview +# In Koma, a phantom is made up of a set of spins (which in many cases are also known as ''isochromats''). +# Each spin is independent of the others in terms of properties, position and state. +# This is a key feature of **KomaMRI**, as it is explained in the [Simulation](6-simulation.md) section. + +# Let's take a look at the definition of the [`Phantom`](@ref) struct +# inside Koma's source code to see what it looks like: +# ```julia +# @with_kw mutable struct Phantom{T<:Real} +# name::String = "spins" +# x::AbstractVector{T} +# y::AbstractVector{T} = zeros(eltype(x), size(x)) +# z::AbstractVector{T} = zeros(eltype(x), size(x)) +# ρ::AbstractVector{T} = ones(eltype(x), size(x)) +# T1::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000 +# T2::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000 +# T2s::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000 +# #Off-resonance related +# Δw::AbstractVector{T} = zeros(eltype(x), size(x)) +# #Diffusion +# Dλ1::AbstractVector{T} = zeros(eltype(x), size(x)) +# Dλ2::AbstractVector{T} = zeros(eltype(x), size(x)) +# Dθ::AbstractVector{T} = zeros(eltype(x), size(x)) +# #Motion +# motion::AbstractMotion{T} = NoMotion{eltype(x)}() +# end +# ``` + +# This structure consists of several elements. Most of them are vectors, except for +# the `name` (self-explanatory) and `motion` (explained below) fields. +# These vectors represent object properties, with each element holding a value associated +# with a single magnetization (i.e. a single spin). +# Specifically, `x`, `y` and `z` are the spatial (starting) coordinates of each spin. +# `ρ` stands for the proton density, and `T1`, `T2` and `T2s` (standing for T2*) +# are the well-known relaxation times. `Δw` accounts for off-resonance effects. +# `Dλ1`, `Dλ2` and `Dθ` are diffusion-related fields which are not in use at the moment. +# Last, the `motion` field stands for spin displacements, which are added to `x`, `y` and `z` +# when simulating in order to obtain the spin positions at each time step. For more information about +# motion, refer to [Motion](2-motion.md) section. + +# To get an even better understanding on how it works, let's look at an example of a brain phantom: + +obj = brain_phantom2D() +# ```julia-repl +# Phantom{Float64} +# name: String "brain2D_axial" +# x: Array{Float64}((6506,)) [-0.084, -0.084, -0.084, -0.084, -0.084, -0.084, -0.084, -0.084, -0.084, -0.084 … 0.084, 0.084, 0.084, 0.084, 0.086, 0.086, 0.086, 0.086, 0.086, 0.086] +# y: Array{Float64}((6506,)) [-0.03, -0.028, -0.026, -0.024, -0.022, -0.02, -0.018, -0.016, -0.014, -0.012 … 0.006, 0.008, 0.01, 0.012, -0.008, -0.006, -0.004, -0.002, 0.0, 0.002] +# z: Array{Float64}((6506,)) [-0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0 … 0.0, 0.0, 0.0, 0.0, -0.0, -0.0, -0.0, -0.0, 0.0, 0.0] +# ρ: Array{Float64}((6506,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] +# T1: Array{Float64}((6506,)) [0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569 … 0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569, 0.569] +# T2: Array{Float64}((6506,)) [0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329 … 0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329, 0.329] +# T2s: Array{Float64}((6506,)) [0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058 … 0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058, 0.058] +# Δw: Array{Float64}((6506,)) [-0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0 … -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0] +# Dλ1: Array{Float64}((6506,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +# Dλ2: Array{Float64}((6506,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +# Dθ: Array{Float64}((6506,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +# motion: NoMotion{Float64} NoMotion{Float64}() +# ``` + +# You can visualize the **Phantom** struct using the [`plot_phantom_map`](@ref) function, +# which is part of the **KomaMRIPlots** subdependency. This function plots the magnitude of a property for +# each magnetization at a specific spatial position. You can observe properties such as proton density +# and relaxation times, so feel free to replace the `:T1` symbol with another property of the phantom in the example below: + +p1 = plot_phantom_map(obj, :T1; height=450) + +#md savefig(p1, "../assets/doc-1-phantom.html") #hide +#jl display(p1) + +#md # ```@raw html +#md #