Isochromats

class pulserver.virtual.Isochromats

Bases: object

Isochromats, and the magnetisation the Bloch equation carries from one block to the next.

The magnetisation turns about the field b = (Re b1, Im b1, bz), in Hz, in the frame rotating at the reference frequency, as the magnetic moment of a nucleus of positive gyromagnetic ratio precesses: dM/dt = 2 pi M x b, clockwise seen from the tip of b. A 90 degree pulse along +x takes +z to +y, and an isochromat at a positive off-resonance or along a positive gradient accrues Mx + i My a negative phase.

Free precession under a piecewise-linear gradient, with relaxation, is integrated exactly and applied when the magnetisation is next needed. An RF pulse is a sequence of steps, each a rotation about the step’s mean field between two half steps of relaxation; isochromats that see the same field during a pulse share its computation, and a pulse that differs from an earlier one by its phase alone, under the same gradient, reuses the earlier one’s computation turned about z by that phase. A pulse played under no gradient or one held throughout it is computed on a grid of the field an isochromat sees and interpolated, to within about 1e-7 of the equilibrium magnetisation, wherever the grid costs fewer maps than the isochromats’ groups: without transmit sensitivities, or with every channel playing one waveform times a weight of its own, when an isochromat’s transmit field is that waveform times one complex drive, whose magnitude the grid spans too and whose phase turns the computation about z. An ADC window under a gradient held throughout it is read by a non-uniform FFT of the isochromats of each T2, to within about 1e-13 of the sum of the magnitudes of their transverse magnetisations times their receive sensitivities, wherever that costs less than turning every isochromat at every sample. A window under any other gradient whose k moves along axes on which the isochromats lie on a lattice is read by FINUFFT, wherever that costs less: each isochromat’s value is summed onto its lattice point for each coil and each of a few Chebyshev points across the window, between which its decay and precession are interpolated, and the lattice is transformed to each sample’s k, to within about 1e-11 of the sum of the magnitudes of the terms each sample sums. play() reads either to within a tolerance instead where given one.

Parameters:

• positions (array_like) – (n, 3) positions, in m, along the axes the gradients are played on.

• proton_density (float or array_like, default=1.0) – Equilibrium longitudinal magnetisation, per isochromat.

• t1 (float or array_like, default=inf) – Relaxation times, in s; inf for none.

• t2 (float or array_like, default=inf) – Relaxation times, in s; inf for none.

• off_resonance (float or array_like, default=0.0) – Precession frequency at rest, in Hz from the reference frequency: field inhomogeneity and chemical shift together.

• transmit (array_like, default=None) – Complex transmit sensitivities, (n,) or (n, channels), scaling the field of each RF channel. By default one channel of unit sensitivity, onto which the channels of a pTx pulse are summed.

• receive (array_like, default=None) – Complex receive sensitivities, (n,) or (n, coils). By default one coil of unit sensitivity. A C-contiguous complex128 array, a memory-mapped one included, is read in place into the isochromats’ own layout, so that it is never copied whole into memory.

• diffusion (float or array_like, default=0.0) – Isotropic diffusion coefficient, in m²/s, per isochromat.

• motion (callable, default=None) – Where the isochromats lie as their subject moves: motion(t, positions) returns the (n, 3) positions, in m, that the positions at rest move to at time t, in s on the isochromats’ clock (elapsed); RigidMotion for a rigid body. Every isochromat keeps its other properties as it moves.

• seed (int or numpy.random.Generator, default=None) – Seed, or generator, of the Brownian walks that diffusion drives.

• threads (int, default=0) – Worker threads; 0 for every core.

• device (str or torch.device, default=None) – A torch device the ADC windows are read and the runs of repetitions() carried on. A window outside a run whose isochromats lie on a lattice is summed onto it by a Triton kernel and transformed by cuFINUFFT, and every other one summed sample by sample by a second Triton kernel; a run is carried and spread onto its windows’ grids by more. A CUDA device, or the CPU under Triton’s interpreter alone (TRITON_INTERPRET=1 before triton is first imported), on which FINUFFT transforms the lattice. Needs the gpu extra. By default the engine does all of it itself.

Raises:

• ValueError – If a property has the wrong shape or is not finite, a relaxation time is not positive, or a diffusion coefficient is negative.

• ImportError – If device is given without torch and Triton, or is a CUDA device without cuFINUFFT.

• RuntimeError – If device is the CPU outside Triton’s interpreter.

Notes

Every isochromat starts at equilibrium, along +z. The engine computes what the sequence does to these isochromats; it does not model the scanner’s hardware.

Isochromats that move or diffuse play block by block. Before each block, motion places them where they are at its start, and they stay there through the block; after it, each is turned by the phase the block’s gradient adds along its path within the block, from the end of the block’s RF pulse on, or from the start of a block without one. A diffusing isochromat keeps its position and follows a Brownian walk of its own, which turns it by the phase the gradient plays along the walk: from the walk so far, and from its increments over the block, drawn with the displacement they make, from their joint normal distribution. A voxel’s diffusion attenuation is the mean of its isochromats’ phase factors, which holds to about one over the square root of their number. Samples within a block are read with every isochromat where it stood at the block’s start.

Examples

>>> import numpy as np
>>> from pulserver.virtual import Isochromats
>>> spins = Isochromats([[0.0, 0.0, 0.0]], t2=0.1)

A 25 kHz field along +x held for 10 us turns +z by 90 degrees, to +y; the transverse magnetisation then decays with t2:

>>> pulse = (0.0, 10e-6, [25e3])
>>> signal = spins.play(10.01e-3, rf=pulse, adc=[10e-6, 10.01e-3])
>>> np.round(signal, 4)
array([[0.+1.j    , 0.+0.9048j]])

Methods

play	Play one block's events and return what each coil receives at each ADC sample.
repetitions	Repetitions of a sequence of blocks, from the magnetisation where it stands.
reset	Return every isochromat to equilibrium, along +z, the clock to zero, and a diffusing one to the start of a new walk.

Attributes

coils	Number of receive coils; one without receive sensitivities.
device_windows	ADC windows read on the device since construction, on a lattice or sample by sample.
elapsed	Time played since construction or the last reset(), in s.
lattice_windows	ADC windows read on a lattice by FINUFFT since construction.
magnetization	(n, 3) magnetisation, a copy; assigning replaces it.
moving	Whether the isochromats move or diffuse, and so play block by block.
positions	(n, 3) positions, in m, as the last block played them, a copy.
ungrouped_pulses	RF pulses played isochromat by isochromat from their tables since construction.