From 7ddb6cfbdca179c131f858854912fc1d451e92d7 Mon Sep 17 00:00:00 2001
From: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
Date: Wed, 31 Jan 2024 11:23:23 +0000
Subject: [PATCH] qval and bval should be in units of rad
---
src/gradients/changing_gradient_blocks.jl | 6 +++---
src/gradients/constant_gradient_blocks.jl | 6 +++---
src/gradients/pulsed_gradients.jl | 2 +-
src/pathways.jl | 8 ++++----
src/variables.jl | 6 +++---
5 files changed, 14 insertions(+), 14 deletions(-)
diff --git a/src/gradients/changing_gradient_blocks.jl b/src/gradients/changing_gradient_blocks.jl
index f631d09..9086770 100644
--- a/src/gradients/changing_gradient_blocks.jl
+++ b/src/gradients/changing_gradient_blocks.jl
@@ -38,7 +38,7 @@ grad_start(cgb::AbstractChangingGradientBlock) = cgb.gradient_strength_start
slew_rate(cgb::AbstractChangingGradientBlock) = cgb.slew_rate
grad_end(cgb::AbstractChangingGradientBlock) = grad_start(cgb) .+ slew_rate(cgb) .* duration(cgb)
gradient_strength(cgb::AbstractChangingGradientBlock) = max.(grad_start(cgb), grad_end(cgb))
-qvec(cgb::AbstractChangingGradientBlock) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb)/2)
+qvec(cgb::AbstractChangingGradientBlock) = (grad_start(cgb) .+ grad_end(cgb)) .* (duration(cgb) * π)
function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)
# grad = (g1 * (duration - t) + g2 * t) / duration
@@ -55,14 +55,14 @@ function bmat_gradient(cgb::AbstractChangingGradientBlock, qstart)
qstart .* permutedims(qstart) .* duration(cgb) .+
duration(cgb)^2 .* qstart .* permutedims(
2 .* grad_start(cgb) .+
- grad_end(cgb)) ./ 3 .+
+ grad_end(cgb)) .* 2Ï€ ./ 3 .+
bmat_gradient(cgb)
)
end
function bmat_gradient(cgb::AbstractChangingGradientBlock)
diff = slew_rate(cgb) .* duration(cgb)
- return (
+ return (2Ï€)^2 .* (
grad_start(cgb) .* permutedims(grad_start(cgb)) ./ 3 .+
grad_start(cgb) .* permutedims(diff) ./ 4 .+
diff .* permutedims(diff) ./ 10
diff --git a/src/gradients/constant_gradient_blocks.jl b/src/gradients/constant_gradient_blocks.jl
index ca3f78c..1f1a8dd 100644
--- a/src/gradients/constant_gradient_blocks.jl
+++ b/src/gradients/constant_gradient_blocks.jl
@@ -35,10 +35,10 @@ const FixedConstantGradientBlock = AbstractConstantGradientBlock{Float64}
duration(cgb::AbstractConstantGradientBlock) = cgb.duration
gradient_strength(cgb::AbstractConstantGradientBlock) = cgb.gradient_strength
slew_rate(::AbstractConstantGradientBlock) = zero(SVector{3, Float64})
-qvec(cgb::AbstractConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb)
+qvec(cgb::AbstractConstantGradientBlock) = duration(cgb) .* gradient_strength(cgb) .* 2Ï€
function bmat_gradient(cgb::AbstractConstantGradientBlock)
- grad = gradient_strength(cgb)
+ grad = 2Ï€ .* gradient_strength(cgb)
return (grad .* permutedims(grad)) .* duration(cgb)^3 ./3
end
@@ -47,7 +47,7 @@ function bmat_gradient(cgb::AbstractConstantGradientBlock, qstart)
# qstart^2 * duration +
# qstart * grad * duration^2 +
# grad * grad * duration^3 / 3 +
- grad = gradient_strength(cgb)
+ grad = 2Ï€ .* gradient_strength(cgb)
return (
qstart .* permutedims(qstart) .* duration(cgb) .+
qstart .* permutedims(grad) .* duration(cgb)^2 .+
diff --git a/src/gradients/pulsed_gradients.jl b/src/gradients/pulsed_gradients.jl
index f081895..bb68539 100644
--- a/src/gradients/pulsed_gradients.jl
+++ b/src/gradients/pulsed_gradients.jl
@@ -107,7 +107,7 @@ gradient_strength(g::PulsedGradient) = gradient_strength(g.flat)
slew_rate(g::PulsedGradient) = g.slew_rate
δ(g::PulsedGradient) = rise_time(g) + flat_time(g)
duration(g::PulsedGradient) = 2 * rise_time(g) + flat_time(g)
-qvec(g::PulsedGradient) = δ(g) .* gradient_strength(g)
+qvec(g::PulsedGradient) = δ(g) .* gradient_strength(g) .* 2π
get_children_indices(::PulsedGradient) = (:rise, :flat, :fall)
Base.getindex(pg::PulsedGradient, symbol::Symbol) = pg[Val(symbol)]
diff --git a/src/pathways.jl b/src/pathways.jl
index f0d0379..4474433 100644
--- a/src/pathways.jl
+++ b/src/pathways.jl
@@ -143,7 +143,7 @@ area_under_curve(pathway::Pathway; scale=nothing, rotate=nothing) = qval(pathway
"""
bmat(pathway::Pathway; scale=nothing, rotate=nothing)
-Return 3x3 diffusion-weighted matrix experienced by the spins following a specific [`Pathway`](@ref).
+Return 3x3 diffusion-weighted matrix experienced by the spins following a specific [`Pathway`](@ref) in rad^2 ms/um^2.
Only gradients active while the spins are in the transverse plane are considered.
@@ -155,12 +155,12 @@ bmat(pathway::Pathway; scale=nothing, rotate=nothing) = get(pathway.bmat, (scal
"""
bval(pathway::Pathway; scale=nothing, rotate=nothing)
-Return size of diffusion-weighting experienced by the spins following a specific [`Pathway`](@ref).
+Return size of diffusion-weighting experienced by the spins following a specific [`Pathway`](@ref) in rad^2 ms/um^2.
-Only gradients active while the spins are in the transverse plane are considered.
+Only gradients active while the spins are in the transverse plane will contribute to the diffusion weighting.
By default gradients that are affected by user-provided `scale` or `rotate` parameters (e.g., bvals/bvecs) are ignored.
-You can set `scale` and/or `rotate` to specific symbols to only consider gradients that are affected by speficic `scale`/`rotate` parameters
+You can set `scale` and/or `rotate` to specific symbols to only consider gradients that are affected by specific `scale`/`rotate` parameters
"""
bval(pathway::Pathway; scale=nothing, rotate=nothing) = tr(bmat(pathway; scale, rotate))
diff --git a/src/variables.jl b/src/variables.jl
index 5d69ac9..0e2d439 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -19,8 +19,8 @@ all_variables_symbols = [
:slice_thickness => (:pulse, "Slice thickness of an RF pulse that is active during a gradient."),
# gradients
- :qvec => (:gradient, "The spatial range and orientation on which the displacements can be detected due to this gradient in 1/um."),
- :qval => (:gradient, "The spatial range on which the displacements can be detected due to this gradient in 1/um."),
+ :qvec => (:gradient, "The spatial range and orientation on which the displacements can be detected due to this gradient in rad/um."),
+ :qval => (:gradient, "The spatial range on which the displacements can be detected due to this gradient in rad/um (i.e., norm of [`qvec`](@ref))."),
:δ => (:gradient, "Effective duration of a gradient pulse ([`rise_time`](@ref) + [`flat_time`](@ref)) in ms."),
:rise_time => (:gradient, "Time for gradient pulse to reach its maximum value in ms."),
:flat_time => (:gradient, "Time of gradient pulse at maximum value in ms."),
@@ -101,7 +101,7 @@ end
"""
bmat_gradient(gradient::GradientBlock, qstart=(0, 0, 0))
-Computes the diffusion-weighting matrix due to a single gradient block.
+Computes the diffusion-weighting matrix due to a single gradient block in rad^2 ms/um^2.
This should be defined for every `GradientBlock`, but not be called directly.
Instead, the `bmat` and `bval` should be constrained for specific `Pathways`
--
GitLab