larmor.simulation#

Simulation of diffusion MRI from a phantom whose ground truth is known.

Use it to make diffusion-weighted images for testing a method, where the true fibre directions and tissue fractions of every voxel are known. A simulation has three parts. A Phantom places tissues on a voxel grid: each Tissue is a sum of water compartments (Compartment), such as sticks for axons, and the tissues fill the voxels through fibre bundles (Bundle, drawn from a BundleGeometry or taken from a tractogram) and tissue maps (TissueMap). A Protocol describes the scan: the gradient table, the pulse timings, the echo and repetition times. simulate then computes the images with the instant engine, which sums the closed-form signal of each compartment in each voxel, and save_simulation writes the images and the ground truth to a folder.

The predefined tissues (standard_model_tissue, sandi_tissue, noddi_tissue, free_water_tissue and csf_tissue) give typical values at 3 T, and example_phantom reads a phantom that ships with Larmor. A phantom and its scan can be saved as a phantom file (write_phantom, read_phantom), which lar simulate runs and the viewer’s Phantom panel edits. signal_breakdown splits the signal of one voxel into its factors, to see where a value comes from.

Two more engines are planned: a k-space engine for the effects of the readout and a Monte Carlo engine for random walkers in a substrate. Their entry points, larmor.simulation.kspace.simulate and larmor.simulation.montecarlo.simulate, raise NotImplementedError until Larmor 1.3 and 1.4.

`EXAMPLE_PHANTOMS`	The names of the bundled phantoms.
`FREE_WATER_DIFFUSIVITY`	Convert a string or number to a floating-point number, if possible.
`Bundle`	A fibre bundle of a phantom, and the tissue that fills it.
`BundleGeometry`	The shape of a synthetic bundle.
`BundleShare`	A bundle's share of the voxel in a `SignalBreakdown`.
`Compartment`	One water pool of a tissue: its diffusion model, fraction and relaxation.
`CompartmentModel`	The diffusion models a compartment can follow.
`CompartmentSignal`	One compartment's share of a voxel's signal, factor by factor.
`Phantom`	An object to simulate: tissues placed on a voxel grid.
`PhantomFile`	What a phantom file holds, as `read_phantom` returns it.
`Protocol`	A diffusion-weighted scan: the gradient table and the pulse timings.
`PulseSequence`	The kind of echo, which decides whether T2 or T2* weights the signal.
`SignalBreakdown`	The noise-free signal of one voxel in one volume, split into its factors.
`SimulatedDWI`	The result of `simulate`: the images, the ground truth and the inputs.
`SimulationFiles`	The files `save_simulation` writes, all inside one folder.
`Tissue`	A named tissue: compartments whose fractions sum to 1.
`TissueMap`	The fraction of one tissue in each voxel, such as a grey matter probability.
`VoxelPartition`	How each voxel of a phantom is split between its tissues and fibres.
`add_noise`	Adds complex Gaussian noise to noise-free signals and takes the magnitude.
`create_bundle`	Creates the streamlines of a synthetic bundle.
`create_gradient_table`	Creates a multi-shell gradient table with evenly spread directions.
`csf_tissue`	Cerebrospinal fluid: `free_water_tissue` with the default name "CSF".
`default_tissues`	The predefined tissues: WM (standard model), GM (SANDI) and CSF (free water).
`example_phantom`	Reads a phantom that ships with Larmor, with its scan.
`free_water_tissue`	Free water: one FREE_WATER compartment, a ball with D = 3e-3 mm^2/s.
`noddi_tissue`	NODDI: Watson-dispersed sticks, a hindered tensor around them, free water.
`partition_voxels`	Splits each voxel of a phantom between its tissues, fibres and background.
`read_phantom`	Reads a phantom file written by `write_phantom`.
`sandi_tissue`	SANDI grey matter: a sphere (soma), a powder-averaged stick and a ball.
`save_simulation`	Writes a simulation and its ground truth into one self-contained folder.
`signal_breakdown`	Splits one voxel's noise-free signal in one volume into its factors.
`simulate`	Simulates the diffusion-weighted images of a phantom, voxel by voxel.
`simulation_files`	The files a simulation of this phantom and protocol is saved as.
`standard_model_tissue`	The standard model of white matter: stick, zeppelin and optional free water.
`write_phantom`	Writes a phantom file (JSON, suggested suffix .phantom.json).

larmor.simulation.EXAMPLE_PHANTOMS = ('crossing',)#: The names of the bundled phantoms.

larmor.simulation.FREE_WATER_DIFFUSIVITY = 0.003#: Convert a string or number to a floating-point number, if possible.

class larmor.simulation.Bundle(source, tissue='WM', path=None, streamline_area=None, densest_fill=0.7)#

A fibre bundle of a phantom, and the tissue that fills it.

A bundle places its tissue in the voxels it passes through, and its streamlines give the fibre directions there. The share of a voxel that a bundle fills is its occupancy, between 0 and 1, and it is found in one of two ways. A bundle drawn from a geometry fills each voxel by the share of the voxel inside its tube (see BundleGeometry), so a straight bundle is full inside and partly full only at its edge. A tractogram, or a geometry with a streamline_area, fills each voxel with the occupancy L a / V, capped at 1. Here L is the length of its streamlines in the voxel (mm), a is the cross-section each streamline stands for (mm^2), and V is the voxel volume (mm^3).

Parameters:

source (BundleGeometry | Tractogram) – A geometry to draw the streamlines from, or a tractogram (for example from load_tractogram), in world coordinates (mm).
tissue (str) – The name of the tissue that fills the bundle. Defaults to “WM”.
path (Path | None) – The file a tractogram source was read from; write_phantom stores it. Defaults to None.
streamline_area (float | None) – The cross-section a that each streamline stands for, in mm^2. Defaults to None: a geometry then fills its tube, and a tractogram gets the value that gives the bundle’s dense voxels an occupancy of densest_fill. Counting streamline length gives random dips where fewer streamlines happen to pass, so a geometry is best left with None.
densest_fill (float) – The occupancy of a tractogram bundle’s dense voxels when streamline_area is None, in (0, 1]. A voxel is “dense” when the bundle’s streamline length in it reaches the 95th percentile of that length over the voxels the bundle reaches. The percentile is used instead of the maximum so that a few outlier voxels (at seeds, endpoints or bottlenecks) do not set the scale; voxels above it are capped at 1. A tractogram carries no physical density, so this value sets it. Defaults to 0.7, a choice that stands for densely packed white matter with room for other tissue. Not used for a geometry.

Raises:

ValueError – If streamline_area is not positive or densest_fill is not in (0, 1].

streamlines(index=0)#

The bundle’s streamlines: drawn from the geometry, or the tractogram.

Parameters:: index (int) – The bundle’s index, stored in streamline_data[“bundle”] of the streamlines drawn from a geometry. Defaults to 0.
Returns:: The streamlines in world coordinates (mm).
Return type:: Tractogram

class larmor.simulation.BundleGeometry(control_points, num_streamlines=1000, radius=2.0, num_points=100, fan_start=0.0, fan_end=0.0, fan_sharpness=2.0, dispersion=0.0, seed=None)#

The shape of a synthetic bundle.

A geometry describes a bundle by a few control points along its centreline, a radius, and how much the bundle fans out at its ends. create_bundle draws its streamlines, and a Bundle made from it places its tissue in a phantom.

The centreline is an interpolating spline (cubic, or of lower degree for fewer than four points) through the control points. Each streamline is the same kind of spline through the control points, shifted by an offset across the centreline. The offset is drawn uniformly in a disc of radius and scaled by the fan factor 1 + fan_start (1 - u)^p + fan_end u^p, where u in [0, 1] is the position along the centreline and p is fan_sharpness. A random wander of standard deviation dispersion is added at the inner control points.

Parameters:

control_points (ndarray) – (K, 3) points the centreline passes through, K >= 2, in world coordinates (RAS+) in mm.
num_streamlines (int) – The number of streamlines. Defaults to 1000.
radius (float) – The bundle’s radius in mm. Defaults to 2.0.
num_points (int) – Points per streamline. Defaults to 100.
fan_start (float) – Extra spread at the first control point, where the disc has the radius radius (1 + fan_start). Defaults to 0.0.
fan_end (float) – Extra spread at the last control point, in the same way. Defaults to 0.0.
fan_sharpness (float) – p, how close to the ends the fanning stays; a larger p keeps it nearer the ends. Defaults to 2.0.
dispersion (float) – Standard deviation of the wander, mm. Defaults to 0.0.
seed (int | None) – The seed of the random offsets and wander, so that the same seed gives the same streamlines; None draws new ones each time. Defaults to None.

Raises:

ValueError – If the control points are not (K, 3) finite values with K >= 2, or a count or size is out of range.

to_dict()#

The geometry as plain values that can be saved as JSON.

Returns:: A dict of the fields, with the control points as lists.
Return type:: dict[str, object]

classmethod from_dict(data)#

Builds a geometry from the values to_dict gives.

Parameters:: data (Mapping[str, object]) – The dict.
Returns:: The geometry.
Raises:: ValueError – If a value is out of range.
Return type:: BundleGeometry

class larmor.simulation.BundleShare(index, name, tissue, share, angle_deg)#

A bundle’s share of the voxel in a SignalBreakdown.

index#

The bundle’s index in the phantom.

Type:: int

name#

The bundle’s name.

Type:: str

tissue#

The tissue that fills it.

Type:: str

share#

The share of the voxel it fills.

Type:: float

angle_deg#

The angle between its main fibre direction in the voxel and the gradient, 0 to 90 degrees; None at b = 0.

Type:: float | None

class larmor.simulation.Compartment(model, fraction, diffusivity, *, t1_ms, t2_ms, t2star_ms, proton_density=1.0, perpendicular_diffusivity=0.0, radius_um=0.0, kappa=0.0)#

One water pool of a tissue: its diffusion model, fraction and relaxation.

A tissue is a sum of compartments, such as the water inside axons and the water around them. Each compartment follows one diffusion model, which gives its signal attenuation. The models and their attenuation for b-value b (s/mm^2), unit gradient direction g and unit fibre direction n are:

STICK: exp(-b D (g.n)^2), diffusion along the fibre only (Behrens et al. 2003).
ZEPPELIN: exp(-b (D_perp + (D - D_perp) (g.n)^2)), an axially symmetric tensor (Kaden et al. 2016, Eq. 5; Novikov et al. 2019, Eq. 3.5). D_perp > D is allowed (an oblate tensor), although white matter has D_perp < D.
BALL: exp(-b D), free isotropic diffusion (Novikov et al. 2019, Eq. 3.5).
POWDER_STICK: sqrt(pi) / (2 sqrt(b D)) erf(sqrt(b D)), the stick averaged over all fibre directions (Kroenke et al. 2004; Kaden et al. 2016).
SPHERE: diffusion restricted to an impermeable sphere of radius R, in the Gaussian phase approximation for rectangular pulses, as SANDI uses for the soma (Palombo et al. 2020, Eq. 8). It uses the pulse timings (see Protocol for the waveform a gradient echo is assumed to have).
WATSON_STICK: sticks of diffusivity D whose directions spread about the fibre direction with the Watson density f(n) = exp(kappa (mu.n)^2) / M(1/2, 3/2, kappa), so that the signal is the stick’s averaged over f (Zhang et al. 2012). kappa = 0 spreads them evenly, which gives the powder-averaged stick, and a large kappa gives the stick. This is the neurite compartment of NODDI.
CYLINDER: diffusion restricted to an impermeable cylinder of radius R along the fibre, in the Gaussian phase approximation for rectangular pulses across it (Van Gelderen et al. 1994; De Santis et al. 2016, Eq. 3) and free with D along it. It uses the pulse timings, and a very thin cylinder is the stick.
DOT: a signal of 1 at every b-value, for stationary water trapped in glial cells and other small structures (Panagiotaki et al. 2012). D is unused.
FREE_WATER: exp(-b D), a ball that the ground truth counts as free water, as in free-water imaging (Pasternak et al. 2009). The presets give it D = 3e-3 mm^2/s.

Each compartment also has its own relaxation times and proton density PD. Its signal at b = 0 is PD (1 - exp(-TR / T1)) exp(-TE / T2) for a spin echo, with T2* in place of T2 for a gradient echo.

Parameters:

model (CompartmentModel) – The diffusion model.
fraction (float) – The compartment’s share of its tissue’s volume, in [0, 1].
diffusivity (float) – D in mm^2/s; for the stick and the zeppelin, the diffusivity along the fibre.
t1_ms (float) – T1 in ms.
t2_ms (float) – T2 in ms.
t2star_ms (float) – T2* in ms, at most T2.
proton_density (float) – The water content relative to pure water, in [0, 1]. Defaults to 1.0.
perpendicular_diffusivity (float) – The zeppelin’s D_perp across the fibre, in mm^2/s. Defaults to 0.0; only the zeppelin uses it.
radius_um (float) – The radius of the sphere or the cylinder in um. Defaults to 0.0; only the sphere and the cylinder use it.
kappa (float) – The Watson concentration of WATSON_STICK, >= 0. Defaults to 0.0; only WATSON_STICK uses it.

Raises:

ValueError – If the model is unknown or a value is out of range.

References

Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61 (2012) 1000-1016. Van Gelderen P, DesPres D, van Zijl PCM, Moonen CTW. Evaluation of restricted diffusion in cylinders. Phosphocreatine in rabbit leg muscle. J Magn Reson B 103 (1994) 255-260. De Santis S, Jones DK, Roebroeck A. Including diffusion time dependence in the extra-axonal space improves in vivo estimates of axonal diameter and density in human white matter. NeuroImage 130 (2016) 91-103. Novikov DS, Fieremans E, Jespersen SN, Kiselev VG. Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation. NMR Biomed 32 (2019) e3998. Panagiotaki E, Schneider T, Siow B, Hall MG, Lythgoe MF, Alexander DC. Compartment models of the diffusion MR signal in brain white matter: A taxonomy and comparison. NeuroImage 59 (2012) 2241-2254. Pasternak O, Sochen N, Gur Y, Intrator N, Assaf Y. Free water elimination and mapping from diffusion MRI. Magn Reson Med 62 (2009) 717-730. Behrens TEJ, Woolrich MW, Jenkinson M, et al. Characterization and propagation of uncertainty in diffusion-weighted MR imaging. Magn Reson Med 50 (2003) 1077-1088. Kroenke CD, Ackerman JJH, Yablonskiy DA. On the nature of the NAA diffusion attenuated MR signal in the central nervous system. Magn Reson Med 52 (2004) 1052-1059. Kaden E, Kelm ND, Carson RP, Does MD, Alexander DC. Multi-compartment microscopic diffusion imaging. NeuroImage 139 (2016) 346-359. Palombo M, Ianus A, Guerreri M, et al. SANDI: a compartment-based model for non-invasive apparent soma and neurite imaging by diffusion MRI. NeuroImage 215 (2020) 116835.

property is_oriented: bool#: Whether the signal depends on the fibre direction.

property engine_parameters: tuple[float, float, float, float]#: D, D_perp, radius_um and kappa, in the order the compiled engine takes.

class larmor.simulation.CompartmentModel#

The diffusion models a compartment can follow.

Each model gives the signal attenuation of one water pool; Compartment lists their formulas.

STICK = 'stick'#

ZEPPELIN = 'zeppelin'#

BALL = 'ball'#

POWDER_STICK = 'powder_stick'#

SPHERE = 'sphere'#

WATSON_STICK = 'watson_stick'#

CYLINDER = 'cylinder'#

DOT = 'dot'#

FREE_WATER = 'free_water'#

property code: int#: The integer that stands for the model in Larmor’s compiled engine.

property is_oriented: bool#: Whether the signal depends on a fibre direction (only bundles have one).

class larmor.simulation.CompartmentSignal(label, tissue, model, bundle, partial_volume, fraction, proton_density, t1_factor, t2_factor, attenuation)#

One compartment’s share of a voxel’s signal, factor by factor.

SignalBreakdown.compartments holds one per compartment in the voxel. The signal of the compartment is the product of its factors, applied in this order. The partial volume, the fraction and the proton density together give the water before any decay. The T1 and T2 factors give the relaxation, and the attenuation gives the diffusion weighting.

label#

The compartment in plain words, such as “WM axons”.

Type:: str

tissue#

The name of its tissue.

Type:: str

model#

Its diffusion model.

Type:: larmor.simulation.tissues.CompartmentModel

bundle#

The index of the bundle whose share of the voxel this is, or None for the share the tissue maps and the background give.

Type:: int | None

partial_volume#

The share of the voxel its tissue fills here (the bundle’s share, for a bundle).

Type:: float

fraction#

The compartment’s fraction of its tissue.

Type:: float

proton_density#

Its proton density.

Type:: float

t1_factor#

1 - exp(-TR / T1).

Type:: float

t2_factor#

exp(-TE / T2), with T2* for a gradient echo.

Type:: float

attenuation#

The diffusion attenuation at the volume’s b-value and direction, averaged over the bundle’s fibre directions in the voxel, each weighted by its share. Oriented compartments (stick, zeppelin, Watson-dispersed sticks, cylinder) occur only in bundles: tissue maps and the background cannot hold an oriented tissue.

Type:: float

property no_decay: float#

partial volume x fraction x proton density.

Type:: The water before decay

property after_decay: float#: The water after the T1 and T2 (or T2*) factors.

property signal: float#: The signal after the diffusion weighting too.

class larmor.simulation.Phantom(grid, tissues=<factory>, bundles=(), tissue_maps=(), background=None)#

An object to simulate: tissues placed on a voxel grid.

A phantom is the input of simulate. Its tissues are placed in the voxels by fibre bundles, by tissue maps, and by a background tissue that fills the rest. Each voxel is split between the tissues in three steps, which never scale a fraction up:

Each bundle takes its occupancy (see Bundle). Where the bundles add up to more than 1 (crossings), they are scaled down to sum to 1.
The tissue maps share what the bundles leave, r = 1 - sum of occupancies. Where the map values add up to more than r, they are scaled down to sum to r.
The rest goes to the background: the tissue named by background, or no signal when it is None.

A voxel filled by one tissue therefore has exactly that tissue’s signal, and a voxel with a grey matter map value of 0.05 and nothing else is 5% grey matter and 95% background.

Parameters:

grid (Grid) – The voxel grid.
tissues (tuple[Tissue, ...]) – The named tissues. Defaults to WM, GM and CSF (default_tissues).
bundles (tuple[Bundle, ...]) – The fibre bundles. Defaults to none.
tissue_maps (tuple[TissueMap, ...]) – At most one map per tissue. Defaults to none.
background (str | None) – The name of the tissue that fills what remains, or None for no signal. Its compartments must not need a fibre direction. Defaults to None.

Raises:

ValueError – If tissue names repeat, a bundle, map or the background names an unknown tissue, a map or the background uses a tissue that needs a fibre direction, two maps share a tissue, or a map is not a 3D volume on the grid with values in [0, 1].

tissue(name)#

The tissue with this name.

Parameters:: name (str) – The tissue’s name.
Returns:: The tissue.
Raises:: ValueError – If the phantom has no tissue of that name.
Return type:: Tissue

class larmor.simulation.PhantomFile(phantom, protocol=None, viewer=<factory>)#

What a phantom file holds, as read_phantom returns it.

phantom#

The phantom.

Type:: larmor.simulation.phantom.Phantom

protocol#

The scan, if the file has one.

Type:: larmor.simulation.protocol.Protocol | None

viewer#

The settings of the viewer’s Phantom panel, such as the bundles’ names, kept as they are; empty if none.

Type:: collections.abc.Mapping[str, Any]

property bundle_names: list[str] | None#

The bundles’ names from the viewer settings, in order.

Returns:: The names, or None unless every bundle has one.

class larmor.simulation.Protocol(gradients, small_delta_ms, big_delta_ms, te_ms, tr_ms, sequence=PulseSequence.SPIN_ECHO)#

A diffusion-weighted scan: the gradient table and the pulse timings.

Every simulation engine takes a protocol with the phantom. The diffusion weighting is a pair of rectangular gradient pulses (pulsed gradient spin echo), so b = gamma^2 G^2 delta^2 (Delta - delta / 3). A gradient echo has no refocusing pulse, so it is assumed to use a bipolar pair instead: a lobe of +G, then one of -G, each delta long and Delta apart. Its effective waveform is the same, and so are its b-value and its restricted signal.

The protocol keeps its own read-only copy of the gradient table, so change a table before building a protocol from it, not after.

Parameters:

gradients (GradientTable) – The b-value (s/mm^2) and direction of each volume.
small_delta_ms (float) – delta, the duration of each gradient pulse, in ms.
big_delta_ms (float) – Delta, the time between the starts of the two pulses, in ms.
te_ms (tuple[float, ...]) – The echo time, or a sequence of echo times for a multi-echo acquisition, in ms. Every echo shares the one diffusion weighting, which comes before the first echo.
tr_ms (float) – The repetition time, in ms.
sequence (PulseSequence) – SPIN_ECHO (relaxation with T2) or GRADIENT_ECHO (with T2*). Defaults to SPIN_ECHO.

Raises:

ValueError – If a timing is not positive and finite, delta > Delta, Delta + delta >= the first TE, or TR <= the last TE.

property num_echoes: int#: The number of echo times.

to_dict()#

The protocol as plain values that can be saved as JSON.

Returns:: A dict with the b-values, the directions and the timings.
Return type:: dict[str, object]

classmethod from_dict(data)#

Builds a protocol from the values to_dict gives.

Parameters:

data (Mapping[str, object]) – The dict.

Returns:

The protocol.

Raises:

KeyError – If a value is missing.
ValueError – If the sequence is unknown or the values do not make a valid protocol.

Return type:

Protocol

class larmor.simulation.PulseSequence#

The kind of echo, which decides whether T2 or T2* weights the signal.

A spin echo refocuses the dephasing from field differences, so its signal decays with T2; a gradient echo does not, so its signal decays with T2*.

SPIN_ECHO = 'spin_echo'#

GRADIENT_ECHO = 'gradient_echo'#

class larmor.simulation.SignalBreakdown(voxel, volume, echo, b_value, direction, te_ms, tr_ms, sequence, tissue_fractions, bundles, compartments, sigma)#

The noise-free signal of one voxel in one volume, split into its factors.

signal_breakdown computes it, and larmor.plot.plot_signal_breakdown draws it. Its properties give the signal at three steps: the water in the voxel before decay, the signal after relaxation, and the signal after the diffusion weighting.

voxel#

The voxel (i, j, k).

Type:: tuple[int, int, int]

volume#

The volume’s index in the protocol’s gradient table.

Type:: int

echo#

The echo’s index in the protocol’s echo times.

Type:: int

b_value#

The volume’s b-value, s/mm^2.

Type:: float

direction#

The volume’s unit gradient direction (zero at b = 0).

Type:: tuple[float, float, float]

te_ms#

The echo time, ms.

Type:: float

tr_ms#

The repetition time, ms.

Type:: float

sequence#

The pulse sequence.

Type:: larmor.simulation.protocol.PulseSequence

tissue_fractions#

A (name, share) pair for each tissue of the phantom in the voxel, the background tissue included.

Type:: tuple[tuple[str, float], …]

bundles#

The bundles in the voxel.

Type:: tuple[larmor.simulation.breakdown.BundleShare, …]

compartments#

The compartments with a share of the voxel, tissue by tissue; a tissue that fills bundles is listed once per bundle.

Type:: tuple[larmor.simulation.breakdown.CompartmentSignal, …]

sigma#

The noise standard deviation; None without noise.

Type:: float | None

property no_decay: float#

the sum of the compartments’ water.

Type:: The water in the voxel

property after_decay: float#

the signal at b = 0.

Type:: The signal after relaxation

property signal: float#: The noise-free signal after the diffusion weighting.

property noise_floor: float | None#: The mean magnitude of pure noise, sigma sqrt(pi / 2); None without noise.

property expected_magnitude: float | None#: The mean of the noisy (Rician) magnitude; None without noise.

class larmor.simulation.SimulatedDWI(echoes, noise_free, gradients, echo_times_ms, partition, sigma, phantom=None, protocol=None, snr=None, noise_reference='WM', seed=None)#

The result of simulate: the images, the ground truth and the inputs.

save_simulation writes it to a folder.

echoes#

One 4D float32 volume (X, Y, Z, N) per echo time, with noise if an SNR was given.

Type:: tuple[larmor.types.volume.Volume, …]

noise_free#

The same without noise.

Type:: tuple[larmor.types.volume.Volume, …]

gradients#

The gradient table of the N volumes.

Type:: larmor.types.gradient_table.GradientTable

echo_times_ms#

The echo time of each volume in echoes, ms.

Type:: tuple[float, …]

partition#

The ground truth: each voxel’s tissue fractions and fibres, and the streamlines drawn for each bundle (see VoxelPartition).

Type:: larmor.simulation.phantom.VoxelPartition

sigma#

The noise standard deviation, or None without noise.

Type:: float | None

phantom#

The phantom that was simulated, with a seed on every bundle geometry, so that it gives the same streamlines again.

Type:: larmor.simulation.phantom.Phantom | None

protocol#

The scan that was simulated.

Type:: larmor.simulation.protocol.Protocol | None

snr#

The signal-to-noise ratio asked for, or None without noise.

Type:: float | None

noise_reference#

The name of the tissue that defines the noise level.

Type:: str

seed#

The seed of the noise and the bundles: the one given, or the one drawn when none was given.

Type:: int | None

class larmor.simulation.SimulationFiles(folder, dwi, bvals, bvecs, noise_free, tissue_fractions, background_fraction, free_water_fraction, bundle_fractions, compartment_fractions, signal_fractions, free_water_signal_fraction, peaks, masks, bundle_mask, bundles, tissue_maps, phantom, info)#

The files save_simulation writes, all inside one folder.

simulation_files gives these paths without writing anything.

folder#

The output folder.

Type:: pathlib.Path

dwi#

The DWI, one 4D NIfTI per echo time.

Type:: tuple[pathlib.Path, …]

bvals#

The b-values, shared by the echoes (FSL format).

Type:: pathlib.Path

bvecs#

The b-vectors in the FSL frame, shared by the echoes.

Type:: pathlib.Path

noise_free#

The noise-free DWI per echo; empty if the simulation had no noise.

Type:: tuple[pathlib.Path, …]

tissue_fractions#

4D, the fraction of each tissue in each voxel, in the order of the phantom’s tissues.

Type:: pathlib.Path

background_fraction#

3D, the fraction of the voxel that no tissue fills.

Type:: pathlib.Path

free_water_fraction#

3D, the volume fraction of each voxel that free water (the FREE_WATER compartments) fills, the background included.

Type:: pathlib.Path

bundle_fractions#

4D, the fraction of each voxel that each bundle fills; None if the phantom has no bundle.

Type:: pathlib.Path | None

peaks#

The ground-truth fibre directions, in the layout lar peaks writes. Each peak takes 3 volumes, the x, y and z of its direction scaled by its amplitude, which is the share of the voxel that the fibres of the peak fill. The peaks are sorted largest first, with NaN where a voxel has fewer peaks.

Type:: pathlib.Path

compartment_fractions#

4D, the volume fraction of each voxel that each compartment fills: its tissue’s fraction times its share of the tissue, where the background tissue’s fraction includes the background. The volumes are in the order of the tissues and, within a tissue, of its compartments, as listed in simulation.json.

Type:: pathlib.Path

signal_fractions#

One 4D file per echo time: each compartment’s share of the voxel’s signal at b = 0, f_c w_c / sum over compartments of f w, with f_c its volume fraction and w_c = PD (1 - exp(-TR / T1)) exp(-TE / T2) (T2* for a gradient echo) its relaxation weight; 0 where the voxel has no signal. The volumes are in the order of compartment_fractions. A model fitted to the images (free-water imaging, NODDI, the standard model) recovers these fractions, not the volume fractions, so compare a fit with them.

Type:: tuple[pathlib.Path, …]

free_water_signal_fraction#

One 3D file per echo time: the share of the voxel’s signal at b = 0 that comes from free water, the sum of the FREE_WATER volumes of signal_fractions.

Type:: tuple[pathlib.Path, …]

masks#

One 3D mask per tissue, by tissue name.

Type:: dict[str, pathlib.Path]

bundle_mask#

One 3D mask of the voxels the tissues of the bundles fill together; None if the phantom has no bundle.

Type:: pathlib.Path | None

bundles#

One .tck per bundle, in the order of the phantom’s bundles.

Type:: tuple[pathlib.Path, …]

tissue_maps#

A copy of each tissue map of the phantom, by tissue name.

Type:: dict[str, pathlib.Path]

phantom#

The phantom file, with the scan, which lar simulate can run again.

Type:: pathlib.Path

info#

simulation.json, which records the protocol, the noise, the seed and what each file holds.

Type:: pathlib.Path

all_paths()#

Every file of the folder, in the order they are written.

class larmor.simulation.Tissue(name, compartments)#

A named tissue: compartments whose fractions sum to 1.

A phantom places tissues in its voxels, and the bundles, tissue maps and background of a phantom refer to a tissue by its name. Functions such as standard_model_tissue and sandi_tissue make tissues with typical values; build a tissue from Compartment objects for any other mix.

Parameters:

name (str) – The tissue’s name, such as “WM”; bundles and tissue maps refer to it.
compartments (tuple[Compartment, ...]) – Its compartments.

Raises:

ValueError – If the name is empty, there are no compartments, or the fractions do not sum to 1 (within 1e-6).

property is_oriented: bool#: Whether a compartment needs a fibre direction, so only bundles can use it.

to_dict()#

The tissue as plain values that can be saved as JSON.

Returns:: A dict with the name and a list of compartments.
Return type:: dict[str, object]

classmethod from_dict(data)#

Builds a tissue from the values to_dict gives.

Parameters:

data (Mapping[str, object]) – The dict.

Returns:

The tissue.

Raises:

TypeError – If the compartments are not a list of objects.
ValueError – If a model is unknown or a value is out of range.

Return type:

Tissue

class larmor.simulation.TissueMap(tissue, volume, path=None)#

The fraction of one tissue in each voxel, such as a grey matter probability.

A tissue map places a tissue without fibres, such as grey matter or CSF, in a phantom; the bundles place the tissues that need a fibre direction.

Parameters:

tissue (str) – The tissue’s name. Its compartments must not need a fibre direction, as a map gives none.
volume (Volume) – 3D values in [0, 1] on the phantom’s grid.
path (Path | None) – The file the map was read from; write_phantom stores it. Defaults to None.

class larmor.simulation.VoxelPartition(grid, tissue_names, tissue_fractions, background, bundle_fractions, offsets, fibre_tissue, fibre_bundle, fibre_weight, fibre_direction, fibre_eigenvalues, streamlines)#

How each voxel of a phantom is split between its tissues and fibres.

This is the ground truth of a simulation: partition_voxels computes it, and simulate returns it in SimulatedDWI.partition. A fibre is one direction of one bundle in one voxel; a bundle that bends or fans out in a voxel has several fibres there. The fibres of all voxels are stored in flat arrays of length F, sorted by voxel, and offsets says where each voxel’s fibres are.

grid#

The phantom’s grid, with V voxels.

Type:: larmor.types.grid.Grid

tissue_names#

The phantom’s tissue names, in the order of the columns.

Type:: tuple[str, …]

tissue_fractions#

(X, Y, Z, T) the fraction of each tissue in each voxel, from bundles and maps; the background is not included.

Type:: numpy.ndarray

background#

(X, Y, Z) the rest, 1 - the sum of tissue_fractions.

Type:: numpy.ndarray

bundle_fractions#

(X, Y, Z, B) each bundle’s occupancy.

Type:: numpy.ndarray

offsets#

(V + 1,) int64. The fibres of voxel v, with voxels numbered in C order, are those from offsets[v] to offsets[v + 1] - 1.

Type:: numpy.ndarray

fibre_tissue#

(F,) int64 the tissue index of each fibre.

Type:: numpy.ndarray

fibre_bundle#

(F,) int64 the index of each fibre’s bundle.

Type:: numpy.ndarray

fibre_weight#

(F,) the fraction of its voxel each fibre fills.

Type:: numpy.ndarray

fibre_direction#

(F, 3) each fibre’s unit direction in world coordinates.

Type:: numpy.ndarray

fibre_eigenvalues#

(F, 3) the eigenvalues of the structure tensor of each fibre’s bundle in its voxel, descending (see partition_voxels).

Type:: numpy.ndarray

streamlines#

Each bundle’s streamlines, as drawn or given.

Type:: tuple[larmor.types.tractogram.Tractogram, …]

larmor.simulation.add_noise(volumes, sigma, *, seed=None)#

Adds complex Gaussian noise to noise-free signals and takes the magnitude.

simulate calls it when given an SNR; call it directly to add noise of another level to the noise-free images. Each value S becomes sqrt((S + n_1)^2 + n_2^2), where n_1 and n_2 are independent and normal with standard deviation sigma, so the magnitude is Rician, as for a single coil. Like thermal noise, the one sigma holds for every voxel, volume and echo. A zero signal therefore gives a mean of sigma sqrt(pi / 2), and the SNR falls with the signal. The noise is drawn in one sequence, volume after volume, so the result depends only on the seed.

Parameters:

volumes (Sequence[Volume]) – The noise-free signal, for example one volume per echo.
sigma (float) – The noise standard deviation, in the signal’s units, >= 0.
seed (int | None) – Seed of the noise. Defaults to None.

Returns:

The noisy magnitudes, float32, one per input volume, on its grid.

Raises:

ValueError – If sigma is negative or not finite.

Return type:

tuple[Volume, …]

References

Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson Med 34 (1995) 910-914.

larmor.simulation.create_bundle(geometry, *, index=0)#

Creates the streamlines of a synthetic bundle.

Parameters:

geometry (BundleGeometry) – The bundle’s shape and seed.
index (int) – The bundle’s index, stored for every streamline in streamline_data[“bundle”]. Defaults to 0.

Returns:

A tractogram of num_streamlines streamlines of num_points points each, in world coordinates (mm).

Return type:

Tractogram

larmor.simulation.create_gradient_table(shells, *, num_b0=1, seed=None)#

Creates a multi-shell gradient table with evenly spread directions.

Use it to build the gradient table of a Protocol. Each shell’s directions start at random and are spread by electrostatic repulsion on the half sphere: every direction and its opposite carry a charge, and projected gradient descent lowers the energy, the sum over pairs of 1 / norm(x_i - x_j) + 1 / norm(x_i + x_j). The directions are then turned to z >= 0. The b = 0 volumes come first, then each shell in the given order.

Parameters:

shells (Sequence[tuple[float, int]]) – One pair (b-value in s/mm^2, number of directions) per shell.
num_b0 (int) – The number of b = 0 volumes. Defaults to 1.
seed (int | None) – Seed of the starting directions; the same seed gives the same table. Defaults to None.

Returns:

The gradient table, directions in world coordinates.

Raises:

ValueError – If a shell has a b-value that is not positive or fewer than one direction, or num_b0 is negative.

Return type:

GradientTable

References

Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn Reson Med 42 (1999) 515-525.

larmor.simulation.csf_tissue(*, name='CSF', **values)#

Cerebrospinal fluid: free_water_tissue with the default name “CSF”.

Parameters:

name (str) – The tissue’s name. Defaults to “CSF”.
**values (float) – The other arguments of free_water_tissue.

Returns:

The tissue.

Return type:

Tissue

larmor.simulation.default_tissues()#

The predefined tissues: WM (standard model), GM (SANDI) and CSF (free water).

A new Phantom starts with these tissues.

Returns:: The three tissues with their default parameters, for 3 T.
Return type:: tuple[Tissue, …]

larmor.simulation.example_phantom(name='crossing')#

Reads a phantom that ships with Larmor, with its scan.

Use it to try a simulation without building a phantom first. The names are in EXAMPLE_PHANTOMS. “crossing” is two straight bundles of 34 mm crossing at 90 degrees, one along x and one along y, with different standard model parameters and T2, on a 40 x 40 x 40 grid of 1 mm voxels. It is scanned with b = 1000 and 2000 s/mm^2 (32 directions each) and two b = 0 volumes.

Parameters:: name (str) – The phantom’s name. Defaults to “crossing”.
Returns:: The phantom and its protocol.
Raises:: ValueError – If there is no bundled phantom of that name.
Return type:: PhantomFile

larmor.simulation.free_water_tissue(*, name='FW', diffusivity=0.003, t1_ms=4000.0, t2_ms=2000.0, t2star_ms=400.0, proton_density=1.0)#

Free water: one FREE_WATER compartment, a ball with D = 3e-3 mm^2/s.

Its signal is exp(-b D). The ground truth of a simulation counts it in two maps. The free-water fraction is the share of the voxel’s volume it fills. The free-water signal fraction, one map per echo time, is its share of the voxel’s signal at b = 0, which weights each compartment’s volume by its proton density and relaxation; this is the quantity that free-water imaging estimates (Pasternak et al. 2009), and with the defaults below it is larger than the volume fraction, because free water has a higher proton density and a longer T2 than tissue. The relaxation defaults assume 3 T. The default diffusivity, 3.0e-3 mm^2/s, is the value NODDI fixes for its CSF compartment (Zhang et al. 2012). The proton density of 1 (pure water) and T1 4000 ms, T2 2000 ms and T2* 400 ms are choices: no verified source for CSF at 3 T is used here.

Parameters:

name (str) – The tissue’s name. Defaults to “FW”.
diffusivity (float) – The diffusivity, mm^2/s. Defaults to 3.0e-3.
t1_ms (float) – T1, ms. Defaults to 4000.
t2_ms (float) – T2, ms. Defaults to 2000.
t2star_ms (float) – T2*, ms. Defaults to 400.
proton_density (float) – The proton density. Defaults to 1.0.

Returns:

The tissue.

Return type:

Tissue

References

Pasternak O, Sochen N, Gur Y, Intrator N, Assaf Y. Free water elimination and mapping from diffusion MRI. Magn Reson Med 62 (2009) 717-730. Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61 (2012) 1000-1016.

larmor.simulation.noddi_tissue(*, name='WM', intra_fraction=0.6, free_water_fraction=0.05, kappa=8.0, diffusivity=0.0017, free_water_diffusivity=0.003, t1_ms=832.0, t2_ms=69.0, t2star_ms=46.55, proton_density=0.7)#

NODDI: Watson-dispersed sticks, a hindered tensor around them, free water.

Use it for white matter whose fibres spread about a main direction. The model’s signal is (1 - v_iso) (v_ic A_ic + (1 - v_ic) A_ec) + v_iso A_iso (Zhang et al. 2012, Eq. 1). Here v_ic is intra_fraction and v_iso free_water_fraction, which are volume fractions, and the simulation weights each compartment’s signal by its proton density and relaxation as well. Eq. 1 then holds with v_ic and v_iso replaced by signal fractions, which are what a NODDI fit recovers and what the simulation writes as the signal-fraction maps; they equal the volume fractions only when every compartment has the same relaxation weight, which the defaults do not give (free water has a higher proton density and a longer T2). A_ic are sticks of diffusivity d spread by a Watson density of concentration kappa (Eqs. 2 and 3; see Compartment). A_ec is the zeppelin of the extra-neurite space: its diffusivity across a neurite is d (1 - v_ic), a tortuosity model, and averaged over the Watson density it is a zeppelin along the fibre with d_par’ = d - d v_ic (1 - tau1) and d_perp’ = d - d v_ic (1 + tau1) / 2 (Eqs. 4 to 6), with tau1 from Eq. 7. A_iso is free water with the relaxation of free_water_tissue, apart from T2*. The orientation dispersion index of the paper is OD = (2 / pi) arctan(1 / kappa) (Eq. 9). The sticks and the zeppelin share T1, T2, T2* and proton density (NODDI itself has no relaxation); the relaxation defaults assume 3 T. Reversible dephasing (T2’) is a property of the field in the voxel, so the free water shares the tissue’s T2’ = 1 / (1/T2* - 1/T2) and gets T2*_fw = 1 / (1/T2_fw + 1/T2’).

Sources of the defaults:

diffusivity 1.7e-3 and free_water_diffusivity 3.0e-3 mm^2/s: the values NODDI fixes (Zhang et al. 2012).
intra_fraction 0.6, free_water_fraction 0.05 and kappa 8 (OD about 0.08): choices for coherent white matter.
T1 832 ms and T2* 46.6 ms: white matter at 3 T in vivo (Wansapura et al. 1999); T2 69 ms: frontal white matter at 3 T in vivo (Lu et al. 2005, Table 1).
proton_density 0.7: an approximate choice for white matter water content (source to be added).

Parameters:

name (str) – The tissue’s name. Defaults to “WM”.
intra_fraction (float) – v_ic, the neurites’ share of the tissue outside free water. Defaults to 0.6.
free_water_fraction (float) – v_iso, the free-water share. Defaults to 0.05.
kappa (float) – The Watson concentration, >= 0. Defaults to 8.
diffusivity (float) – d, the diffusivity along the neurites, mm^2/s. Defaults to 1.7e-3.
free_water_diffusivity (float) – The free water’s D, mm^2/s. Defaults to 3.0e-3.
t1_ms (float) – T1 of the sticks and the zeppelin, ms. Defaults to 832.
t2_ms (float) – T2 of the sticks and the zeppelin, ms. Defaults to 69.
t2star_ms (float) – T2* of the sticks and the zeppelin, ms. Defaults to 46.6.
proton_density (float) – The proton density of the sticks and the zeppelin. Defaults to 0.7.

Returns:

The tissue – the Watson sticks, the zeppelin, and free water if free_water_fraction > 0.

Raises:

ValueError – If a fraction is outside [0, 1] or a value is out of range.

Return type:

Tissue

References

Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61 (2012) 1000-1016. Wansapura JP, Holland SK, Dunn RS, Ball WS. NMR relaxation times in the human brain at 3.0 tesla. J Magn Reson Imaging 9 (1999) 531-538. Lu H, Nagae-Poetscher LM, Golay X, Lin D, Pomper M, van Zijl PCM. Routine clinical brain MRI sequences for use at 3.0 Tesla. J Magn Reson Imaging 22 (2005) 13-22. doi:10.1002/jmri.20356

larmor.simulation.partition_voxels(phantom, *, seed=None, n_jobs=-1, progress=None)#

Splits each voxel of a phantom between its tissues, fibres and background.

simulate calls it; call it directly to look at the ground truth of a phantom without simulating the images. The split follows the rules of Phantom. A bundle’s fibres in a voxel come from its streamline directions there. The tangents of the streamlines are sampled every tenth of a voxel and pooled in cells of about 7 degrees, and each cell gives one fibre, the length-weighted mean direction of its tangents. Each fibre’s weight is its share of the bundle’s length in the voxel times the bundle’s occupancy. Each fibre also keeps the eigenvalues of the bundle’s length-weighted structure tensor in the voxel, T = sum_i l_i t_i t_i^T / sum_i l_i, which describe how much the bundle’s directions spread there.

Parameters:

phantom (Phantom) – The phantom.
seed (int | None) – The seed for the bundle geometries that have none of their own; each of them gets a seed drawn from this one, bundle by bundle. Defaults to None.
n_jobs (int) – Threads for resampling the streamlines, -1 for all cores. Defaults to -1.
progress (Callable[[int, int, str], None] | None) – Called with (done, total, “voxelise”) after each bundle. Defaults to None.

Returns:

The partition.

Return type:

VoxelPartition

larmor.simulation.read_phantom(path)#

Reads a phantom file written by write_phantom.

The tractograms and tissue maps the file refers to are loaded too, from their path relative to the file if it exists, else from their absolute path, so a folder of files can be moved together.

Parameters:

path (str | Path) – The file.

Returns:

The phantom, the protocol if the file has one, and the viewer block.

Raises:

ValueError – If the file is not a phantom file, was written by a newer Larmor, or holds invalid values.
FileNotFoundError – If a tractogram or tissue map cannot be found.

Return type:

PhantomFile

larmor.simulation.sandi_tissue(*, name='GM', neurite_fraction=0.5, soma_fraction=0.15, neurite_diffusivity=0.0023, soma_diffusivity=0.003, soma_radius_um=6.0, extracellular_diffusivity=0.0015, t1_ms=1331.0, t2_ms=88.0, t2star_ms=46.7, proton_density=0.8)#

SANDI grey matter: a sphere (soma), a powder-averaged stick and a ball.

The neurites are sticks averaged over all directions, the cell bodies (somas) are impermeable spheres, and the extracellular space is a ball with the fraction 1 - neurite_fraction - soma_fraction. The signal therefore does not depend on the gradient direction, and the tissue needs no fibre direction. The three compartments share T1, T2, T2* and proton density, as SANDI itself has no relaxation. The relaxation defaults assume 3 T. With a gradient-echo protocol, the sphere sees the bipolar gradient pair that Protocol describes.

Sources of the defaults:

soma_diffusivity 3.0e-3 mm^2/s: the value SANDI fixes (Palombo et al. 2020).
soma_radius_um 6, the fractions 0.5 (neurite), 0.15 (soma) and 0.35 (extracellular), extracellular_diffusivity 1.5e-3 mm^2/s: choices; the paper gives no single default.
neurite_diffusivity 2.3e-3 mm^2/s: a choice, the same as the white matter stick, close to the in vivo intra-axonal diffusivity of 2.25 um^2/ms measured by Dhital et al. (2019).
T1 1331 ms: grey matter at 3 T in vivo (Wansapura et al. 1999).
T2 88 ms: frontal grey matter at 3 T in vivo (Lu et al. 2005, Table 1).
T2* 46.7 ms: the mean of the frontal and occipital grey matter values at 3 T (Wansapura et al. 1999).
proton_density 0.8: an approximate choice for grey matter water content (source to be added).

Parameters:

name (str) – The tissue’s name. Defaults to “GM”.
neurite_fraction (float) – The powder-averaged stick fraction. Defaults to 0.5.
soma_fraction (float) – The sphere fraction. Defaults to 0.15.
neurite_diffusivity (float) – The stick’s D, mm^2/s. Defaults to 2.3e-3.
soma_diffusivity (float) – The D of the water inside the sphere, mm^2/s. Defaults to 3.0e-3.
soma_radius_um (float) – The sphere radius, um. Defaults to 6.
extracellular_diffusivity (float) – The ball’s D, mm^2/s. Defaults to 1.5e-3.
t1_ms (float) – T1 of every compartment, ms. Defaults to 1331.
t2_ms (float) – T2 of every compartment, ms. Defaults to 88.
t2star_ms (float) – T2* of every compartment, ms. Defaults to 46.7.
proton_density (float) – The proton density of every compartment. Defaults to 0.8.

Returns:

The tissue.

Raises:

ValueError – If the fractions are negative or exceed 1 together.

Return type:

Tissue

References

Palombo M, Ianus A, Guerreri M, et al. SANDI: a compartment-based model for non-invasive apparent soma and neurite imaging by diffusion MRI. NeuroImage 215 (2020) 116835. Dhital B, Reisert M, Kellner E, Kiselev VG. Intra-axonal diffusivity in brain white matter. NeuroImage 189 (2019) 543-550. Wansapura JP, Holland SK, Dunn RS, Ball WS. NMR relaxation times in the human brain at 3.0 tesla. J Magn Reson Imaging 9 (1999) 531-538. Lu H, Nagae-Poetscher LM, Golay X, Lin D, Pomper M, van Zijl PCM. Routine clinical brain MRI sequences for use at 3.0 Tesla. J Magn Reson Imaging 22 (2005) 13-22. doi:10.1002/jmri.20356

larmor.simulation.save_simulation(simulated, folder, *, force=False, bundle_names=None)#

Writes a simulation and its ground truth into one self-contained folder.

The folder holds:

the DWI as NIfTI, one file per echo, with .bval and .bvec in the FSL frame;
the ground truth: the tissue, bundle, compartment and free-water volume fractions, the compartments’ and free water’s signal fractions per echo, the fibre directions with their amplitudes, and one mask per tissue;
the bundles’ streamlines as .tck files;
the phantom file with the scan;
simulation.json, with the noise, the seed and the name of every file.

Running the phantom file again with the same seed gives the same images. When no seed was given, the seed that simulate drew is in simulation.json. The folder is created if needed. If it holds an earlier simulation, the files that its simulation.json lists are removed first, so no file of the earlier run stays beside the new ones. A folder that holds other files is refused unless force is True, which then writes into it and keeps them. The paths are those of SimulationFiles.

The ground-truth fibre directions come from the voxel’s own fibres, not from a fit. The fibres are taken from the largest down, and each joins the first peak within 20 degrees of it; a direction and its opposite count as the same. A peak’s direction is the weighted mean of its fibres, and its amplitude is their summed share of the voxel. Peaks below 10% of the voxel’s largest are dropped. A tissue mask holds the voxels where the tissue, background included, fills at least half.

The tissue, bundle, compartment and free-water fractions are volume fractions, the share of the voxel’s volume that each fills. A model fitted to the images recovers signal fractions, which differ, because each compartment’s signal is weighted by its proton density and relaxation: at TE 80 ms and TR 3000 ms, the NODDI preset’s free-water volume fraction of 0.05 is about 0.11 of the b = 0 signal. The signal fractions are written too, one file per echo time (see SimulationFiles); compare a fit with those.

Parameters:

simulated (SimulatedDWI) – The result of simulate.
folder (str | Path) – The output folder.
force (bool) – Whether to write into a folder that is not empty and holds no simulation.json. Defaults to False.
bundle_names (Sequence[str] | None) – A name for each bundle, kept in the phantom file for the viewer, which shows the bundles by these names (for example PhantomFile.bundle_names of the file the phantom was read from). Defaults to None: the viewer numbers them.

Returns:

The paths written.

Raises:

FileExistsError – If the folder is not empty, has no simulation.json and force is False.
ValueError – If simulated does not carry its phantom and protocol (it was not made by simulate), two tissue names give the same file name, or bundle_names does not name every bundle.

Return type:

SimulationFiles

larmor.simulation.signal_breakdown(phantom, protocol, partition, voxel, volume, *, echo=0, snr=None, noise_reference='WM', bundle_names=None)#

Splits one voxel’s noise-free signal in one volume into its factors.

Use it to see where a simulated value comes from. The signal is the one simulate computes for the voxel, from the same compartment models. For each compartment of each tissue in the voxel, it is the product of the tissue’s partial volume, the compartment’s fraction and proton density, the relaxation factors 1 - exp(-TR / T1) and exp(-TE / T2) (with T2* for a gradient echo), and the diffusion attenuation. A tissue that fills bundles is split by bundle, since each bundle has its own fibre directions. Only the one voxel is evaluated, so the call is fast enough to follow the mouse in the viewer.

With snr, the noise is that of simulate, sigma = S0_ref / snr, and the result also gives the noise floor and the mean of the Rician magnitude.

Parameters:

phantom (Phantom) – The phantom.
protocol (Protocol) – The scan.
partition (VoxelPartition) – The phantom’s voxel partition, as partition_voxels or simulate returns it.
voxel (Sequence[int]) – The voxel (i, j, k).
volume (int) – The volume’s index in the protocol’s gradient table.
echo (int) – The echo’s index in the protocol’s echo times. Defaults to 0.
snr (float | None) – The signal-to-noise ratio of S0_ref, as simulate takes it. Defaults to None: no noise.
noise_reference (str) – The name of the tissue that defines S0_ref. Defaults to “WM”.
bundle_names (Sequence[str] | None) – A name for each bundle of the phantom. Defaults to None: “bundle 1”, “bundle 2” and so on.

Returns:

The breakdown.

Raises:

ValueError – If the voxel is outside the grid, the volume or echo does not exist, the partition is not the phantom’s, snr is not positive, or the noise reference is not a tissue of the phantom.

Return type:

SignalBreakdown

References

Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson Med 34 (1995) 910-914.

larmor.simulation.simulate(phantom, protocol, *, snr=None, noise_reference='WM', seed=None, n_jobs=-1, progress=None)#

Simulates the diffusion-weighted images of a phantom, voxel by voxel.

This is the instant engine: it sums the closed-form signal of each compartment in each voxel, as if the whole acquisition happened at once, so there is no k-space, readout, distortion or motion. The signal of a voxel at echo e and gradient n is S = sum_c f_c A_c(b_n, g_n) w_c(e), a sum over the compartments c of the tissues in the voxel (see partition_voxels) of three factors:

f_c, the compartment’s fraction of the voxel: the tissue’s fraction in the voxel times the compartment’s fraction in the tissue.
A_c, the compartment’s attenuation (see Compartment). A compartment that needs a fibre direction is evaluated along each fibre direction in the voxel, weighted by the fibre’s share.
w_c, the relaxation weight PD (1 - exp(-TR / T1)) exp(-TE / T2), with T2* in place of T2 for a gradient echo (whose diffusion weighting is the bipolar pair that Protocol describes).

The signal is computed on n_jobs threads and does not depend on their number. With snr, noise is added by add_noise with sigma = S0_ref / snr, where S0_ref is the b = 0 signal, at the first echo time, of a voxel filled with the noise_reference tissue.

Parameters:

phantom (Phantom) – The object to simulate.
protocol (Protocol) – The scan.
snr (float | None) – The signal-to-noise ratio of S0_ref. Defaults to None: no noise.
noise_reference (str) – The name of the tissue that defines S0_ref. Defaults to “WM”.
seed (int | None) – The seed of the noise and of the bundle geometries that have no seed of their own (see partition_voxels); the same seed gives the same result. Defaults to None: a seed is drawn and returned in the result.
n_jobs (int) – The number of threads, -1 for all cores. Defaults to -1.
progress (Callable[[int, int, str], None] | None) – A function called with (done, total, stage) as the run goes on, where stage is “voxelise”, “signal” or “noise”; raise an exception in it to cancel the run. Defaults to None.

Returns:

The images, the noise-free images and the ground truth.

Raises:

ValueError – If snr is not positive, the noise reference is not a tissue of the phantom, or n_jobs is 0 or below -1.

Return type:

SimulatedDWI

larmor.simulation.simulation_files(phantom, protocol, folder, *, noise)#

The files a simulation of this phantom and protocol is saved as.

Use it to know the paths save_simulation will write before running a simulation.

Parameters:

phantom (Phantom) – The phantom.
protocol (Protocol) – The scan; its echo times name the DWI files.
folder (str | Path) – The output folder.
noise (bool) – Whether the simulation has noise, which adds the noise-free DWI.

Returns:

The paths; nothing is written.

Raises:

ValueError – If two tissues have names that give the same file name.

Return type:

SimulationFiles

larmor.simulation.standard_model_tissue(*, name='WM', intra_fraction=0.7, free_water_fraction=0.0, intra_diffusivity=0.0023, extra_parallel_diffusivity=0.0017, extra_perpendicular_diffusivity=0.0005, free_water_diffusivity=0.003, intra_t2_ms=80.0, extra_t2_ms=45.0, t1_ms=832.0, t2star_ms=46.55, proton_density=0.7)#

The standard model of white matter: stick, zeppelin and optional free water.

The water inside the axons is a stick, and the water in the space around them is a zeppelin. Free water is an optional FREE_WATER compartment with the parameters of csf_tissue apart from T2*. The extra-axonal fraction is 1 - intra_fraction - free_water_fraction. The relaxation defaults assume 3 T.

t2star_ms is the tissue’s T2*. Reversible dephasing (T2’) is a property of the field in the voxel, not of one compartment, so the stick and the zeppelin share one T2’ = 1 / (1/T2* - 1/T2_mean). T2_mean is the mean T2 of the stick and the zeppelin, weighted by their signal at TE = 0 (fraction times proton density). Each compartment c, free water included, then gets T2*_c = 1 / (1/T2_c + 1/T2’).

Sources of the defaults:

intra_diffusivity 2.3e-3 mm^2/s: a choice close to the in vivo intra-axonal diffusivity of 2.25 um^2/ms measured by Dhital et al. (2019); NODDI fixes 1.7e-3.
extra-axonal diffusivities 1.7e-3 and 0.5e-3 mm^2/s: they follow the NODDI constraints listed by Novikov et al. (2019, Sec. 3.4.2), D_e_par = D_a = 1.7 um^2/ms and the tortuosity D_e_perp = (1 - f) D_e_par, about 0.5 um^2/ms at f = 0.7.
intra_fraction 0.7: a choice.
T2 of 80 ms (intra) and 45 ms (extra): choices around the T2 of 69 ms of frontal white matter at 3 T in vivo (Lu et al. 2005, Table 1); that the intra-axonal T2 is the longer agrees with Veraart et al. (2018).
T1 832 ms: white matter at 3 T in vivo (Wansapura et al. 1999).
t2star_ms 46.6 ms: white matter at 3 T in vivo, the mean of the frontal and occipital values of Wansapura et al. (1999).
proton_density 0.7: an approximate choice for white matter water content (source to be added).

Parameters:

name (str) – The tissue’s name. Defaults to “WM”.
intra_fraction (float) – The stick fraction. Defaults to 0.7.
free_water_fraction (float) – The free-water fraction. Defaults to 0.0 (no free water).
intra_diffusivity (float) – The stick’s D, mm^2/s. Defaults to 2.3e-3.
extra_parallel_diffusivity (float) – The zeppelin’s D along the fibre, mm^2/s. Defaults to 1.7e-3.
extra_perpendicular_diffusivity (float) – The zeppelin’s D across the fibre, mm^2/s. Defaults to 0.5e-3.
free_water_diffusivity (float) – The free water’s D, mm^2/s. Defaults to 3.0e-3.
intra_t2_ms (float) – The stick’s T2, ms. Defaults to 80.
extra_t2_ms (float) – The zeppelin’s T2, ms. Defaults to 45.
t1_ms (float) – T1 of the stick and the zeppelin, ms. Defaults to 832.
t2star_ms (float) – The tissue’s T2*, ms, at most T2_mean. Defaults to 46.6.
proton_density (float) – The proton density of the stick and the zeppelin. Defaults to 0.7.

Returns:

The tissue.

Raises:

ValueError – If the fractions are negative or exceed 1 together, or t2star_ms is not positive or exceeds T2_mean.

Return type:

Tissue

References

Novikov DS, Fieremans E, Jespersen SN, Kiselev VG. Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation. NMR Biomed 32 (2019) e3998. Dhital B, Reisert M, Kellner E, Kiselev VG. Intra-axonal diffusivity in brain white matter. NeuroImage 189 (2019) 543-550. Veraart J, Novikov DS, Fieremans E. TE dependent Diffusion Imaging (TEdDI) distinguishes between compartmental T2 relaxation times. NeuroImage 182 (2018) 360-369. Wansapura JP, Holland SK, Dunn RS, Ball WS. NMR relaxation times in the human brain at 3.0 tesla. J Magn Reson Imaging 9 (1999) 531-538. Lu H, Nagae-Poetscher LM, Golay X, Lin D, Pomper M, van Zijl PCM. Routine clinical brain MRI sequences for use at 3.0 Tesla. J Magn Reson Imaging 22 (2005) 13-22. doi:10.1002/jmri.20356

larmor.simulation.write_phantom(phantom, path, *, protocol=None, viewer=None)#

Writes a phantom file (JSON, suggested suffix .phantom.json).

A phantom file stores a phantom, and optionally its scan, so that read_phantom, lar simulate and the viewer’s Phantom panel can open it. The file holds the format name “larmor-phantom” and version 1, the grid, the tissues, the bundles (each a geometry or a tractogram’s path), the tissue maps (their paths) and the background. An optional “scan” block holds the protocol, and an optional “viewer” block holds the viewer’s settings as given. Each path is stored both as an absolute path and relative to the file’s folder. A bundle geometry without a seed gets one drawn at random and written, so the file always gives the same streamlines.

Parameters:

phantom (Phantom) – The phantom.
path (str | Path) – The file to write.
protocol (Protocol | None) – The scan to store with it. Defaults to None.
viewer (Mapping[str, Any] | None) – The viewer’s settings, any JSON-ready values. Defaults to None.

Raises:

ValueError – If a tractogram bundle or a tissue map has no path.