Common features#
Units#
MRI simulations deal with various quantities: times, frequencies and angular frequencies, magnetic field strength, gradient moments, etc. All those quantities use different usual prefixes, but not always the same: relaxation time and sequence timing are usually expressed in milliseconds, while the gyromagnetic ratio is expressed in MHz/T (i.e. inverse of microseconds) gradients moments are expressed in mT/m⋅ms and slice thickness in millimeters (on clinical scanners) or micrometers (on pre-clinical scanners). This wealth of units makes it very easy to get quantities wrong by a factor of 1000 or more.
Sycomore provides a unit system so that users do not have to convert their quantities to a specific unit. Units may be declared by multiplying or dividing by the unit name (e.g. 500*ms). Those two syntaxes can be mixed in order to use more complex units (e.g. 267.522*MHz/T). Unit objects follow the usual arithmetic rules, and all SI base units, derived units and prefixes are defined.
Quantity (C++, Python) objects contain their value in the base SI unit and may be converted to a compatible unit. Common arithmetic operations (addition, subtraction, multiplication, division, power) are implemented, as well as most of numpy numeric functions. The following code sample summarizes these features.
Note
The micro prefix is u, not μ, in order to keep ASCII names
import sycomore
from sycomore.units import *
# Combination of base units
diffusion_coefficient = 0.89*um**2/ms
# Magnitude of the quantity, in SI unit
length = 180*cm
length_in_meters = length.magnitude # Equals to 1.8
# Conversion: magnitude of the quantity, in specified unit
duration = 1*h
duration_in_seconds = duration.convert_to(s) # Equals to 3600
#include <sycomore/units.h>
int main()
{
using namespace sycomore::units;
// Combination of base units
auto diffusion_coefficient = 0.89*std::pow(um, 2)/ms;
// Magnitude of the quantity, in SI unit
auto length = 180*cm;
auto length_in_meters = length.magnitude; // Equals to 1.8
// Conversion: magnitude of the quantity, in specified unit
auto duration = 1*h;
auto duration_in_seconds = duration.convert_to(s); // Equals to 3600
return 0;
}
Note
Take caution when importing all units name (from sycomore.units import *), as some names may clash with your code. In a long module with many variables, it is better to import only the required units (e.g. from sycomore.units import mT, m, ms)
Species#
A species is characterized by its relaxation rates (R1, R2 and R2’), its diffusivity D and its relative resonance frequency Δω. It is described by the Species (C++, Python) class. R1 and R2 are mandatory parameters, all others are optional and are equal to 0 if unspecified. Using the unit system, relaxation rates and relaxation times may be used interchangeably.
import sycomore
from sycomore.units import *
# Create a Species from either relaxation times, relaxation rates or both
species = sycomore.Species(1000*ms, 100*ms)
species = sycomore.Species(1*Hz, 10*Hz)
species = sycomore.Species(1000*ms, 10*Hz)
#include <sycomore/Species.h>
#include <sycomore/units.h>
int main()
{
using namespace sycomore::units;
// Create a Species from either relaxation times, relaxation rates or both
sycomore::Species species_1(1000*ms, 100*ms);
sycomore::Species species_2(1*Hz, 10*Hz);
sycomore::Species species_3(1000*ms, 10*Hz);
return 0;
}
The diffusivity can be assigned either as a scalar (for isotropic diffusion) or as a tensor (for anisotropic diffusion), but will always be returned as a tensor:
import sycomore
from sycomore.units import *
species = sycomore.Species(1000*ms, 100*ms)
# Assign the diffusion coefficient as a scalar
species.D = 3*um**2/s
# The diffusion coefficient is stored on the diagonal of the tensor
print(species.D[0,0])
# Assign the diffusion coefficient as a tensor
species.D = [
[3*um**2/s, 0*um**2/s, 0*um**2/s],
[0*um**2/s, 2*um**2/s, 0*um**2/s],
[0*um**2/s, 0*um**2/s, 1*um**2/s]]
print(species.D)
3e-12 [ L^2 T^-1 ]
[[3e-12 [ L^2 T^-1 ] 0 [ L^2 T^-1 ] 0 [ L^2 T^-1 ]]
[0 [ L^2 T^-1 ] 2e-12 [ L^2 T^-1 ] 0 [ L^2 T^-1 ]]
[0 [ L^2 T^-1 ] 0 [ L^2 T^-1 ] 1e-12 [ L^2 T^-1 ]]]
#include <iostream>
#include <sycomore/Species.h>
#include <sycomore/units.h>
#include <xtensor/xio.hpp>
int main()
{
using namespace sycomore::units;
sycomore::Species species(1000*ms, 100*ms);
// Assign the diffusion coefficient as a scalar
species.set_D(3*std::pow(um, 2)/s);
// The diffusion coefficient is stored on the diagonal of the tensor
std::cout << species.D()(0,0) << "\n";
// Assign the diffusion coefficient as a tensor
species.set_D( {
{3*std::pow(um, 2)/s, 0*std::pow(um, 2)/s, 0*std::pow(um, 2)/s },
{0*std::pow(um, 2)/s, 2*std::pow(um, 2)/s, 0*std::pow(um, 2)/s },
{0*std::pow(um, 2)/s, 0*std::pow(um, 2)/s, 1*std::pow(um, 2)/s } });
std::cout << species.D() << "\n";
return 0;
}
3e-12 [ L^2 T^-1 ]
{{3e-12 [ L^2 T^-1 ], 0 [ L^2 T^-1 ], 0 [ L^2 T^-1 ]},
{ 0 [ L^2 T^-1 ], 2e-12 [ L^2 T^-1 ], 0 [ L^2 T^-1 ]},
{ 0 [ L^2 T^-1 ], 0 [ L^2 T^-1 ], 1e-12 [ L^2 T^-1 ]}}
Time intervals#
A time interval (C++, Python) is specified by its duration and an optional magnetic field gradient. The gradient can be either as a scalar or as a 3D array, and can describe the amplitude (in T/m), the area (in T/m*s) or the dephasing (in rad/m). A sycomore.TimeInterval.set_gradient() function is available for generic modification of the gradient.
import sycomore
from sycomore.units import *
# Scalar gradient, defined by its amplitude
interval = sycomore.TimeInterval(1*ms, 20*mT/m)
print(
interval.duration,
interval.gradient_amplitude,
interval.gradient_area/interval.duration,
interval.gradient_dephasing/(sycomore.gamma*interval.duration),
sep="\n")
0.001 [ T ]
[0.02 [ L^-1 M T^-2 I^-1 ] 0.02 [ L^-1 M T^-2 I^-1 ]
0.02 [ L^-1 M T^-2 I^-1 ]]
[0.02 [ L^-1 M T^-2 I^-1 ] 0.02 [ L^-1 M T^-2 I^-1 ]
0.02 [ L^-1 M T^-2 I^-1 ]]
[0.02 [ L^-1 M T^-2 I^-1 ] 0.02 [ L^-1 M T^-2 I^-1 ]
0.02 [ L^-1 M T^-2 I^-1 ]]
#include <iostream>
#include <sycomore/sycomore.h>
#include <sycomore/TimeInterval.h>
#include <sycomore/units.h>
#include <xtensor/xio.hpp>
int main()
{
using namespace sycomore::units;
// Scalar gradient, defined by its amplitude
sycomore::TimeInterval interval(1*ms, 20*mT/m);
std::cout
<< interval.duration() << "\n"
<< interval.gradient_amplitude() << "\n"
<< interval.gradient_area()/interval.duration() << "\n"
<< interval.gradient_dephasing()/(sycomore::gamma*interval.duration())
<< "\n";
return 0;
}
0.001 [ T ]
{0.02 [ L^-1 M T^-2 I^-1 ], 0.02 [ L^-1 M T^-2 I^-1 ], 0.02 [ L^-1 M T^-2 I^-1 ]}
{0.02 [ L^-1 M T^-2 I^-1 ], 0.02 [ L^-1 M T^-2 I^-1 ], 0.02 [ L^-1 M T^-2 I^-1 ]}
{0.02 [ L^-1 M T^-2 I^-1 ], 0.02 [ L^-1 M T^-2 I^-1 ], 0.02 [ L^-1 M T^-2 I^-1 ]}