MCMRSimulator.jl public API

Simulation(
    sequences; geometry=[], diffusivity=3.,
    R1=0, T1=Inf, R2=0, T2=Inf, off_resonance=0, MT_fraction=0, permeability=0,, 
    timestep=<default parameters>,
)

Defines the setup of the simulation and stores the output of the run.

Argument

General parameters:

• sequences: One or more sequences. During the spin random walk the simulation will keep track of the spin magnetisations for all of the provided sequences. Sequences can be filenames of pulseq files or KomaMRIFiles.Sequence objects.
• geometry: Set of obstructions, which can be used to restrict the diffusion, produce off-resonance fields, alter the local T1/T2 relaxation, and as sources of magnetisation transfer.
• diffusivity: Rate of the random motion of the spins in um^2/ms (default: 3).
• verbose: set to false to silence descriptions of the simulation parameters (default: true).

MRI properties

These parameters determine the evolution and relaxation of the spin magnetisation.

• R1/T1: sets the longitudinal relaxation rate (R1 in kHz) or relaxation time (T1=1/R1 in ms). This determines how fast the longitudinal magnetisation returns to its equilibrium value of 1.
• R2/T2: sets the transverse relaxation rate (R2 in kHz) or relaxation time (T2=1/R2 in ms). This determines how fast the transverse magnetisation is lost.
• off_resonance: Size of the off-resonance field in this voxel in kHz.

These MRI properties can be overriden for spins inside the ObstructionGroup objects of the geometry.

Collision parameters

These parameters determine how parameters behave when hitting the ObstructionGroup objects of the geometry. They can be overriden for individual objects for each ObstructionGroup.

• MT_fraction: the fraction of magnetisation transfered between the obstruction and the water spin at each collision.
• permeability: the rate of spins passing through the surface in arbitrary units (set to infinity for fully permeable surface).
• surface_density: Density of spins stuck on the surface relative to the volume density of hte free water.
• dwell_time: Typical time that spins will stay at the surface after getting stuck.

Note that MT_fraction and permeability are internally adjusted to make their effect independent of the timestep.

Timestep parameters

timestep controls the timepoints at which the simulation is evaluated. By default, the maximum allowable timestep will be determined by the geometry and biophysical parameters as described in MCMRSimulator.TimeStep. That documentation also describes how to adjust these settings. The timestep can also be set to a number to ignore any of these parameters. Note that a too large timestep will lead to inaccurate results.

Running the simulation

To run a Snapshot of spins through the simulations you can use one of the following functions:

• evolve: evolves the spins in the snapshot until a single given time and returns that state in a new Snapshot.
• readout: evolves the spins to particular times in each TR and return the total signal at that time (or a Snapshot).

readout(spins, simulation[, readout_times]; bounding_box=<1x1x1 mm box>, skip_TR=0, nTR=1, return_snapshot=false, subset=<all>)

Evolves a set of spins through the Simulation. Returns the total signal or a full Snapshot at every readout time in the simulated sequences over one or more repetition times (TRs).

Positional arguments:

• spins: Number of spins to simulate or an already existing Snapshot.
• simulation: Simulation object defining the environment, scanner, and sequence(s).
• times (optional): time of the readouts relative to the start of the TR (in ms). If not provided, the times of any MRIBuilder.ADC objects in the sequence will be used (see get_readouts for details).

Keyword arguments:

• bounding_box: size of the voxel in which the spins are initiated in um (default is 1000, corresponding to a 1x1x1 mm box centered on zero). Can be set to a BoundingBox object for more control.
• skip_TR: Number of repetition times to skip before starting the readout. Even if set to zero (the default), the simulator will still skip the current TR before starting the readout if the starting snapshot is from a time past one of the sequence readouts. See get_readouts for details.
• nTR: number of TRs for which to store the output. See get_readouts for details.
• return_snapshot: set to true to output the state of all the spins as a Snapshot at each readout instead of a SpinOrientationSum with the total signal.
• subset: Return the signal/snapshot for a subset of all spins. Can be set to a single or a vector of Subset objects. If set to a vector, this will add an attional dimension to the output.

Returns

The function returns an up to 3-dimensional (KxLxMxN) array, with the following dimensions:

• K: the number of sequences. This dimension is not included if the simulation only contains a single sequencen (and this single sequence is not passed into the Simulation as a vector).
• L: the number of readout times with a single TR. This dimension is skipped if the readout_times is set to a scalar number. This dimension might contain nothings for sequences that contain fewer Readout.ADC objects than the maximum (M).
• M: the number of TRs (controlled by the nTR keyword). If nTR is not explicitly set by the user, this dimension is skipped.
• N: the number of subsets (controlled by the subset keyword). If subset is set to a single value (<all> by default), this dimension is skipped.

By default each element of this matrix is either a SpinOrientationSum with the total signal. If return_snapshot=true is set, each element is the full Snapshot instead.

evolve(snapshot, simulation[, new_time]; bounding_box=<1x1x1 mm box>)

Evolves the Snapshot through the Simulation to a new time. Returns a Snapshot at the new time, which can be used as a basis for further simulation.

Subset(; bound=nothing, inside=nothing, geometry_index=nothing, obstruction_index=nothing)

This creates a helper object to extract a subset of a Snapshot from the total snapshot. It defines which spins should be kept. This definition is determined by:

• bound: set to true to return only bound spins, to false to return only free spins (default: whether spins are bound is not relevant).
• inside: set to true to return only spins inside the geometry, to false to return only spins outside of the geometry (default: whether spins are inside or outside is not relevant). It can also be set to a positive integer number. Only spins that are inside that exact number of obstructions will be returned.
• geometry_index: set to an index to only consider that obstruction group within the total geometry (default: consider the full geometry).
• obstruction_index: set to an index to only consider that obstruction group within the total geometry (default: consider the full geometry).

Represents the positions and orientations of multiple Spin objects at a specific time.

Note that times are in milliseconds and positions in micrometer. The equilibrium longitudinal spin (after T1 relaxation) is always 1.

Useful constructors

Snapshot(positions; time=0., longitudinal=1., transverse=0., phase=0., nsequences=1)
Snapshot(nspins[, bounding_box[, geometry]]; time=0., longitudinal=1., transverse=0., phase=0., nsequences=1)
Snapshot(nspins, simulation[, bounding_box; time=0., longitudinal=1., transverse=0., phase=0., nsequences=1)

Creates a new Snapshot at the given time with spins initialised for simulating nsequences sequences. All spins will start out in equilibrium, but that can be changed using the longitudinal, transverse, and/or phase flags. This initial spin locations are given by positions (Nx3 matrix or sequence of vectors of size 3). Alternatively the number of spins can be given in which case the spins are randomly distributed in the given bounding_box (default: 1x1x1 mm box centered on origin). The bounding_box can be a BoundingBox object, a tuple with the lower and upper bounds (i.e., two vectors of length 3) or a number r (resulting in spins filling a cube from -r to +r)

Snapshot(snap::Snapshot{1}, nsequences)

Replicates the positions and orientations for a single sequence in the input snapshot across nsequences.

Extracting summary information

• longitudinal(snapshot) to get the nsequences spin magnitudes in the z-direction (equilibrium of 1) averaged over all spins
• transverse(snapshot) to get the nsequences spin magnitudes in the x-y-plane averaged over all spins
• phase(snapshot) to get the nsequences spin angles in x-y plane (in degrees) averaged over all spins
• orientation(snapshot) to get a (nsequencesx3) matrix with the spin orientations in 3D space averaged over all spins
• SpinOrientation(snapshot) to get a nsequences vector of [SpinOrientation] objects with the average spin orientation across all spins
• position.(snapshot) to get a the position for each spin in a vector (no averaging applied)

Information for a single sequence can be extracted by calling get_sequence first.

BoundingBox(lower::Vector, upper::Vector)
BoundingBox([center::Vector, ]radius::Number)

Represents a box in N-dimensional space.

Check whether particles are inside using isinside.

This can be used to define where spins should be generated in a Snapshot (see readout).

Spin particle with a position and nsequences spin orientations (stored as SpinOrientation).

A random number generator is stored in the Spin object as well, which will be used for evolving the spin into the future in a reproducible manner.

Constructors

Spin(;nsequences=1, position=[0, 0, 0], longitudinal=1., transverse=0., phase=0.)

Creates a new spin with nsequences identical spin orientations (given by longitudinal, transverse, and phase flags). The spin will start at given position.

Spin(reference_spin::Spin{1}, nsequences)

Create a new spin with the same position as reference_spin with the orientation of reference_spin replicated nsequences times.

Extracting spin information

• longitudinal to get the nsequences spin magnitudes in the z-direction (equilibrium of 1)
• transverse to get the nsequences spin magnitudes in the x-y-plane
• phase to get the nsequences spin angles in x-y plane (in degrees)
• orientation to get a (nsequencesx3) matrix with the spin orientations in 3D space
• position to get a length-3 vector with spin location

get_subset(snapshot, simulation, subset)
get_subset(snapshot, simulation; bound=nothing, inside=nothing, geometry_index=nothing, obstruction_index=nothing)

Returns a subset of the Snapshot from the Simulation that obey some specific properties. These properties can be either defined by a Subset object or a set of keyword arguments.

These keyword arguments are:

• bound: set to true to return only bound spins, to false to return only free spins (default: whether spins are bound is not relevant).
• inside: set to true to return only spins inside the geometry, to false to return only spins outside of the geometry (default: whether spins are inside or outside is not relevant). It can also be set to a positive integer number. Only spins that are inside that exact number of obstructions will be returned.
• geometry_index: set to an index to only consider that obstruction group within the total geometry (default: consider the full geometry).
• obstruction_index: set to an index to only consider that obstruction group within the total geometry (default: consider the full geometry).

SpinOrientation(longitudinal, transverse, phase)

The spin orientation. Usually created as part of a Spin object.

SpinOrientation(snapshot::Snapshot)

Returns the average spin orientations of all Spin objects in the Snapshot.

This information can be extracted using:

• longitudinal to get the spin in the z-direction (equilibrium of 1)
• transverse to get the spin in the x-y-plane
• phase to get the spin angle in x-y plane (in degrees)
• orientation to get the spin orientation as a length-3 vector

SpinOrientationSum(snapshot)

Computes the total signal and the number of spins in a Snapshot. The number of spins can be found by running length(spin_orientation_sum). The spin orientation information can be found in the same way as for SpinOrientation, namely by calling transverse, longitudinal, or phase.

Simulation(timestep=(tortuosity=3e-2, gradient=1e-4, permeability=0.5, surface_relaxation=0.01, transfer_rate=0.01, dwell_time=0.1))

Creates an object controlling the timestep of the MCMR simulation.

It can be set by supplying a named tuple to the timestep keyword when creating a Simulation.

At any time the timestep is guaranteed to be shorter than:

1. tortuosity * size_scale(geometry)^2 / D, where size_scale is the average size of the obstructions and D is the diffusivity.
2. timestep greater than permeability times 1 / (maximum permeability parameter)^2
3. timestep that would allow surface relaxation rate at single collision to be greater than surface_relaxation.
4. timestep that would allow magnetisation transfer rate at single collision to be greater than transfer_rate.
5. the minimum dwell time of the bound pool times dwell_time.
6. (gradient /( D * \gamma^2 * G^2))^(1//3), where \gamma is the gyromagnetic_ratio and G is the current gradient_strength.

MCMRSimulator.Geometries.User.Obstructions.ObstructionGroups.ObstructionGroup{:Annuli})(; fields...)

Creates a set of MCMRSimulator.Annulus objects. Fields can be set using keyword arguments. The following fields are available:

• inner: Radius of the inner cylinder. Field is required. Expected type: Float64.
• outer: Radius of the outer cylinder. Field is required. Expected type: Float64.
• myelin: Whether the annulus is myelinated. Field is required. Expected type: Bool. default value: false
• rotation: Rotation applied to all obstructions in group. Can be set to a matrix or one of :x, :y, or, :z (see MCMRSimulator.get_rotation). Field is required. Expected type: StaticArraysCore.SMatrix{3, 2, Float64, 6}. default value: [1.0 0.0; 0.0 1.0; 0.0 0.0]
• R1_inner_surface: Additional longitudinal relaxation rate (kHz). Inner_surface property. Field is required. Expected type: Float64. default value: 0.0
• R1_outer_surface: Additional longitudinal relaxation rate (kHz). Outer_surface property. Field is required. Expected type: Float64. default value: 0.0
• R1_inner_volume: Additional longitudinal relaxation rate (kHz). Inner_volume property. Field is required. Expected type: Float64. default value: 0.0
• R1_outer_volume: Additional longitudinal relaxation rate (kHz). Outer_volume property. Field is required. Expected type: Float64. default value: 0.0
• R2_inner_surface: Additional transverse relaxation rate (kHz). Inner_surface property. Field is required. Expected type: Float64. default value: 0.0
• R2_outer_surface: Additional transverse relaxation rate (kHz). Outer_surface property. Field is required. Expected type: Float64. default value: 0.0
• R2_inner_volume: Additional transverse relaxation rate (kHz). Inner_volume property. Field is required. Expected type: Float64. default value: 0.0
• R2_outer_volume: Additional transverse relaxation rate (kHz). Outer_volume property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_inner_surface: Additional off-resonance field offset (kHz). Inner_surface property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_outer_surface: Additional off-resonance field offset (kHz). Outer_surface property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_inner_volume: Additional off-resonance field offset (kHz). Inner_volume property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_outer_volume: Additional off-resonance field offset (kHz). Outer_volume property. Field is required. Expected type: Float64. default value: 0.0
• position: Spatial offset of obstruction from origin. Field is required. Expected type: StaticArraysCore.MVector{2, Float64}. default value: [0.0, 0.0]
• susceptibility_iso: Isotropic component of the myelin susceptibility (in ppm). Field can be null. Expected type: Float64. default value: -0.1
• susceptibility_aniso: Anisotropic component of the myelin susceptibility (in ppm). Field can be null. Expected type: Float64. default value: -0.1
• lorentz_radius: Only compute field explicitly for a annuli with this Lorentz radius. Field can be null. Expected type: Float64. default value: 5.0
• repeats: Length scale on which the obstructions are repeated (um). Field can be null. Expected type: StaticArraysCore.MVector{2, Float64}.
• use_boundingbox: Use bounding boxes for an initial filtering of possible intersections. Field can be null. Expected type: Bool. default value: true
• grid_resolution: Resolution of the grid that the volume is split up into (um). Defaults to roughly one grid element per obstruction. Field can be null. Expected type: Float64.
• dwell_time_inner_surface: Average time a particle stays stuck to the surface (ms). Inner_surface property. Field can be null. Expected type: Float64.
• dwell_time_outer_surface: Average time a particle stays stuck to the surface (ms). Outer_surface property. Field can be null. Expected type: Float64.
• density_inner_surface: Surface density of stuck particles relative to the volume density (um). Inner_surface property. Field can be null. Expected type: Float64.
• density_outer_surface: Surface density of stuck particles relative to the volume density (um). Outer_surface property. Field can be null. Expected type: Float64.
• permeability_inner_surface: Rate of particle passing through the obstruction in arbitrary units. Inner_surface property. Field can be null. Expected type: Float64.
• permeability_outer_surface: Rate of particle passing through the obstruction in arbitrary units. Outer_surface property. Field can be null. Expected type: Float64.
• relaxation_inner_surface: Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Inner_surface property. Field can be null. Expected type: Float64.
• relaxation_outer_surface: Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Outer_surface property. Field can be null. Expected type: Float64.
• size_scale: Size of the smallest obstructions. If not set explicitly, this will be determined by the minimum radius or distance between objects (see size_scale). Field can be null. Expected type: Float64.

Singular Annulus object that is obtained by indexing a Annuli object.

MCMRSimulator.Geometries.User.Obstructions.ObstructionGroups.ObstructionGroup{:Cylinders})(; fields...)

Creates a set of MCMRSimulator.Cylinder objects. Fields can be set using keyword arguments. The following fields are available:

• radius: Radius of the cylinder. Field is required. Expected type: Float64.
• rotation: Rotation applied to all obstructions in group. Can be set to a matrix or one of :x, :y, or, :z (see MCMRSimulator.get_rotation). Field is required. Expected type: StaticArraysCore.SMatrix{3, 2, Float64, 6}. default value: [1.0 0.0; 0.0 1.0; 0.0 0.0]
• R1_surface: Additional longitudinal relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R1_inside: Additional longitudinal relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• R2_surface: Additional transverse relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R2_inside: Additional transverse relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_surface: Additional off-resonance field offset (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_inside: Additional off-resonance field offset (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• position: Spatial offset of obstruction from origin. Field is required. Expected type: StaticArraysCore.MVector{2, Float64}. default value: [0.0, 0.0]
• g_ratio: Inner/outer radius used for susceptibility calculation Field can be null. Expected type: Float64. default value: 1.0
• susceptibility_iso: Isotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the isotropic susceptibility of the simulated tissue by the thickness. Field can be null. Expected type: Float64. default value: -0.1
• susceptibility_aniso: Anisotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the anisotropic susceptibility of the simulated tissue by the thickness Field can be null. Expected type: Float64. default value: -0.1
• lorentz_radius: Only compute field explicitly for cylinders with this Lorentz radius. Field can be null. Expected type: Float64. default value: 5.0
• repeats: Length scale on which the obstructions are repeated (um). Field can be null. Expected type: StaticArraysCore.MVector{2, Float64}.
• use_boundingbox: Use bounding boxes for an initial filtering of possible intersections. Field can be null. Expected type: Bool. default value: true
• grid_resolution: Resolution of the grid that the volume is split up into (um). Defaults to roughly one grid element per obstruction. Field can be null. Expected type: Float64.
• dwell_time: Average time a particle stays stuck to the surface (ms). Surface property. Field can be null. Expected type: Float64.
• density: Surface density of stuck particles relative to the volume density (um). Surface property. Field can be null. Expected type: Float64.
• permeability: Rate of particle passing through the obstruction in arbitrary units. Surface property. Field can be null. Expected type: Float64.
• relaxation: Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Surface property. Field can be null. Expected type: Float64.
• size_scale: Size of the smallest obstructions. If not set explicitly, this will be determined by the minimum radius or distance between objects (see size_scale). Field can be null. Expected type: Float64.

Singular Cylinder object that is obtained by indexing a Cylinders object.

MCMRSimulator.Geometries.User.Obstructions.ObstructionGroups.ObstructionGroup{:Walls})(; fields...)

Creates a set of MCMRSimulator.Wall objects. Fields can be set using keyword arguments. The following fields are available:

• rotation: Rotation applied to all obstructions in group. Can be set to a matrix or one of :x, :y, or, :z (see MCMRSimulator.get_rotation). Field is required. Expected type: StaticArraysCore.SMatrix{3, 1, Float64, 3}. default value: [1.0; 0.0; 0.0;;]
• R1: Additional longitudinal relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R2: Additional transverse relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance: Additional off-resonance field offset (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• position: Spatial offset of obstruction from origin. Field is required. Expected type: Float64. default value: 0.0
• repeats: Length scale on which the obstructions are repeated (um). Field can be null. Expected type: StaticArraysCore.MVector{1, Float64}.
• use_boundingbox: Use bounding boxes for an initial filtering of possible intersections. Field can be null. Expected type: Bool. default value: true
• grid_resolution: Resolution of the grid that the volume is split up into (um). Defaults to roughly one grid element per obstruction. Field can be null. Expected type: Float64.
• dwell_time: Average time a particle stays stuck to the surface (ms). Surface property. Field can be null. Expected type: Float64.
• density: Surface density of stuck particles relative to the volume density (um). Surface property. Field can be null. Expected type: Float64.
• permeability: Rate of particle passing through the obstruction in arbitrary units. Surface property. Field can be null. Expected type: Float64.
• relaxation: Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Surface property. Field can be null. Expected type: Float64.
• size_scale: Size of the smallest obstructions. If not set explicitly, this will be determined by the minimum radius or distance between objects (see size_scale). Field can be null. Expected type: Float64.

Singular Wall object that is obtained by indexing a Walls object.

MCMRSimulator.Geometries.User.Obstructions.ObstructionGroups.ObstructionGroup{:Spheres})(; fields...)

Creates a set of MCMRSimulator.Sphere objects. Fields can be set using keyword arguments. The following fields are available:

• radius: Radius of the sphere. Field is required. Expected type: Float64.
• rotation: Rotation applied to all obstructions in group. Can be set to a matrix or one of :x, :y, or, :z (see MCMRSimulator.get_rotation). Field is required. Expected type: StaticArraysCore.SMatrix{3, 3, Float64, 9}. default value: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
• R1_surface: Additional longitudinal relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R1_inside: Additional longitudinal relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• R2_surface: Additional transverse relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R2_inside: Additional transverse relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_surface: Additional off-resonance field offset (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_inside: Additional off-resonance field offset (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• position: Spatial offset of obstruction from origin. Field is required. Expected type: StaticArraysCore.MVector{3, Float64}. default value: [0.0, 0.0, 0.0]
• overlapping: Whether overlapping spheres should be treated as permeable in the overlapping region. Field can be null. Expected type: Bool. default value: false
• repeats: Length scale on which the obstructions are repeated (um). Field can be null. Expected type: StaticArraysCore.MVector{3, Float64}.
• use_boundingbox: Use bounding boxes for an initial filtering of possible intersections. Field can be null. Expected type: Bool. default value: true
• grid_resolution: Resolution of the grid that the volume is split up into (um). Defaults to roughly one grid element per obstruction. Field can be null. Expected type: Float64.
• dwell_time: Average time a particle stays stuck to the surface (ms). Surface property. Field can be null. Expected type: Float64.
• density: Surface density of stuck particles relative to the volume density (um). Surface property. Field can be null. Expected type: Float64.
• permeability: Rate of particle passing through the obstruction in arbitrary units. Surface property. Field can be null. Expected type: Float64.
• relaxation: Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Surface property. Field can be null. Expected type: Float64.
• size_scale: Size of the smallest obstructions. If not set explicitly, this will be determined by the minimum radius or distance between objects (see size_scale). Field can be null. Expected type: Float64.

Singular Sphere object that is obtained by indexing a Spheres object.

MCMRSimulator.Geometries.User.Obstructions.ObstructionGroups.ObstructionGroup{:Mesh})(; fields...)

Creates a set of MCMRSimulator.Triangle objects. Fields can be set using keyword arguments. The following fields are available:

• triangles: Each triangle is defined by 3 vertices into the mesh. Field is required. Expected type: StaticArraysCore.MVector{3, Int64}.
• vertices: Positions of the corners of the triangular mesh. Field is required. Expected type: Vector{StaticArraysCore.MVector{3, Float64}}.
• myelin: Whether the mesh is myelinated. Field is required. Expected type: Bool. default value: false
• rotation: Rotation applied to all obstructions in group. Can be set to a matrix or one of :x, :y, or, :z (see MCMRSimulator.get_rotation). Field is required. Expected type: StaticArraysCore.SMatrix{3, 3, Float64, 9}. default value: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
• R1_surface: Additional longitudinal relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R1_inside: Additional longitudinal relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• R2_surface: Additional transverse relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R2_inside: Additional transverse relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_surface: Additional off-resonance field offset (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_inside: Additional off-resonance field offset (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• susceptibility_iso: Isotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the isotropic susceptibility of the simulated tissue by the thickness. Field can be null. Expected type: Float64. default value: -0.1
• susceptibility_aniso: Anisotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the anisotropic susceptibility of the simulated tissue by the thickness. Field can be null. Expected type: Float64. default value: -0.1
• lorentz_radius: Only compute field explicitly for triangles with this Lorentz radius. Field can be null. Expected type: Float64. default value: 5.0
• components: Which component this triangle belongs to. If not provided explicitly, this will be determined based on connectivity. Field can be null. Expected type: Int64.
• repeats: Length scale on which the obstructions are repeated (um). Field can be null. Expected type: StaticArraysCore.MVector{3, Float64}.
• use_boundingbox: Use bounding boxes for an initial filtering of possible intersections. Field can be null. Expected type: Bool. default value: true
• grid_resolution: Resolution of the grid that the volume is split up into (um). Defaults to roughly one grid element per obstruction. Field can be null. Expected type: Float64.
• dwell_time: Average time a particle stays stuck to the surface (ms). Surface property. Field can be null. Expected type: Float64.
• density: Surface density of stuck particles relative to the volume density (um). Surface property. Field can be null. Expected type: Float64.
• permeability: Rate of particle passing through the obstruction in arbitrary units. Surface property. Field can be null. Expected type: Float64.
• relaxation: Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Surface property. Field can be null. Expected type: Float64.
• size_scale: Size of the smallest obstructions. If not set explicitly, this will be determined by the minimum radius or distance between objects (see size_scale). Field can be null. Expected type: Float64.

Singular Triangle object that is obtained by indexing a Mesh object.

MCMRSimulator.Geometries.User.Obstructions.ObstructionGroups.ObstructionGroup{:BendyCylinder})(; fields...)

Creates a set of MCMRSimulator.Ring objects. Fields can be set using keyword arguments. The following fields are available:

• control_point: Control points defining the path of the cylinder. Field is required. Expected type: StaticArraysCore.MVector{3, Float64}.
• radius: Radius at each control point. Field is required. Expected type: Float64.
• spline_order: Sets the order of the b-spine interpolating between control points. Field is required. Expected type: Int64. default value: 3
• myelin: Whether the cylinder is myelinated. Field is required. Expected type: Bool. default value: false
• rotation: Rotation applied to all obstructions in group. Can be set to a matrix or one of :x, :y, or, :z (see MCMRSimulator.get_rotation). Field is required. Expected type: StaticArraysCore.SMatrix{3, 3, Float64, 9}. default value: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
• R1_surface: Additional longitudinal relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R1_inside: Additional longitudinal relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• R2_surface: Additional transverse relaxation rate (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• R2_inside: Additional transverse relaxation rate (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_surface: Additional off-resonance field offset (kHz). Surface property. Field is required. Expected type: Float64. default value: 0.0
• off_resonance_inside: Additional off-resonance field offset (kHz). Inside property. Field is required. Expected type: Float64. default value: 0.0
• nsamples: Number of mesh vertices along each ring. Field can be null. Expected type: Int64. default value: 100
• closed: After how many repeats in each dimension does the cylinder connect with itself. If not set the cylinder is not closed. Field can be null. Expected type: StaticArraysCore.MVector{3, Int64}.
• susceptibility_iso: Isotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the isotropic susceptibility of the simulated tissue by the thickness. Field can be null. Expected type: Float64. default value: -0.1
• susceptibility_aniso: Anisotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the anisotropic susceptibility of the simulated tissue by the thickness. Field can be null. Expected type: Float64. default value: -0.1
• lorentz_radius: Only compute field explicitly for triangles with this Lorentz radius. Field can be null. Expected type: Float64. default value: 5.0
• repeats: Length scale on which the obstructions are repeated (um). Field can be null. Expected type: StaticArraysCore.MVector{3, Float64}.
• use_boundingbox: Use bounding boxes for an initial filtering of possible intersections. Field can be null. Expected type: Bool. default value: true
• grid_resolution: Resolution of the grid that the volume is split up into (um). Defaults to roughly one grid element per obstruction. Field can be null. Expected type: Float64.
• dwell_time: Average time a particle stays stuck to the surface (ms). Surface property. Field can be null. Expected type: Float64.
• density: Surface density of stuck particles relative to the volume density (um). Surface property. Field can be null. Expected type: Float64.
• permeability: Rate of particle passing through the obstruction in arbitrary units. Surface property. Field can be null. Expected type: Float64.
• relaxation: Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Surface property. Field can be null. Expected type: Float64.
• size_scale: Size of the smallest obstructions. If not set explicitly, this will be determined by the minimum radius or distance between objects (see size_scale). Field can be null. Expected type: Float64.

Singular Ring object that is obtained by indexing a BendyCylinder object.

load_mesh(file)

Loads a Mesh from a file.

Currently only PLY files are supported.

PLY format

PLY stands for Polygon File Format (http://paulbourke.net/dataformats/ply/). PLY IO is handled by PlyIO.jl (https://github.com/JuliaGeometry/PlyIO.jl).

read_geometry(io::IO; format=:json, kwargs...)

Read geometry from an open stream. Since streams do not have a filename extension, format must be provided for non-JSON input.

read_geometry(filename::AbstractString; format=nothing, kwargs...)

Read geometry from a JSON, PLY, or SWC file. When format is omitted, the format is inferred from the filename extension.

read_geometry_json(filename/io)

Read geometry from a JSON file.

read_swc(swc_file; swc_as_spheres=false, kwargs...)

Read an SWC file and return a Spheres object with the node positions and radii.

See Spheres for the available keyword arguments.

In the future, this function may be extended to return a more complete representation of the SWC file, including the links between spheres. For now, it only returns the spheres themselves with the overlapping keyword set to true by default.

For forwards compatibility, the swc_as_spheres keyword argument is provided, but it must be set to true to avoid an error. In the future, this argument may be removed when the function is extended to return a more complete representation of the SWC file.

random_positions_radii(box_size, target_density, n_dimensions; distribution=Gamma, mean=1., variance=1., max_iter=1000, min_radius=0.1, max_radius=Inf)

Randomly distributes circles or spheres in space.

Arguments:

• box_size: Size of the infinitely repeating box of random positions
• target_density: Final density of the circles/spheres. This density will only be approximately reached
• n_dimensions: dimensionality of the space (2 for cicles; 3 for spheres)
• distribution: distribution from which the radii are drawn (from Distributions.jl)
• mean: mean of the gamma distribution (ignored if distribution explicitly set)
• variance: variance of the gamma distribution (ignored if distribution explicitly set)
• max_iter: maximum number of iterations to try to prevent the circles/spheres from overlapping. An error is raised if they still overlap after this number of iterations.
• repulsion_strength: strength of the repulsion in each step (defaults to 0.001).
• max_iter_repulse: maximum number of iterations that circles/spheres will be repulsed from each other
• min_radius: samples from the distribution below this radius will be rejected (in um).
• max_radius: samples from the distribution above this radius will be rejected (in um).

read_pulseq(IO)

Reads a sequence from a pulseq file (http://pulseq.github.io/). Pulseq files can be produced using matlab (http://pulseq.github.io/) or python (https://pypulseq.readthedocs.io/en/master/).

write_pulseq(IO, sequence)

Writes a sequence to an output IO file.

Scanner(;B0=3., gradient=Inf, slew_rate=Inf, units=:kHz)

Properties of an MRI scanner relevant for the MR signal simulations.

• B0: magnetic field strength (in Tesla)
• gradient_strength: maximum gradient strength long each axis.
• slew_rate: maximum rate of change in the gradient strength

By default gradient and slew_rate are expected to be provided in units of, respectively, kHz/um and kHz/um/ms. However, if the keyword units=:Tesla is set, the gradient and slew_rate should be provided in units of, respectively, mT/m and T/m/s.

A default 1.5T scanner.

Matches the one used in pulseq (https://github.com/pulseq/pulseq/blob/master/matlab/%2Bmr/opts.m).

Siemens MAGNETOM 3T Prisma MRI scanner (https://www.siemens-healthineers.com/en-uk/magnetic-resonance-imaging/3t-mri-scanner/magnetom-prisma).

Siemens MAGNETOM 7T Terra MRI scanner (https://www.siemens-healthineers.com/en-uk/magnetic-resonance-imaging/7t-mri-scanner/magnetom-terra)

Siemens 3T Connectom MRI scanner ((Fan et al., 2022)).

Siemens 3T Connectom 2.0 MRI scanner ((Huang et al., 2021))

position(s::Spin)

Returns the position of the spin particle as a vector of length 3.

position.(s::Snapshot)

Returns all the positions of the spin particles as a vector of length-3 vectors.

longitudinal(spin)
longitudinal(snapshot)

Returns the longitudinal magnitude of the spin (i.e., magnitude aligned with the magnetic field) for a single particle (Spin) or averaged across a group of particles in a [Snapshot]. When orientations for multiple sequences are available an array of longitudinal values is returned with a value for each sequence.

transverse(spin)
transverse(snapshot)

Returns the transverse spin (i.e., magnitude in the plane perpendicular to the magnetic field) for a single particle (Spin) or averaged across a group of particles in a [Snapshot]. When orientations for multiple sequences are available an array of transverse values is returned with a value for each sequence.

phase(spin)
phase(snapshot)

Returns the phase in the x-y plane of the spin for a single particle (Spin) or averaged across a group of particles in a [Snapshot]. When orientations for multiple sequences are available an array of phase values is returned with a value for each sequence.

orientation(spin)
orientation(snapshot)

Returns the spin orientation as a length-3 vector for a single particle (Spin) or averaged across a group of particles in a [Snapshot]. When orientations for multiple sequences are available an array of vectors is returned with a value for each sequence.

get_sequence(spin, sequence_index)
get_sequence(snapshot, sequence_index)

Extracts the spin orientation corresponding to a specific sequence, where the sequence index uses the order in which the sequences where provided in the Simulation.

isinside(obstruction, position)

Returns true if the position is inside the obstruction. This will be inaccurate if the particle with that position is stuck on the surface of the obstruction.

isinside([geometry, ]spin)

Returns vector of obstructions that the spin is inside. If geometry is not provided, will return a vector of indices instead. If a non-fixed geometry is provided, will return the number of obstructions that the spin is inside.

isinside(grid, position)

Get the indices of obstructions that contain the position.

isinside(obstruction_group, position[, stuck_to])

Returns a vector of indices with all the obstructions in FixedObstructionGroup containing the position (in order). For obstructions with only a single inside, will return an empty vector ("[]") if the particle is outside and a "[0]" if inside.

isinside(geometry, position[, stuck_reflection])

Returns a vector of pairs of indices with all the obstructions in FixedGeometry containing the position (in order). The first index indicates the index of the FixedObstructionGroup within the geometry. The second index is the index of the specific FixedObstruction within the obstruction group.

stuck(spin)

Returns true if the spin is stuck on the surface. This can be used to filter a Snapshot using:

only_stuck = filter(stuck, snapshot)
only_free = filter(s -> !stuck(s), snapshot)

stuck_to(spin)

Return the indices of the obstruction the spin is stuck to. Will return (0, 0) for a free spin.

stuck_to(spin, geometry)

Return the internal representation of the obstruction the spin is stuck to. Raises an error if the spin is free.

get_time(snapshot)
get_time(sequence_component)
get_time(sequence, sequence_index)

Returns the time in milliseconds that a snapshot was taken or that a sequence component will have effect.

Defines a finite plane in the 3D space used for plotting.

Constructor

PlotPlane(normal=[0, 0, 1], position=[0, 0, 0]; size=10., sizex=<size>, sizey=<size>, ngrid=100)

Arguments:

• normal: length-3 vector with the orientation perpendicular to the plane (default to z-direction).
• position: position of plane as a length-3 vector (defaults to origin).
• sizex: size of the plane in the x-direction (before rotating to normal).
• sizey: size of the plane in the y-direction (before rotating to normal).
• size: set default value for sizex and sizey.
• ngrid: number of grid elements to split the plane up into for plotting.

plot([plot_plane], snapshot; kind=:scatter, kwargs...)
plot!([scene,] [plot_plane], snapshot; kind=:scatter, kwargs...)
plot_snapshot([plot_plane], snapshot; kind=:scatter, kwargs...)
plot_snapshot!([scene,] [plot_plane], snapshot; kind=:scatter, kwargs...)

Plots a Snapshot in a new plot.

The spins are plotted in 2D projected onto the PlotPlane if one is provided. Otherwise, the spins are plotted in 3D (does not work for kind=:image).

There are three kinds of snapshot plots available:

Scatter plot

Default (or set using kind=:scatter). Each spin is plotted as a point with the colour set by the transverse magnetisation. Additional keywords are passed on to Makie.scatter (namely, marker, markersize, strokecolor, strokewidth, glowcolor, glowwidth).

Dyad plot

Set using kind=:dyad. Each spin is plotted as a dyad. For a 2D dyad the orientation is set by the transverse magnetisation. For a 3D dyad the full magnetisation is used to set the orientation. Additional keywords are passed on to Makie.arrows (namely, arrowsize, arrowhead, arrowtail, linestyle, lengthscale, quality, markerspace, diffuse, specular, shininess).

Image

Set using kind=:image. The average magnetisation is plotted across the PlotPlane. The colour in each pixel is set by the average transverse magnetisation of the local spins. Set interpolate=false to disable interpolating the colors. This plot will not work in 3D (i.e., a PlotPlane is required).

This function will only work if a Makie backend is imported.

plot([plot_plane,] geometry; kwargs...)
plot!([scene,] [plot_plane,] geometry; kwargs...)
plot_geometry([plot_plane,] geometry; kwargs...)
plot_geometry!([scene,] [plot_plane,] geometry; kwargs...)

Plots a given geometry. If a PlotPlane is provided the 2D projection of the geometry onto this plane is plotted. Otherwise, the geometry is plotted in 3D.

If you want to overlay the off-resonance field, call plot_off_resonance first before calling this function.

This function will only work if a Makie backend is imported.

plot([plot_plane], snapshots; kwargs...)
plot!([scene,] [plot_plane], snapshots; kwargs...)
plot_trajectory([plot_plane], snapshots; kwargs...)
plot_trajectory!([scene,] [plot_plane], snapshots; kwargs...)

Plots the spin trajectory in a vector of Snapshot on an existing plot scene.

The spins are plotted in 2D projected onto the PlotPlane if one is provided. Otherwise, the spins are plotted in 3D. At each location along the trajectory, the colour is set by the transverse magnetisation. Additional keywords are passed on to Makie.lines!.

This function will only work if a Makie backend is imported.

plot_off_resonance(plot_plane, geometry; kwargs...)
plot_off_resonance!([scene,] plot_plane, geometry; kwargs...)

Plots the off-resonance field generated by a given geometry within the PlotPlane.

This function will only work if a Makie backend is imported.

simulator_movie(filename, simulation, times, size; resolution=(1600, 800), trajectory_init=30, signal_init=10000, framerate=50, plane_orientation=:z, kwargs...)

Writes a movie of the Simulation to the given filename.

Each frame of the movie shows the Snapshot at given times. size is a tuple with the size of the plotted region in the x- and y-direction. If there is a repeating geometry, then it is strongly recommended to use the size of the repeat for size.

Keyword arguments:

• resolution: pixel resolution of each frame in the movie.
• trajectory_init: how many spins to plot on each frame.
• signal_init: how many spins to use to evaluate the signal evolution.
• framerate: wait time between each subsequent frame in the movie.
• plane_orientation: orienation of the plane on which the spins are projected (see PlotPlane).

Additional keyword arguments are passed on to plot_snapshot.

This function will only work if Makie is installed and imported.

gyromagnetic ratio of a proton (1H) in water in kHz/T

get_rotation(rotation_mat, ndim)

Returns a (3, ndim) rotation matrix, that is the relevant part of the full 3x3 rotation_mat to map to the x-axis (if ndim is 1) or the x-y plance (if ndim is 2). If ndim is 3, the full rotation matrix rotation_mat is returned.

get_rotation(vector, ndim; reference_dimension)

Returns the (3, ndim) rotation matrix mapping the vector to the reference_dimension. By default, the reference_dimension is the x-direction (ifndimis 1 or 3) or the z-direction (ifndimis 2).vectorandreference_dimension` can be a length 3 array or one of the symbols :x, :y, or :z (representing vectors in those cardinal directions).

get_readouts(sequence[, start_time]; readouts=nothing, nTR=1, skip_TR=0)
get_readouts(adc_sample_times, TR[, start_time]; readouts=nothing, nTR=1, skip_TR=0)

Returns whether this sequence repeats and when it will be readout.

The return in a tuple with:

1. a boolean indicating whether the sequence is repeating or not, and
2. a vector of IndexedReadouts indicating the readouts that will be used in the simulation.

The sequence adc_sample_times and TR can be given directly, or they can be extracted from the sequence using sequencee_waveform(sequence, start_time; kwargs...).

This can be used to identify which readouts will be (or have been) used in the simulation by running: collect(get_readouts(sequence, snapshot.current_time; kwargs...)) where snapshot is the starting snapshot (which has a current_time of 0 by default) and kwargs are the keyword arguments used in readout (i.e., readouts,nTR, andskip_TR`).

By default the readout/ADC objects with the actual sequence definition are used. These can be overriden by readouts, which can be set to a vector of the timings of the readouts within each TR.

Non-repeating sequences

This is the behaviour if both nTR and skip_TR are not set by the user. Any readouts before the start_time are ignored. Any readouts after the start_time (whether from the readouts keyword or within the sequence definition) are returned.

Repeating sequences

This is the behaviour if either nTR or skip_TR or both are set by the user. If at start_time any of the readouts in the current TR have already passed, then the readout will only start in the next TR (unless skip_TR is set to -1). We will skip an additional number of TRs given by skip_TR (default: 0). Then readouts will continue for the number of TRs given by nTR (default: 1).

IndexedReadout(time, TR, readout)

Represents a readout in an MR sequence.

Generate these for a sequence of interest using get_readouts.

All indices are integers. They refer to:

• time: time since beginning of this repetition of the sequence in ms.
• TR: which TR the simulation is in (defaults to 0 if not set).
• readout: which readout within the TR (or total sequence) this is.

susceptibility_off_resonance(susceptibility_grid, position[, inside])

Computes the susceptibility off-resonance caused by a SusceptibilityGrid at given position.

The field is computed in ppm. Knowledge of the scanner B0 is needed to convert it into KHz.

susceptibility_off_resonance(simulation, spin)

Computes the susceptibility off-resonance caused by all susceptibility sources in the Simulation affecting the Spin

The field is computed in ppm. Knowledge of the scanner B0 is needed to convert it into KHz.

Parent type of all geometry objects like Cylinders or Mesh.

A single instance of an ObstructionGroup obtained by indexing.

nvolumes(group)

For any ObstructionType with group_volumes, this computes the number of distinct volumes.

This is mainly used to compute the number of distinct components in a mesh.

fix(user_geometry; permeability=0., density=0., dwell_time=0., relaxation=0.)

Creates a fixed version of the user-created geometry that will be used internally by the simulator.

fix_susceptibility(geometry)

Create a new FixedSusceptibility based on the user-defined geometry settings.

FixedGeometry([obstruction_groups...])

A collection of FixedObstructionGroup objects each reperesenting part of the geometry.

A fixed version of ObstructionGroup that is used internally within the simulator.

This is the main internal representation of a group of identical FixedObstruction objects.

Properties:

• obstructions: vector of the actual FixedObstruction objects.
• parent_index: Index of this group within the larger FixedGeometry.
• original_index: Index of this group within the original user-provided geometry.
• rotation: rotation from global 3-dimensional space to the 1, 2, or 3-dimensional space of the obstructions.
• inv_rotation: inverse of the rotation above
• grid: Grid object on which the obstruction intersections have been precomputed. This speeds up the detection of intersections.
• bounding_boxes: vector of BoundingBox objects for each obstruction. These are used to predect whether a spin could intersect with the obstruction.
• volume: R1, R2, and off-resonance properties of the spins inside the obstructions.
• surface: R1, R2, off-resonance, surfacedensity and dwelltime properties of particles stuck to the surface. Also, contains the permeability and surface relaxation to process collsions.
• vertices: vector of vertices (only used for a mesh).

surface_relaxation(geometry, reflection)

Returns the surface_relaxation experienced by the spin hitting the surface represented by a Reflection. The geometry has to be a FixedGeometry.

surface_density(geometry, reflection)

Returns the surface_density experienced by the spin hitting the surface represented by a Reflection. The geometry has to be a FixedGeometry.

dwell_time(geometry, reflection)

Returns the dwell_time experienced by the spin hitting the surface represented by a Reflection. The geometry has to be a FixedGeometry.

permeability(geometry, reflection)

Returns the permeability experienced by the spin hitting the surface represented by a Reflection. The geometry has to be a FixedGeometry.

GlobalProperties(; R1=0, R2=0, off_resonance=0)

Stores global MRI and collision properties.

R1(snapshot, geometry, global_properties)
R1(spins, geometry, global_properties)
R1(positions, geometry, global_properties)

Returns the R1 experienced by the Spin objects given the surface and volume properties of the FixedGeometry. Alternatively, the position of the spins can be provided. In that case the spins will be presumed to be free.

R2(snapshot, geometry, global_properties)
R2(spins, geometry, global_properties)
R2(positions, geometry, global_properties)

Returns the R2 experienced by the Spin objects given the surface and volume properties of the FixedGeometry. Alternatively, the position of the spins can be provided. In that case the spins will be presumed to be free.

off_resonance(snapshot, geometry, global_properties)
off_resonance(spin(s), geometry, global_properties)
off_resonance(position(s), geometry, global_properties)

Computes the off-resonance field experienced by each spin in the Snapshot.

A tuple is returned with:

1. the off-resonance field due to susceptibility sources in ppm (e.g., myelin/iron).
2. the off-resonance field due to other sources in kHz (i.e., off_resonance=... set by the user).

run_main(args=ARGS; kwargs...)

Main run command for the command line interface.

Any keyword arguments are passed on to ArgParse.ArgParseSettings.

This package supports the running of MR Monte Carlo simulations.

In these simulations hundreds of thousands or millions of Spin particles randomly diffuse through some tissue microstructure. At each timepoint these spins are represented as a Snapshot object. The spin diffusion is constrained by an ObstructionGroup (represented internally as FixedGeometry) The spins of these particles will be evolved based on the Bloch equations with the field strength and relaxation rates set by the local geometry and the effect of one or more MRIBuilder.Sequence objects. All these variables are combined into a single Simulation object. See Simulation for how to run the simulation.

Plotting support for the sequence and resulting signal is also available based on Makie.jl.

Defines biophysical constants used in simulator.

• gyromagnetic_ratio

Defines methods shared across multiple sub-modules.

• get_time
• norm_angle
• MCMRSimulator.get_rotation
• off_resonance

Types:

• Snapshot
• Spin
• SpinOrientation
• FixedXoshiro

Methods:

• longitudinal
• transverse
• phase
• orientation
• isinside
• stuck
• stuck_to
• get_sequence

Defines the functions that run the actual simulation:

• readout: get total signal or Snapshot at any MRIBuilder.ADC objects in the sequences.
• evolve: Return a single Snapshot with the state of the simulation at a given time. This snapshot can be used as initialisation for further runs.

All of these functions call evolve_to_time under the hood to actually run the simulation.

Implements the Bloch equations to model the spin magnetisation under influence of a sequence during a single timestep.

Basic interface for setting the simulation properties

For the main interface see Simulation.

Types:

• MCMRSimulator.GlobalProperties

Methods:

• R1
• R2
• off_resonance
• stick_probability

This module contains the base plotting functions.

If Makie is not installed, none of the plotting functions will work. If a Makie backend is installed and importend, concrete methods will be added to these functions by the MakieMCMRSimulator extension. This makes Makie an optional dependency of MCMRSimulator, which will only be required if you want to use the MCMRSimulator plotting capabilities.

In addition to these empty plotting functions, this module defines the PlotPlane for any 2D-projections and helper functions to project onto this plane.

Support for selecting a subset of the total signal.

Types:

• Subset

Functions:

• get_subset

Breaks up one or more MRIBuilder.Sequence objects into individual timesteps.

The main interface is by calling iter_parts, which will result into a vector of some subtypes of SequencePart.

To just get the readouts call MCMRSimulator.get_readouts.

Defines the main Simulation object.

Defines the main interface to set the TimeStep of the simulator.

Defines the command line interface to MCMRSimulator.jl (mcmr).

Functions:

• MCMRSimulator.run_main

Defines command line interface for mcmr geometry

Defines the command line interface to mcmr run.

Defines the microstructural environment with which the spins interact.

User interface for defining geometry.

Defines the types and main methods for user-defined geometries.

Defines a flexible Field interface to set the properties of user-defined geometries.

Defines the properties of a certain geometry type.

Types:

• ObstructionType

Functions:

• fields
• key_value_pairs

Defines ObstructionGroup, which is the parent type of all user-defined geometries.

Defines random_positions_radii.

Defines MCMRSimulator.fix.

Computes a default value for the size scale of different geometries used for timestep calculations.

Defines MCMRSimulator.fix_susceptibility.

Defines load_mesh.

Supports I/O to and from json

Functions:

• write_geometry
• read_geometry_json

Methods that extend the Mesh constructor to convert different types into meshes.

Mesh-specific operations to normalise (i.e., ensure normals point outwards) and split a mesh.

Internal representation of the tissue geometry.

All objects are immutable and optimised for both memory usage and computational speed.

Methods to extract properties from the FixedGeometry.

MRI properties (using MRIProperties):

• R1
• R2
• off_resonance

Collision properties:

• permeability
• surface_relaxation
• surface_density
• dwell_time

Other functions:

• max_timestep_sticking

Defines ray_grid_intersections, which computes the intersections between a spin path and some grid.

Defines the bounding box and its interface.

Types:

• BoundingBox: N-dimensional box fully containing an obstruction.

Functions

• MCMRSimulator.isinside: returns whether point/spin is inside a bounding box (or other obstruction)
• could_intersect: returns whether spin trajectory could intersect with obstruction in bounding box

Defines the Intersection between a spin path and the geometry.

Defines all the base obstructions supported by MCMRSimulator.

They are all sub-types of FixedObstruction.

Base types for the API with individually fixed obstructions.

Types:

• FixedObstruction

Methods:

• detect_intersection
• has_inside
• isinside
• radius
• BoundingBox
• obstruction_type
• random_surface_positions

Defines the intersection between a spin path and an obstruction.

Types:

• ObstructionIntersection

Attributes:

• empty_obstruction_intersections

Methods:

• has_intersection

Defines cylinders and spheres and how to intersect with them.

All base obstructions are assumed to be centered on the origin (except for triangles in the mesh).

This module defines the Shift from this origin to the actual position in space.

Triangle within a mesh and how to intersect with them.

Infinite walls and how to intersect with them.

Defines the fixed versions of MCMRSimulator.ObstructionGroup and its methods (including intersection detection).

Types:

• MCMRSimulator.FixedGeometry
• MCMRSimulator.FixedObstructionGroup

Methods:

• has_inside
• isinside
• detect_intersection
• size_scale

Computes the reflection of the spin off some geometry.

This is also used internally to represent a bound spin.

Types:

• Reflection

Methods:

• direction
• has_hit
• previous_hit

Base module for computing the off-resonance field from magnetic susceptibility sources.

Defines parent type and functions for all magnetic susceptibility sources.

Types:

• BaseSusceptibility

Functions:

• single_susceptibility
• single_susceptibility_gradient

Off-resonance field due to magnetic susceptibility of an hollow cylinder, which is approximated as being infinitely thin.

Computes the off-resonance field produced by a triangular mesh.

Off-resonance field due to magnetic susceptibility of an hollow cylinder.

Defines how to compute the off-resonance field for distant magnetic susceptibility sources.

The contribution from very distant sources are pre-computed on a SusceptibilityGrid. Intermediate sources are represented as point magnetic susceptibility elements (SusceptibilityGridElement).

Defines a HitGrid used to precompute the locations of each obstruction within a grid.

This helps limit the number of intersection calculations required for each spin movement.