diff --git a/.gitignore b/.gitignore index c2d15588c5c91608146516ca8fc1b9b357776d69..fbc8936473fe27ba37cc4ee2d4709401cbd35776 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,8 @@ docs/source/example_gallery/**/*.ipynb build/ dist/ tests/test_output/ + + +pylint/ +htmlcov/ +.coverage \ No newline at end of file diff --git a/README.rst b/README.rst index 086f06ff165acd60b6906617b55c8f6c24fcaaca..d71064633c58c137806d2a4428c4c3a136b39bf1 100644 --- a/README.rst +++ b/README.rst @@ -1,4 +1,4 @@ -.. image:: docs/source/_static/Logo_cmrsim.svg +.. image:: https://people.ee.ethz.ch/~jweine/cmrsim/latest/_images/Logo_cmrsim_moving_light.svg :align: right :scale: 150 % @@ -22,14 +22,15 @@ What makes CMRsim special? - Reproducibility by using versioned and documented open-source software. - Easily scalable to detailed digital phantoms + Installation ============== - Pip ^^^^ CMRsim can be installed using the python package manager :code:`pip`. The package can be found -in the `registry`_ of the repository hosted on the ETH-Gitlab instance +in the `registry`_ of the repository hosted on the ETH-Gitlab instance as well as on the python +package index (PyPI). .. _registry: https://gitlab.ethz.ch/jweine/cmrsim/-/packages @@ -38,30 +39,17 @@ Stable release versions are denoted with the ususal :code:`major.minor` notation development versions are constantly updated on solved issued and can be installed using the :code:`major.minor.devX` notation -Installation in command line: +Release versions are published on PyPI, which allows installing cmrsim with: .. code-block:: - pip install cmrsim --extra-index-url https://__token__:<your_personal_token>@gitlab.ethz.ch/api/v4/projects/23798/packages/pypi/simple - + pip install cmrsim -Requirements -^^^^^^^^^^^^^^^^ -CMRsim does rely on functionality of 3 non-standard python stack packages: +Development versions are only available using the extra url: -1. `TensorFlow`_: Pretty much all computationally heavy calculations. -2. `Pyvista`_: For processing and visualizing graphs (e.g. 3D Velocity fields). -3. `CMRseq`_: For defining MR-sequences used for Fourier-Encoding calculation and Bloch-simulation. - -.. _TensorFlow: https://www.tensorflow.org/ -.. _Pyvista: https://docs.pyvista.org/ -.. _CMRseq: https://people.ee.ethz.ch/~jweine/cmrseq/index.html - -Further (already indirectly included) requirements are: --h5py~=2.10.0 --scipy>=1.4.1 --tqdm +.. code-block:: + pip install cmrsim --extra-index-url https://gitlab.ethz.ch/api/v4/projects/23798/packages/pypi/simple Docker ^^^^^^^^ diff --git a/ci_cd_config/.run-test-stage-ci.yml b/ci_cd_config/.run-test-stage-ci.yml index 3712ad07518560c083f3e459a8280411c4116566..b6eb17b041818a86a6be1f981dbfd8a843dd0db8 100644 --- a/ci_cd_config/.run-test-stage-ci.yml +++ b/ci_cd_config/.run-test-stage-ci.yml @@ -5,11 +5,17 @@ run_unit_test: image: registry.ethz.ch/jweine/cmrsim:unittest stage: test + variables: + BATCH_TAG: "" + rules: + - if: $CI_COMMIT_BRANCH + variables: + BATCH_TAG: ":dev" script: - coverage run --source cmrsim -m unittest discover tests - coverage report - coverage xml - - genbadge coverage -i coverage.xml -o coverage.svg + - genbadge coverage -i coverage.xml -o coverage.svg --name "coverage${BATCH_TAG}" artifacts: paths: - coverage.svg diff --git a/cmrsim/analytic/contrast/base.py b/cmrsim/analytic/contrast/base.py index f16c8c0ca6837c2bdce1ac76b75ae0636296c7cb..4e6b59fe1fca9e920f3ba67ece6583fca0363bbf 100644 --- a/cmrsim/analytic/contrast/base.py +++ b/cmrsim/analytic/contrast/base.py @@ -250,7 +250,7 @@ class LookUpTableModule(BaseSignalModel): :param r_vectors: (-1, 3) in meter :param batch_index: (#batch) indices of the first dimension (T) used for spatial lookup - :return: (#batch_idx, -1, self._n_channels) trilinear interpolated map values for specifiec + :return: (#batch_idx, -1, self._n_channels) trilinear interpolated map values for specified locations """ with tf.name_scope("linear_interpolation_lookup") as scope: diff --git a/cmrsim/analytic/contrast/coil_sensitivities.py b/cmrsim/analytic/contrast/coil_sensitivities.py index 8a0ec0d393f7397b69cdb415a41b70953f3ffe73..162bf1f1c483963a8a19371a51781e239af6050d 100644 --- a/cmrsim/analytic/contrast/coil_sensitivities.py +++ b/cmrsim/analytic/contrast/coil_sensitivities.py @@ -116,6 +116,21 @@ class CoilSensitivity(LookUpTableModule): result = tf.reshape(product, new_shape) return result + def lookup_complex_coil_factors(self, r_vectors: tf.Tensor): + """ Looks the sensitivity values for all coils at specified locations. + + :param r_vectors: [..., 3] dtype: tf.float32 + :return: coil_factors [..., n_coils] dtype: tf.complex64 + """ + old_shape = tf.shape(r_vectors) + new_shape = tf.concat([old_shape[:-1], [self.expansion_factor, ]], axis=0) + r_flat = tf.reshape(r_vectors, [-1, 3]) + sense_factors = self.linear_interpolation(r_vectors=r_flat) + sense_factors = tf.reshape(sense_factors, [-1, self.expansion_factor, 2]) + sense_factors = tf.complex(sense_factors[..., 0], sense_factors[..., 1]) + sense_factors = tf.reshape(sense_factors, new_shape) + return sense_factors + @property # pylint: disable=invalid-name def coil_maps(self): @@ -137,6 +152,46 @@ class CoilSensitivity(LookUpTableModule): coils = tf.reshape(coils, [*coil_sensitivities.shape[1:-1], -1]) self.look_up_map3d.assign(coils) + @classmethod + def from_3d_birdcage(cls, map_shape: Union['np.ndarray', tf.Tensor, Tuple], + map_dimensions: np.ndarray, relative_coil_radius: float = 1.5, + device: str = 'CPU:0', **kwargs): + """Lookup for ideal bird-cage coil sensitivities. + + .. note:: + + expand_repetitions is always set to True for coil-sensitivities, which means that + the output hast a #repetions * #coils dimension on the repetitions axis + + :param map_shape: (X, Y, Z, #coils) array shape of the resulting coil-maps + :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)] Extend of + :param relative_coil_radius: distance of the birdcage ring from the map center + the map in scanner coordinate system. Used to sort iso-chromat positions into bins. + This is assumed to be the *length* in meter in each spatial direction + :param device: (str) Device on which all operations except for the look-up are + placed (default: GPU:0) + :param device_lookup: (str) Device where the lookup table is places (defaults to CPU:0) + + """ + nx, ny, nz, nc = map_shape + c, z, y, x = np.mgrid[:nc, :nz, :ny, :nx] + coilx = relative_coil_radius * np.cos(c * (2 * np.pi / nc)) + coily = relative_coil_radius * np.sin(c * (2 * np.pi / nc)) + coilz = np.floor(c / nc) - 0.25 # * nc + coil_phs = -(c + np.floor(c / nc)) * (2 * np.pi / nc) + + x_co = (x - nx / 2.0) / (nx / 2.0) - coilx + y_co = (y - ny / 2.0) / (ny / 2.0) - coily + z_co = (z - nz / 2.0) / (nz / 2.0) - coilz + rr = (x_co ** 2 + y_co ** 2 + z_co ** 2) ** 0.5 + phi = np.arctan2(x_co, -y_co) + coil_phs + out = (1 / rr) * np.exp(1j * phi) + + rss = sum(abs(out) ** 2, 0) ** 0.5 + out /= rss + coil_maps3d = tf.constant(out.transpose(3, 2, 1, 0), dtype=tf.complex64) + return cls(coil_maps3d, map_dimensions, device, **kwargs) + class BirdcageCoilSensitivities(CoilSensitivity): """ Lookup for ideal bird-cage coil sensitivities. diff --git a/cmrsim/analytic/contrast/phase_tracking.py b/cmrsim/analytic/contrast/phase_tracking.py index adcd8938edfeb75ad780dcf80360691154e13dd2..81b704a5aa078f50b63d6d84dbb8b026e6bfe08b 100644 --- a/cmrsim/analytic/contrast/phase_tracking.py +++ b/cmrsim/analytic/contrast/phase_tracking.py @@ -78,7 +78,8 @@ class PhaseTracking(BaseSignalModel): # Case 2: repetition-axis of signal_tensor must match self.expansion_factor else: tf.Assert(input_shape[1] == self.expansion_factor, - ["Shape missmatch for input argument signal_tensor for case no-expand"]) + ["Shape missmatch for input argument signal_tensor for case no-expand", + input_shape, self.expansion_factor]) result = signal_tensor * complex_phase return result diff --git a/cmrsim/analytic/contrast/sequences.py b/cmrsim/analytic/contrast/sequences.py index c10e6fa2fccf1e45f4bea7f13a7b04f5c9ccc7ef..3720180e0f2fa9e7be468fe0f980d5887ab9759f 100644 --- a/cmrsim/analytic/contrast/sequences.py +++ b/cmrsim/analytic/contrast/sequences.py @@ -124,7 +124,7 @@ class BSSFP(BaseSignalModel): :param expand_repetitions: See BaseSignalModel for explanation """ # pylint: disable=invalid-name - required_quantities = ('T1', 'T2', 'T2star', 'theta') + required_quantities = ('T1', 'T2', 'T2star', 'off_res') #: Echo time in ms TE: tf.Variable = None #: Repetition time in ms @@ -161,7 +161,7 @@ class BSSFP(BaseSignalModel): # pylint: disable=too-many-arguments def __call__(self, signal_tensor: tf.Tensor, T1: tf.Tensor, T2: tf.Tensor, - T2star: tf.Tensor, theta: tf.Tensor, **kwargs): + T2star: tf.Tensor, off_res: tf.Tensor, **kwargs): # pylint: disable=invalid-name, too-many-locals, """ Evaluates the signal Equation. @@ -169,10 +169,13 @@ class BSSFP(BaseSignalModel): :param T1: in ms :param T2: in ms :param T2star: in ms - :param theta: phase in radians acquired over one TR due to off-resonance + :param off_res: b0 inhomogeneity in rad/ms (angular frequency difference) :return: (#voxel, #repetitions, #k-space-samples) """ + # phase in radians acquired over one TR due to off-resonance + theta = off_res * tf.reshape(self.TR, [1, -1, 1]) + T2p = 1. / (1. / T2star - 1./T2) E1 = tf.exp(-tf.reshape(self.TR, [1, -1, 1]) / T1) E2 = tf.exp(-tf.reshape(self.TR, [1, -1, 1]) / T2) @@ -183,7 +186,7 @@ class BSSFP(BaseSignalModel): b = tf.math.cos(tf.reshape(self.FA, [1, -1, 1])) - E1 # Equation 5 - Lambda = (a - E2**2 * b - tf.sqrt((a**2 - E2**2 * b**2) * (1 - E2**2))) / (a - b) + Lambda = (a - E2**2 * b - tf.sqrt((a**2 - E2**2 * b**2) * (1 - E2**2))) / (a - b) # Equation 3 A = ((a + b) * E2 / 2.) / (a - Lambda * (a - b) / 2.) diff --git a/cmrsim/analytic/contrast/t2_star.py b/cmrsim/analytic/contrast/t2_star.py index 9122e87e011d373e2afe676d721061eeb97aef93..72d5507852c5deda3cc350d7d1a338b90c32f083 100644 --- a/cmrsim/analytic/contrast/t2_star.py +++ b/cmrsim/analytic/contrast/t2_star.py @@ -94,7 +94,7 @@ class StaticT2starDecay(BaseSignalModel): if self.expand_repetitions or self.expansion_factor == 1: result = signal_tensor[:, :, tf.newaxis, :] * decay_factor[:, tf.newaxis, :, :] - new_shape = tf.concat([input_shape[0], input_shape[1] * self.expansion_factor, + new_shape = tf.stack([input_shape[0], input_shape[1] * self.expansion_factor, n_samples], axis=0) result = tf.reshape(result, new_shape) else: diff --git a/cmrsim/analytic/simulation.py b/cmrsim/analytic/simulation.py index a3de5877f3acc6d0695a9a324b47ba1c21e83cb0..d046de0710bcd3aca05fe369cf8f6a42e45cda43 100644 --- a/cmrsim/analytic/simulation.py +++ b/cmrsim/analytic/simulation.py @@ -162,6 +162,8 @@ class AnalyticSimulation(tf.Module): s_of_k_segments = tf.TensorArray(dtype=tf.complex64, size=n_segments, dynamic_size=False) self.progress_bar.reset_segments() + m_transverse = self.forward_model(batch_dict["M0"], segment_index=0, **batch_dict) + # Loop over k-space segments, defined by k-space trajectory in self.encoding_module m_transverse = self.forward_model(batch_dict["M0"], segment_index=0, **batch_dict) for segment_index in tf.range(n_segments): diff --git a/cmrsim/bloch/__init__.py b/cmrsim/bloch/__init__.py index 7b4c1ad5639b6a04b0318fe4e2db4176f5bcd540..949228d9fa806dda8cffee6aca717c8062532615 100644 --- a/cmrsim/bloch/__init__.py +++ b/cmrsim/bloch/__init__.py @@ -3,10 +3,11 @@ from cmrsim.bloch._base import * from cmrsim.bloch._generic import * from cmrsim.bloch._ideal import * +from cmrsim.bloch._multi_coil import * import cmrsim.bloch.submodules import importlib -_submodules = ["_base", "_generic", "_ideal", "_ideal"] +_submodules = ["_base", "_generic", "_ideal", "_multi_coil"] mod_handles = [importlib.import_module(f"cmrsim.bloch.{m}") for m in _submodules] __all__ = [item for m in mod_handles for item in getattr(m, '__all__')] + ["submodules"] diff --git a/cmrsim/bloch/_generic.py b/cmrsim/bloch/_generic.py index e7ff204d9b25c8ef8cf1ce894279456ae3b95481..6e7494560b7b67bd5c331a51fdb5d3707c8b45cb 100644 --- a/cmrsim/bloch/_generic.py +++ b/cmrsim/bloch/_generic.py @@ -328,118 +328,158 @@ class GeneralBlochOperator(BaseBlochOperator): :param n_repetitions: actual n_reps for par, 1 for sequential :param n_samples: :param trajectory_module: - :param initial_position: (n_repetitions/1, #batch, 3) + :param initial_position: (#batch, 3) :param magnetization: (n_repetitions/1, #batch, 3) :param T1: (n_repetitions/1, #batch) :param T2: (n_repetitions/1, #batch) :param M0: (n_repetitions/1, #batch) :param kwargs: - :return: magnetization, r_next, signal at adc_on==True + :return: - magnetization + - r_next + - signal at adc_on==True with shape (n_repetitions, n_samples) + """ with tf.device(self.device): - batch_signal = tf.TensorArray(dtype=tf.complex64, size=n_samples + tf.constant(1), - element_shape=(n_repetitions,)) - result_writing_idx = 0 - r_next = initial_position - repetitions = tf.range(repetition_index, repetition_index + n_repetitions, 1) - - # If gradient waveforms are specified, gather all required repetitions (either all or 1) - # and calculate the left value of for trapezoidal phase integration - if self.gradient_waveforms is not None: - current_gradient_waveforms = tf.gather(self.gradient_waveforms, repetitions, - axis=0)[:, tf.newaxis] - gdotr_l = tf.reduce_sum(current_gradient_waveforms[..., 0, :] * r_next, axis=-1) - else: - gdotr_l = tf.constant(0, dtype=tf.float32) - gdotr_r = tf.zeros_like(gdotr_l) - - # If rf waveforms are specified, gather all required repetitions (either all or 1) - if self.rf_waveforms is not None: - current_rf_waveforms = tf.gather(self.rf_waveforms, repetitions, axis=0) - - # If adc events are specified, gather all required repetitions (either all or 1) - # or otherwise create a default zero return value for the batch signal as this it is - # required autographing this function by tensorflow - if self.adc_events is None: - for i in tf.range(n_samples): - batch_signal = batch_signal.write(i, tf.zeros(n_repetitions, - dtype=tf.complex64)) - else: - current_adc_events = tf.gather(self.adc_events, repetitions, axis=0) - current_adc_phase = tf.gather(self.adc_phase, repetitions, axis=0) + rf_waveforms, grad_waveforms, adc_on, adc_phase, adc_widx, batch_signal = \ + self._init_tensors(n_samples, n_repetitions, repetition_index) + + r_prev = tf.reshape(initial_position, [1, -1, 3]) + r_next = r_prev # Iterate over all time intervals to obtain the magnetization evolution - for t_idx in tf.range(0, tf.shape(self.dt_grid)[0]): + for t_idx in tf.range(1, tf.shape(self.dt_grid)[0]+1): # The shape of the magnetization tensor must always be (#reps, #batch, 3), with # needs to be specified for tf.autograph as: tf.autograph.experimental.set_loop_options( shape_invariants=[(magnetization, tf.TensorShape([None, None, 3])), - (r_next, tf.TensorShape([None, None, 3])), - (gdotr_r, tf.TensorShape([None, None])), - (gdotr_l, tf.TensorShape([None, None])) - ]) + (r_next, tf.TensorShape([None, None, 3]))]) # Increment particle positions and if available field definitions at given locations # r_next.shape = (#particles, #nsteps=1) - r_next, additional_fields = trajectory_module.increment_particles( - r_next[0], self.dt_grid[t_idx]) + r_next, additional_fields = trajectory_module.increment_particles(r_prev[0], self.dt_grid[t_idx-1]) r_next = r_next[tf.newaxis] # add dummy #reps-dim for subsequent broadcasting # Evaluate effective flip angle for given interval and apply the hard-pulse rotation if self.rf_waveforms is not None: - rf = current_rf_waveforms[:, t_idx:t_idx + 2] - alpha = (tf.complex(self.dt_grid[t_idx] * self.gamma, 0.) + rf = rf_waveforms[:, t_idx-1:t_idx+1] + alpha = (tf.complex(self.dt_grid[t_idx-1] * self.gamma, 0.) * (rf[:, 1] + rf[:, 0]) / 2) magnetization = self.hard_pulse(alpha[:, tf.newaxis], magnetization) # Evaluate effective phase increment for given interval including the submodule # effects and apply the corresponding rotation of the magnetization vectors + delta_phi = tf.zeros_like(M0, dtype=tf.float32) if self.gradient_waveforms is not None: - gdotr_r = tf.reduce_sum(current_gradient_waveforms[..., t_idx + 1, :] * r_next, - axis=-1) - delta_phi = (gdotr_l + gdotr_r) / 2. * self.gamma * self.dt_grid[t_idx] - gdotr_l = gdotr_r - for submod in self._submodules: - phi_s = submod( - gradient_wave_form=current_gradient_waveforms[:, 0, t_idx:t_idx + 1, :], - trajectories=r_next[:, :, tf.newaxis], dt=self.dt_grid[t_idx], - **additional_fields, **kwargs) - if tf.shape(phi_s)[0] < tf.shape(delta_phi)[0]: - phi_s = tf.tile(phi_s, [tf.shape(delta_phi)[0], 1, 1]) - # phi_s shape = (#reps, #particles, #nsteps=1) - tf.ensure_shape(phi_s, (None, None, 1)) - delta_phi = tf.reshape(delta_phi + phi_s[:, :, 0], [n_repetitions, -1]) - - # delta phi for current interval with (nreps, nbatch) dimension - tf.ensure_shape(delta_phi, (None, None)) - magnetization = self.phase_increment(delta_phi[:, tf.newaxis], magnetization) + grad = grad_waveforms[..., t_idx-1:t_idx+1, :] + gdotr_l = tf.reduce_sum(grad[..., 0, :] * r_prev, axis=-1) + gdotr_r = tf.reduce_sum(grad[..., 1, :] * r_next, axis=-1) + delta_phi += (gdotr_l + gdotr_r) / 2. * self.gamma * self.dt_grid[t_idx-1] + + for submod in self._submodules: + phi_s = submod( + gradient_wave_form=grad_waveforms[:, 0, t_idx:t_idx+1, :], + trajectories=r_next[:, :, tf.newaxis], dt=self.dt_grid[t_idx-1], + **additional_fields, **kwargs) + if tf.shape(phi_s)[0] < tf.shape(delta_phi)[0]: + phi_s = tf.tile(phi_s, [tf.shape(delta_phi)[0], 1, 1]) + # phi_s shape = (#reps, #particles, #nsteps=1) + tf.ensure_shape(phi_s, (None, None, 1)) + delta_phi = tf.reshape(delta_phi + phi_s[:, :, 0], [n_repetitions, -1]) + + # delta phi for current interval with (nreps, nbatch) dimension + tf.ensure_shape(delta_phi, (None, None)) + magnetization = self.phase_increment(delta_phi[:, tf.newaxis], magnetization) # Relax magnetization for the time interval - magnetization = self.relax(T1, T2, self.dt_grid[t_idx], magnetization) + magnetization = self.relax(T1, T2, self.dt_grid[t_idx-1], magnetization) # If adc-events are specified and the current interval is set to 'on' calculate # the weighed sum over all particles and demodulate according to the current ADC # phase. The result is then written into the signal return array. - # Note: To make this work with tf.function, the conditional needs an 'else' that - # also writes into the batch_signal return array at the last position. Therefore, - # the last entry is cut-off in the end, which implicitly assumes the last interval - # is never an active ADC-sample if self.adc_events is not None: - n_active_adc_channels = tf.reduce_sum(current_adc_events[:, t_idx + 1], axis=0) - if n_active_adc_channels >= 1.: - signal = tf.reduce_sum(magnetization[..., 0] * tf.complex(M0, 0.), axis=1) - signal = signal * tf.complex(current_adc_events[:, t_idx + 1], 0.) - demodulated_signal = signal * tf.exp( - tf.complex(0., -current_adc_phase[:, t_idx + 1])) - batch_signal = batch_signal.write(result_writing_idx, demodulated_signal) - result_writing_idx += 1 - else: - batch_signal = batch_signal.write(result_writing_idx, - tf.zeros(n_repetitions, - dtype=tf.complex64)) - batch_signal = batch_signal.stack() - return magnetization, r_next, tf.transpose(batch_signal[:-1], [1, 0]) + batch_signal = self._sample(adc_on[:, t_idx], adc_phase[:, t_idx], + adc_widx[:, t_idx], magnetization, M0, + r_next, batch_signal) + + r_prev = r_next + + transpose_indices = tf.range(tf.shape(tf.shape(batch_signal))[0]) + transpose_indices = tf.tensor_scatter_nd_update(transpose_indices, [[0], [1]], [1, 0]) + return magnetization, r_next, tf.transpose(batch_signal, transpose_indices) + + def _init_tensors(self, n_samples: int, n_repetitions: int, repetition_index: int): + """ Initializes / slices the waveforms according to the specified current repetitions and + number of samples. + + :param n_samples: + :param n_repetitions: + :param repetition_index: + :return: - rf_waveform (#reps, #steps) or None + - grad_waveform (#reps, #steps, 3) or None + - adc_events (#reps, #steps) or None + - adc_phase (#reps, #steps) or None + - adc_writing_idx (#reps, #steps) or None + - batch_signal (#reps, #samples) -> used for signal accumulation + """ + batch_signal = tf.zeros(shape=(n_samples, n_repetitions), dtype=tf.complex64) + repetitions = tf.range(repetition_index, repetition_index + n_repetitions, 1) + + # If rf waveforms are specified, gather all required repetitions (either all or 1) + if self.rf_waveforms is not None: + current_rf_waveforms = tf.gather(self.rf_waveforms, repetitions, axis=0) + else: + current_rf_waveforms = None + + # If gradient waveforms are specified, gather all required repetitions (either all or 1) + # and calculate the left value of for trapezoidal phase integration + if self.gradient_waveforms is not None: + current_gradient_waveforms = tf.gather(self.gradient_waveforms, repetitions, + axis=0)[:, tf.newaxis] + else: + current_gradient_waveforms = None + + # If adc events are specified, gather all required repetitions (either all or 1) + # or otherwise create a default zero return value for the batch signal as this it is + # required autographing this function by tensorflow + if self.adc_events is not None: + current_adc_events = tf.gather(self.adc_events, repetitions, axis=0) + current_adc_phase = tf.gather(self.adc_phase, repetitions, axis=0) + adc_writing_indices = tf.cast(tf.math.cumsum(current_adc_events, axis=1) - 1, tf.int32) + else: + current_adc_events, current_adc_phase, adc_writing_indices = None, None, None + + return (current_rf_waveforms, current_gradient_waveforms, current_adc_events, + current_adc_phase, adc_writing_indices, batch_signal) + + @staticmethod + def _sample(adc_events: tf.Tensor, adc_phase: tf.Tensor, adc_writing_idx: tf.Tensor, + magnetization: tf.Tensor, M0: tf.Tensor, r: tf.Tensor, batch_signal: tf.Tensor): + """Sums up the :math:`m_{xy}^{+}` magnetization weighted with the associated proton-density + of the particles and adds the signal to the batch-signal accumulator ad index specified + in adc_writing_index. + + :param adc_events: (#reps, ) + :param adc_phase: (#reps, ) + :param adc_writing_idx: (#reps, ) + :param magnetization: (#batch, ) + :param M0: (n_repetitions/1, #batch) + :param r: (n_repetitions/1, #batch, 3) + :param batch_signal: (#samples, #reps) + :return: + """ + n_active_adc_channels = tf.reduce_sum(adc_events, axis=0) + if n_active_adc_channels >= 1.: + signal = tf.reduce_sum(magnetization[..., 0] * tf.complex(M0, 0.), axis=1) + signal = signal * tf.complex(adc_events, 0.) + demodulated_signal = signal * tf.exp(tf.complex(0., -adc_phase)) + n_reps = tf.shape(adc_writing_idx)[0] + indices = tf.stack([adc_writing_idx, tf.range(n_reps, dtype=tf.int32)], axis=1) + in_bounds = tf.where(adc_writing_idx > -1) + indices = tf.gather(indices, in_bounds)[0] + demodulated_signal = tf.gather(demodulated_signal, in_bounds)[0] + batch_signal = tf.tensor_scatter_nd_add(batch_signal, indices, demodulated_signal) + return batch_signal def reset(self): """ Resets the time-signal accumulators, i.e. sets all entries to 0""" diff --git a/cmrsim/bloch/_multi_coil.py b/cmrsim/bloch/_multi_coil.py new file mode 100644 index 0000000000000000000000000000000000000000..8157cbe288c758e890fb2cccc74c6489f6b36398 --- /dev/null +++ b/cmrsim/bloch/_multi_coil.py @@ -0,0 +1,121 @@ +"""Contains implementation of Bloch-solving operator including multiple receive coils""" + +__all__ = ["ParallelReceiveBlochOperator"] + +from typing import List, Tuple, Union, TYPE_CHECKING + +import tensorflow as tf +import numpy as np + +from cmrsim.bloch._generic import GeneralBlochOperator +from cmrsim.bloch.submodules import BaseSubmodule +from cmrsim.trajectory import StaticParticlesTrajectory + +if TYPE_CHECKING: + import cmrsim.analytic.contrast.coil_sensitivities + + +class ParallelReceiveBlochOperator(GeneralBlochOperator): + """Extends the GeneralBlochOperator with parallel receive functionality. The only additionally + required input is an instance of a CoilSensitivity contrast module from the + cmrsim.analytic.contrast submodule, which is used internally to apply a spatially dependent + weight to the transverse magnetization before summing, in the process of acquiring signals. + + The resulting signal is also stored in the 'time_signal_acc' attribute, where for the + parallel receive case each entry has an additional axis corresponding to the number of coils. + + :param name: Verbose name of the module (non-functionally relevant) + :param gamma: gyromagnetic ratio in rad/ms/m + :param coil_module: + :param time_grid: (n_steps) - tf.float32 - time definition for numerical integration steps + :param gradient_waveforms: (#repetitions, #steps, 3) - tf.float32 - Definition of gradients + in mT/m + :param rf_waveforms: (#repetitions, #steps) - tf.complex64 - + :param adc_events: (#repetition, #steps, 2) - tf.float32 - last axis contains the adc-on + definitions and the adc-phase + :param submodules: List[tf.Modules] implementing the __call__ function returning the + resulting phase increment contribution + :param device: str + """ + #: Instance of the cmrsim.analytic.contrast.CoilSensitivity class that is used to compute + #: the positionally dependent sensitivity factors prior to sampling + coil_lookup_module: 'cmrsim.analytic.contrast.coil_sensitivities.CoilSensitivity' + + #: + n_coils: int + + def __init__(self, name: str, + gamma: float, + coil_module: Union['cmrsim.analytic.contrast.coil_sensitivities.CoilSensitivity'], + time_grid: tf.Tensor, + adc_events: tf.Tensor, + gradient_waveforms: tf.Tensor = None, + rf_waveforms: tf.Tensor = None, + submodules: Tuple[BaseSubmodule, ...] = (), + device: str = None): + """ + + + """ + super().__init__(name=name, gamma=gamma, time_grid=time_grid, + gradient_waveforms=gradient_waveforms, rf_waveforms=rf_waveforms, + adc_events=adc_events, submodules=submodules, device=device) + self.coil_lookup_module = coil_module + self.n_coils = self.coil_lookup_module.expansion_factor + + # Expand the dimension of the signal accumulator to incorporate coil-maps + with tf.device("CPU"): + self.time_signal_acc = [ + tf.Variable(initial_value=tf.zeros([v.read_value().shape[0], self.n_coils], dtype=tf.complex64), + shape=(None, self.n_coils), dtype=tf.complex64, trainable=False, + name=f"k_segment_{idx}") for idx, v in enumerate(self.time_signal_acc)] + + def _init_tensors(self, n_samples: int, n_repetitions: int, repetition_index: int): + """ Initializes / slices the waveforms according to the specified current repetitions and + number of samples. + + :param n_samples: + :param n_repetitions: + :param repetition_index: + :return: - rf_waveform (#reps, #steps) or None + - grad_waveform (#reps, #steps, 3) or None + - adc_events (#reps, #steps) or None + - adc_phase (#reps, #steps) or None + - adc_writing_idx (#reps, #steps) or None + - batch_signal (#reps, #samples) -> used for signal accumulation + """ + (rf_waveforms, gradient_waveforms, adc_events, adc_phase, adc_writing_indices, + batch_signal) = super()._init_tensors(n_samples, n_repetitions, repetition_index) + batch_signal = tf.zeros(shape=(n_samples, n_repetitions, self.n_coils), dtype=tf.complex64) + return (rf_waveforms, gradient_waveforms, adc_events, adc_phase, + adc_writing_indices, batch_signal) + + def _sample(self, adc_events: tf.Tensor, adc_phase: tf.Tensor, adc_writing_idx: tf.Tensor, + magnetization: tf.Tensor, M0: tf.Tensor, r: tf.Tensor, batch_signal: tf.Tensor): + """Sums up the :math:`m_{xy}^{+}` magnetization weighted with the associated proton-density + of the particles and adds the signal to the batch-signal accumulator ad index specified + in adc_writing_index. Before summation, a coil-sensitivity weighting is applied + + :param adc_events: (#reps, ) + :param adc_phase: (#reps, ) + :param adc_writing_idx: (#reps, ) + :param magnetization: (#batch, ) + :param M0: (n_repetitions/1, #batch) + :param r: (n_repetitions/1, #batch, 3) particle positions + :param batch_signal: (#samples, #reps, #coils) + :return: + """ + n_active_adc_channels = tf.reduce_sum(adc_events, axis=0) + if n_active_adc_channels >= 1.: + sense_factors = self.coil_lookup_module.lookup_complex_coil_factors(r) + weighted_magnetization = magnetization[..., 0, tf.newaxis] * sense_factors + signal = tf.reduce_sum(weighted_magnetization * tf.complex(M0[..., tf.newaxis], 0.), axis=1) + signal = signal * tf.complex(adc_events[..., tf.newaxis], 0.) + demodulated_signal = signal * tf.exp(tf.complex(0., -adc_phase[..., tf.newaxis])) + n_reps = tf.shape(adc_writing_idx)[0] + indices = tf.stack([adc_writing_idx, tf.range(n_reps, dtype=tf.int32)], axis=1) + in_bounds = tf.where(adc_writing_idx > -1) + indices = tf.gather(indices, in_bounds)[0] + demodulated_signal = tf.gather(demodulated_signal, in_bounds)[0] + batch_signal = tf.tensor_scatter_nd_add(batch_signal, indices, demodulated_signal) + return batch_signal diff --git a/cmrsim/bloch/submodules.py b/cmrsim/bloch/submodules.py index 9df756633425ba0145e2516f602376bb3637b47c..bc656184b66416f9436c4bddd4b14906e3d53fce 100644 --- a/cmrsim/bloch/submodules.py +++ b/cmrsim/bloch/submodules.py @@ -122,14 +122,14 @@ class ConcomitantFields(BaseSubmodule): """ with tf.device(self.device): tf.Assert(tf.shape(tf.shape(gradient_wave_form))[0] == 3, - ["Shape missmatch in ConcomittantFields for specified gradient", + ["Shape missmatch in ConcomitantFields for specified gradient", tf.shape(gradient_wave_form)]) gradient_wave_form = tf.expand_dims(gradient_wave_form, axis=1) if tf.shape(tf.shape(trajectories))[0] == 3: trajectories = tf.expand_dims(trajectories, axis=0) tf.Assert(tf.shape(tf.shape(gradient_wave_form))[0] == 4, - ["Shape missmatch in ConcomittantFields for trajectories", + ["Shape missmatch in ConcomitantFields for trajectories", tf.shape(trajectories)]) max_shape = tf.reduce_max(tf.stack([tf.shape(gradient_wave_form), diff --git a/cmrsim/datasets/__init__.py b/cmrsim/datasets/__init__.py index 549dc1007dc69cfb84b89a43fcab33c0d1247ccd..b6b1c70a49a585466a6d7150a2a3c4e8e1c69c51 100644 --- a/cmrsim/datasets/__init__.py +++ b/cmrsim/datasets/__init__.py @@ -10,4 +10,5 @@ from cmrsim.datasets._analytic import AnalyticDataset from cmrsim.datasets._flow import RefillingFlowDataset from cmrsim.datasets._bloch import BlochDataset from cmrsim.datasets._cardiac_mesh import MeshDataset, CardiacMeshDataset +from cmrsim.datasets._regular_grid import RegularGridDataset diff --git a/cmrsim/datasets/_flow.py b/cmrsim/datasets/_flow.py index faf8e9e46c75a4cdf315ff1adde321f0828ac8db..9526c19b53c77460db43dfaaa3a62709206fb1d0 100644 --- a/cmrsim/datasets/_flow.py +++ b/cmrsim/datasets/_flow.py @@ -12,6 +12,7 @@ import pyvista from scipy.ndimage import gaussian_filter + # pylint: disable=too-many-instance-attributes class RefillingFlowDataset(tf.Module): """ Wraps a vtk-mesh of containing a velocity field. Functionality tackles the problem of in-/ diff --git a/cmrsim/datasets/_regular_grid.py b/cmrsim/datasets/_regular_grid.py new file mode 100644 index 0000000000000000000000000000000000000000..5b9c185ce10bd3bdc346475ee532efd14d000b32 --- /dev/null +++ b/cmrsim/datasets/_regular_grid.py @@ -0,0 +1,299 @@ +__all__ = ["RegularGridDataset", ] + +from typing import Dict, Tuple, Union, Sequence + +import tensorflow as tf +import numpy as np +import pyvista +from pint import Quantity + + +class RegularGridDataset: + """Dataset class containing functionality, that requires a regular-grid representation of the + data. Multiple class-methods define inputs from which an instance of this class can be + created. + + :param mesh: Instance of a pyvista.UniformGrid which is copied on instantiation + """ + #: Instance of a pyvista.UniformGrid used as basis for all contained computations + mesh: pyvista.UniformGrid + + def __init__(self, mesh: pyvista.UniformGrid): + self.mesh = mesh.copy() + + def get_phantom_def(self, filter_by: str = None, keys: Sequence[str] = None, + perturb_positions_std: float = None) -> Dict[str, np.ndarray]: + """Extracts points and their corresponding physical properties from the contained mesh + and returns them as dictionary of flattened array + + :param filter_by: used to determine which points to extract from the mesh by using the + indices = np.where(self.mesh[filter_by]) + :param keys: sequence of strings defining the field arrays returned as dictionary keys + :param perturb_positions_std: if specified the positions of the extracted mesh-points are + perturbed by sampling a normal distribution with zero mean and given standard + deviation to reduce descretization artifacts in simulation + :return: + """ + + if filter_by is not None: + indices = np.where(self.mesh[filter_by]) + else: + indices = np.where(np.ones_like(self.mesh.points[:, 0])) + + if keys is None: + keys = self.mesh.array_names + + r_vectors = np.array(self.mesh.points[indices], dtype=np.float32).reshape(-1, 3) + if perturb_positions_std is not None: + r_vectors += np.random.normal(0, perturb_positions_std, size=r_vectors.shape) + + phantom = dict(r_vectors=r_vectors) + phantom.update({k: np.squeeze(self.mesh[k][indices].reshape(r_vectors.shape[0], -1)) + for k in keys}) + return phantom + + def add_field(self, key: str, data: np.ndarray): + """Adds field to self.mesh from a numpy array (takes care of correct ordering from + default numpy C-order to F-order used in pyvista). + + :param key: name of the new field + :param data: (x, y, z, N) dimensional array containing arbitrary data + """ + n_channels = np.prod(data.shape) // np.prod(self.mesh.dimensions) + flat_data = data.reshape([-1, n_channels], order="F") + self.mesh[key] = flat_data + + def resample_spacing(self, spacing: Quantity): + """Resamples the contained self.mesh to a new spacing, which might be necessary to reduce + discretization artifacts on simulation. + + :param spacing: New target spacing (3, ) + """ + r_max = np.max(self.mesh.points, axis=0) + r_min = np.min(self.mesh.points, axis=0) + nbins = (Quantity(r_max - r_min, "m") / spacing.to("m")).m.astype(int) + new_mesh = pyvista.UniformGrid(dimensions=nbins, spacing=spacing.m_as("m"), origin=r_min) + self.mesh = new_mesh.sample(self.mesh) + + def select_slice(self, slice_normal: np.array, spacing: Quantity, + slice_position: Quantity, field_of_view: Quantity, + readout_direction: np.array = None, + in_mps: bool = False) -> pyvista.UniformGrid: + """Creates a uniform grid (slice) according to the specified parameters, and uses it to + probe the self.mesh (see example). Slice (coordinates) can either be returned in mps + or in global coordinates. + + .. Dropdown:: Example Visualization + + .. image:: ../../_static/api/dataset_regular_grid_select_slice.png + + :param slice_normal: (3, ) normal vector of the extracted slice + :param spacing: (3, ) dimension of voxels in all directions + :param slice_position: (3, ) vector pointing to the slice-center + :param field_of_view: (3, ) total extent of the extracted slice in 3 dimensions + :param readout_direction: (3, ) if not provided, is assumed to be x (part orthogonal to + specified slice normal). Used to rotate the extracted slice about the slice-normal + :param in_mps: if true the slice is returned in MPS coordinates + :return: + """ + from scipy.spatial.transform import Rotation + slice_normal /= np.linalg.norm(slice_normal) + + if readout_direction is None: + if np.allclose(slice_normal, np.array([1., 0., 0.]), atol=1e-3): + readout_direction = np.array([0., 1., 0.]) + else: + readout_direction = np.array([1., 0., 0.]) + else: + + if np.isclose( + np.dot(slice_normal, readout_direction / np.linalg.norm(readout_direction)), + 1., atol=1e-3): + raise ValueError("Slice normal and readout-direction are parallel") + + readout_direction = readout_direction - np.dot(slice_normal, + readout_direction) * slice_normal + readout_direction /= np.linalg.norm(readout_direction) + pe_direction = np.cross(readout_direction, slice_normal) + pe_direction /= np.linalg.norm(pe_direction) + original_basis = np.eye(3, 3) + new_basis = np.stack([readout_direction, pe_direction, slice_normal]) + + rot, _ = Rotation.align_vectors(original_basis, new_basis) + + total_transform = np.eye(4, 4) + total_transform[:3, :3] = rot.as_matrix() + total_transform[:3, 3] = slice_position.m_as("m") + nbins = (field_of_view / spacing).m_as("dimensionless").astype(int) + new_slice = pyvista.UniformGrid(dimensions=nbins, spacing=spacing.m_as("m"), + origin=-field_of_view.m_as("m") / 2) + new_slice = new_slice.transform(total_transform, inplace=False) + new_slice = self.mesh.probe(new_slice) + + if in_mps: + inv_transform = np.eye(4, 4) + inv_transform[:3, :3] = rot.as_matrix().T + inv_transform[:3, 3] = -slice_position.m_as("m") + new_slice = new_slice.transform(inv_transform, inplace=False, + transform_all_input_vectors=True) + return new_slice + + def compute_offresonance(self, b0: Quantity, susceptibility_key: str) -> Quantity: + """Approximates the magnetic field variations in z-direction due to interfaces with varying + susceptibilities by applying a convolution (Fourier Product) with a dipole kernel to the + susceptibility field contained in self.mesh. The name of the mesh array is specified by + argument. + + Assumes that the main magnetic field is pointing in z-direction of the contained mesh. + The susceptibility os assumed to be specified as parts per million + (e.g. :math:`\chi_{air} = 0.36ppm`). + + The result is store in self.mesh["offres"] and is also returned as Quantity. + + .. Dropdown:: Literature Reference + + This function is based on the work of Frank-Zijlstra, published on Mathworks: + + Frank Zijlstra (2022). MRI simulation using FORECAST: Fourier-based Off-REsonanCe + Artifact simulation in the STeady-state + (https://www.mathworks.com/matlabcentral/fileexchange/56680-mri-simulation-using- + forecast-fourier-based-off-resonance-artifact-simulation-in-the-steady-state), + + MATLAB Central File Exchange. Retrieved December 8, 2022. + + :param b0: Field strength of main magnetic field pointing in z-direction in Tesla + :param susceptibility_key: name of the array contained in self.mesh, used as susceptibility + :return: (X, Y, Z) array wrapped as Quantity with dimension ("T") + """ + chi_3d = self.mesh[susceptibility_key].reshape(self.mesh.dimensions, order="F") + # Pad the chi-map to a multiple of four: + residual_per_dim = np.mod(self.mesh.dimensions, 4) + chi_3d = np.pad(chi_3d, + np.stack([np.zeros_like(residual_per_dim), residual_per_dim], axis=1), + mode="constant", constant_values=0.) + voxel_size = Quantity(self.mesh.spacing, "m") + field_of_view = voxel_size * Quantity(chi_3d.shape, "dimensionless") + + # The fourier based convolution (from original implementation): + # 1) The forward Fourier Transform (FT): Instead of doing an extra convolution to calculate + # the effect of aliasing (as proposed in the article), here we use the 'imaginary channel' + # to do this. This is done by downscaling the artificial environment, and adding this as + # an imaginary component to the original (unpadded) susceptibility distribution. + # Together this is referred to as the "DUAL" distribution: + # down_sampled_chi = tf.nn.avg_pool3d(chi_3d[np.newaxis, ..., np.newaxis], ksize=(2, 2, 2), + # strides=(2, 2, 2), padding="VALID") + # imaginary_channel = self._zc(down_sampled_chi[0, ..., 0]) + # ft_chi_3d_dual = tf.signal.fft3d(chi_3d + 1j * imaginary_channel) + ft_chi_3d_dual = tf.signal.fft3d(chi_3d) + # 2) multiplication with the dipole function + dipole_kernel = self._compute_fourier_dipole_kernel(field_of_view, + np.array(ft_chi_3d_dual.shape)) + ft_field_3d = ft_chi_3d_dual * dipole_kernel + # 3) inverse FT to spatial domain + field_3d_dual = tf.signal.ifft3d(ft_field_3d) + # 4) Subtract the up-scaled cropped center of the aliasing + # field_3d = tf.math.real(field_3d_dual) - UC(tf.math.imag(field_3d_dual)) + field_3d = tf.math.real(field_3d_dual) + i, j, k = [_ if _ != 0 else None for _ in -residual_per_dim] + field_3d = field_3d.numpy()[:i, :j, :k] + self.mesh["offres"] = b0.to("mT") * field_3d.reshape(-1, order="f") + return field_3d * b0.to("mT") + + @staticmethod + def _compute_fourier_dipole_kernel(fov: Quantity, n_samples: np.ndarray): + """ Computes the dipole kernel for a given FOV and 3D array + :param fov: (3, ) length - quantity + :param n_samples: (3, ) integer number of regular points per dimension + :return: + """ + k_squared = [np.arange(-n/2, n/2, 1) / fov[i].m_as("m") for i, n in enumerate(n_samples)] + k_squared = [np.fft.ifftshift(k)**2 for k in k_squared] + [kx2, ky2, kz2] = np.meshgrid(*k_squared, indexing="ij") + kernel = 1/3 - kz2 / (kx2 + ky2 + kz2) + kernel[0, 0, 0] = 0 + return kernel + + @classmethod + def from_unstructured_grid(cls, input_mesh: pyvista.DataSet, pixel_spacing: Quantity, + padding: Quantity = None) -> 'RegularGridDataset': + """Resamples a given unstructured mesh onto a regular grid and instantiates a + RegularGridDataset from it. + + :param input_mesh: pyvista.Dataset + :param pixel_spacing: Quantity (3,) specifying the 3D spacing of the uniform grid + :param padding: Space added symmetrically around the bounding box of the input_mesh, + defaults to (0, 0, 0) "meter" + """ + if padding is None: + padding = Quantity([0., 0., 0.], "m") + r_max = np.max(input_mesh.points, axis=0) + padding.m_as("m") / 2 + r_min = np.min(input_mesh.points, axis=0) - padding.m_as("m") / 2 + nbins = (Quantity(r_max - r_min, "m") / pixel_spacing.to("m")).m.astype(int) + mesh = pyvista.UniformGrid(dimensions=nbins, spacing=pixel_spacing.m_as("m"), origin=r_min) + mesh = input_mesh.probe(mesh) + return cls(mesh) + + @classmethod + def from_numpy(cls, data: Union[np.ndarray, Dict[str, np.ndarray]], + origin_offset: Quantity = Quantity([0, 0, 0], "m"), extent: Quantity = None, + pixel_spacing: Quantity = None) -> 'RegularGridDataset': + """Constructs a uniform regular 3D grid of the same shape as the data provided as argument. + + Assumes that the origin is at the center of the provided array. + + :raises: ValueError: + - if neither extent nor pixel_spacing is specified + - data-argument is a dictionary containing arrays with non-identical shapes + + :param data: (X, Y, Z, ...) + :param origin_offset: (3, ) If not provided, the grid is positioned symmetrically around the + point (0., 0., 0.). Otherwise this is vector is added onto the + center-point. + :param extent: (3, ) Length-Quantity defining the extend of the provided data in 3 dimen + :param pixel_spacing: (3, ) + """ + + if (all(v is None for v in [pixel_spacing, extent]) or + all(v is not None for v in [pixel_spacing, extent])): + raise ValueError("Either 'extent'/'pixel-spacing' keyword argument must be specified") + + if isinstance(data, dict): + shapes = [v.shape for v in data.values()] + if not all(shapes[0][0:3] == x[0:3] for x in shapes): + raise ValueError("Data arrays specified in argument 'data' are not identical:\n" + + "\n".join([str(s) for s in shapes])) + else: + shapes = [data.shape, ] + + if extent is not None: + pixel_spacing = extent.to("m") / np.array(shapes[0]) + + shape = np.array(shapes[0][0:3]) + origin = - pixel_spacing * shape / 2 + origin_offset + + mesh = pyvista.UniformGrid(dimensions=shape, spacing=pixel_spacing.m_as("m"), + origin=origin.m_as("m")) + + for name, arr in data.items(): + mesh[name] = arr.reshape([shape.prod(), *arr.shape[3:]], order="F") + + return cls(mesh) + + @classmethod + def from_shepp_logan(cls, map_size: Tuple[int, int, int], extent: Quantity) \ + -> 'RegularGridDataset': + """Creates a 3D Shepp-Logan phantom and stores it in a RegularGridDataset. + + :param map_size: Target number of voxels per dimension + :param extent: Target spatial extent of the resulting shepp-logan phantom. + """ + try: + import phantominator + except ImportError: + raise ImportError("The python-package phantominator is necessary to create a shepp-" + "logan phantom. To install it, type:\n pip install phantominator") + + M0, T1, T2 = phantominator.shepp_logan(map_size, MR=True) + M0, T1, T2 = [t.transpose(2, 1, 0) for t in [M0, T1, T2]] + data = dict(rho=M0, T1=T1 * 1000, T2=T2*1000) + return cls.from_numpy(data=data, extent=extent) diff --git a/cmrsim/trajectory/_flow.py b/cmrsim/trajectory/_flow.py index c5449b9804cd6e821e3f2924e79509112cb0db43..a1d73649791f768db250a3441b3ae52acb8fa613 100644 --- a/cmrsim/trajectory/_flow.py +++ b/cmrsim/trajectory/_flow.py @@ -49,11 +49,11 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule): :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)] Physical extend of the gridded velocity fields. :param additional_dimension_mapping: List of name/len pairs to unstack the looked up - values. If not None this mapping must explain all - 'n' additional dimension in the parameter - velocity-field e.g. ('off_res', 1) meaning that - the 1 index (starting from 3) in velocity field - describes the off-resonance at all time points. + values. If not None this mapping must explain all + 'n' additional dimension in the parameter + velocity-field e.g. [('off_res', 1)] meaning that + off-resonance for all time points has a size of 1, + starting from index 3 in velocity field. :param kwargs: - device: str defaults to CPU:0 """ @@ -70,9 +70,9 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule): :param additional_dimension_mapping: List of name/len pairs to unstack the looked up values. If not None this mapping must explain all 'n' additional dimension in the parameter - velocity-field e.g. ('off_res', 1) meaning that - the 1 index (starting from 3) in velocity field - describes the off-resonance at all time points. + velocity-field e.g. [('off_res', 1)] meaning that + off-resonance for all time points has a size of 1, + starting from index 3 in velocity field. :param kwargs: - device: str defaults to CPU:0 """ @@ -88,6 +88,8 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule): field_dims.insert(0, 0) self._field_names = {k: (field_dims[i], field_dims[i + 1]) for i, k in enumerate(field_labels)} + else: + self._field_names = {} super().__init__(velocity_field, map_dimensions, name="flow_trajectory_evaluation", **kwargs) @@ -175,19 +177,75 @@ class FlowTrajectory(BaseTrajectoryModule, LookUpTableModule): # pylint: disable=abstract-method class TurbulentTrajectory(FlowTrajectory): - """ Implements the path integration for particles in a temporally static velocity field + """ Numerically integrates the path of particles in a mean velocity field and incorporates + random turbulent velocity + fluctuations based on Langevin updates from sampled RST values. + Looks up and returns additional fields for the previous particle positions. + .. dropdown:: Example usage + .. code-block:: python + :caption: Instantiation - """ - def __init__(self, velocity_field: np.ndarray, map_dimensions: np.ndarray, **kwargs): - """ + velocity_field_3d = ... # shape - (X, Y, Z, 10) + map_dimensions = [(x_lo, x_hi), (y_lo, y_hi), (z_lo, z_hi)] + field_list = [] # ... no addition a field mapping + trajectory_module = TurbulentTrajectory(velocity_field_3d, map_dimensions, field_list) + + .. code-block:: python + :caption: Call function + + r_init = ... # initial_positions of shape (N, 3) + timing = np.arange(0., 1.01, 0.01) + r_at_timing, fields = trajectory_module(r_init, timing, dt_max=0.01, + return_velocities=True) + + .. code-block:: python + :caption: increment particle positions - :param flow_field: (X, Y, Z, 10) 3D map storing gridded velocities in m/s at positions - X, Y, Z, Chol11, Chol12, Chol13, Chol22, Chol23, Chol33, Tau. + r_init = ... # initial_positions of shape (N, 3) + delta_t = tf.constant(0.01, tf.float32) + r_new, fields = self.trajectory_module.increment_particles(r_init, dt=delta_t) + + + :param velocity_field: (X, Y, Z, 10+n) 3D map storing gridded fields at positions X, Y, Z. + Fields 0-2 : Mean velocities (X, Y, Z). Unit should be m/ms + Fields 3-8 : Lower triangular Cholesky decomposition of the turbulent + Reynolds stress tensor, L. Ordering is + (L11, L21, L22, L31, L32, L33) + Field 9 : Langevin timescale tau, in Hz. + Additional dimensions past 9 are treated as additional + data fields and returned :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)] Physical extend of the gridded velocity fields. + :param additional_dimension_mapping: List of name/len pairs to unstack the looked up + values. If not None this mapping must explain all + 'n' additional dimension in the parameter + velocity-field e.g. [('off_res', 1)] meaning that + off-resonance for all time points has a size of 1, + starting from index 10 in velocity field. + :param kwargs: - device: str defaults to CPU:0 + """ + def __init__(self, velocity_field: np.ndarray, map_dimensions: np.ndarray, + additional_dimension_mapping: List[Tuple[str, int]] = None, **kwargs): + """ + :param velocity_field: (X, Y, Z, 10+n) 3D map storing gridded fields at positions X, Y, Z. + Fields 0-2 : Mean velocities (X, Y, Z). Unit should be m/ms + Fields 3-8 : Lower triangular Cholesky decomposition of the turbulent + Reynolds stress tensor, L. Ordering is + (L11, L21, L22, L31, L32, L33) + Field 9 : Langevin timescale tau, in Hz. + Additional dimensions past 9 are treated as additional + data fields and returned + :param map_dimensions: (3, 2) -> [(xlow, xhigh), (ylow, yhigh), (zlow, zhigh)] + Physical extend of the gridded velocity fields. + :param additional_dimension_mapping: List of name/len pairs to unstack the looked up + values. If not None this mapping must explain all + 'n' additional dimension in the parameter + velocity-field e.g. [('off_res', 1)] meaning that + off-resonance for all time points has a size of 1, + starting from index 10 in velocity field. :param kwargs: - device: str defaults to CPU:0 """ self.prev_v_turb = tf.Variable(tf.zeros([1, 3]), shape=(None, 3), dtype=tf.float32) @@ -206,11 +264,29 @@ class TurbulentTrajectory(FlowTrajectory): @tf.function(jit_compile=False, reduce_retracing=True) def increment_particles(self, particle_positions: tf.Tensor, dt: tf.Tensor, return_velocities: bool = False): - data = self.linear_interpolation(r_vectors=particle_positions) - impuls = tf.random.normal([tf.shape(particle_positions)[0], 3]) + """ For given particle positions and specified time-interval returns the values of all + additional fields at the initial particle position and moves the particles according + to: + + .. math:: + + r_{t+\\delta t} = r_{t} + \\delta t (U(r_{t}) + u_{t}) + + u_{t} = u_{t-\\delta t} e^{-\\frac{\\delta t}{\\tau}} + \\zeta \\sqrt{1-e^{-2 \\frac{\\delta t}{\\tau}}} - # Langevin equation + where U is the mean velocity in m/ms, u is fluctuating velocity in m/ms, :math:`\\tau` is the Langevin timescale + in ms, and :math:`\\zeta` is a random Gaussian sample from the Reynolds stress tensor. + + :param particle_positions: (-1, 3) batch of 3D particle coordinates + :param dt: scalar value in ms + :param return_velocities: if True the additional field lookup contains the field "velocity" + :return: new particle positions (-1, 3), dict(\*\*aditional_fields[-1, c]) + """ + data = self.linear_interpolation(r_vectors=particle_positions) v_mean = data[:, 0:3] + + # Random sample from RST (Cholesky decomposition) + impuls = tf.random.normal([tf.shape(particle_positions)[0], 3]) vtu = impuls[:, 0] * data[:, 3] vtv = tf.reduce_sum(impuls[:, :2] * data[:, 4:6], axis=1) vtw = tf.reduce_sum(impuls * data[:, 6:9], axis=1) @@ -227,12 +303,12 @@ class TurbulentTrajectory(FlowTrajectory): # New turbulence velocity v_turb = damping_update[:, tf.newaxis] * tf.stack([vtu, vtv, vtw], axis=1) - refmap = tf.reshape(tf.floor(data[:, 10]), [tf.shape(v_turb)[0], 1]) + #refmap = tf.reshape(tf.floor(data[:, 10]), [tf.shape(v_turb)[0], 1]) # Old + new turbulence velocity v_turb = self.prev_v_turb * tf.reshape(damping, [tf.shape(v_turb)[0], 1]) + v_turb - v_turb = refmap * v_turb + (refmap - 1) * self.prev_v_turb + #v_turb = refmap * v_turb + (refmap - 1) * self.prev_v_turb v_total = v_mean + v_turb additional_fields = self._unstack_add_dims(data[:, 10:]) @@ -243,8 +319,8 @@ class TurbulentTrajectory(FlowTrajectory): return particle_positions + v_total * dt, additional_fields def warmup_step(self, particle_positions: tf.Tensor): - """ Updates turbulence velocity of all particles based on sampled RST values at their - positions. + """ Updates turbulence velocity of all particles based on sampled RST values at their positions. + Used to initialize turbulence velocity to prevent transient states due to long Langevin timescales :param particle_positions: (-1, 3) batch of 3D particle coordinates """ diff --git a/cmrsim/utils/__init__.py b/cmrsim/utils/__init__.py index c517baa1a665162e66f4ba1507cadaa570b2aff6..ba1478cabec16d7dade84f31d4291b0ab380cd57 100644 --- a/cmrsim/utils/__init__.py +++ b/cmrsim/utils/__init__.py @@ -6,5 +6,4 @@ __all__ = ["coordinates", "snr", "particle_properties", "display"] import cmrsim.utils.coordinates import cmrsim.utils.snr import cmrsim.utils.particle_properties -import cmrsim.utils.offresonance import cmrsim.utils.display diff --git a/cmrsim/utils/offresonance.py b/cmrsim/utils/offresonance.py deleted file mode 100644 index 2266aa9898a22268da827818f82e2683784bd28c..0000000000000000000000000000000000000000 --- a/cmrsim/utils/offresonance.py +++ /dev/null @@ -1,44 +0,0 @@ -""" This modules contains utility to create offresonance related parameters""" -__all__ = ["generate_inhomogeneity_map"] - -import numpy as np -from pint import Quantity - - -def generate_inhomogeneity_map(labels, B0: Quantity, sus_maps: dict, padsize): - """ Calculates inhomogeneity maps on discrete 3D maps of suseptibility values - - :param labels: 3D volume oof labels - :param B0: Field strength in - :param sus_maps: dictionary that maps a suszeptibility for all values in labels - :param padsize: - :return: - """ - labels_pad = np.pad(labels, padsize, constant_values=0) - - phantom = np.ones_like(labels_pad).astype(np.float64) * -7e-7 - for key, val in sus_maps.items(): - phantom[np.where(labels_pad == key)] = val - # phantom = Quantity(phantom, "") - - mu0 = Quantity(4. * np.pi * 1e-7, "H/m") - - # Equation 5 - center = np.ceil((np.array(phantom.shape) + 1) / 2) - x, z, y = np.meshgrid(*[np.arange(phantom.shape[i]) - center[i] for i in range(3)], indexing="ij") - corterm = 3 * z ** 2 / (x ** 2 + y ** 2 + z ** 2) - 1 - - # Equation 4 - Mz = B0.to("T") * phantom / mu0 - Mk = np.fft.fftshift(np.fft.fftn(Mz.m)) - Bk = (-mu0 * Mk / 3. * corterm).m - Bk[int(center[0]), int(center[1]), int(center[2])] = -B0.m_as("T") * 2.9e-8 / 3 - - result = np.fft.ifftn(np.fft.ifftshift(Bk)) - inhomogeneity = np.real(result[padsize:-padsize, padsize:-padsize, padsize:-padsize]) - background_mean = np.mean(inhomogeneity[np.where(labels == 0)]) - inhomogeneity -= background_mean - - pmax = np.percentile(np.abs(inhomogeneity[np.where(labels != 0)]), 95) - inhomogeneity /= pmax - return inhomogeneity diff --git a/docker/Dockerfile_cpu b/docker/Dockerfile_cpu index 3093a3fcd8d8f540de1dd9a8e8b8cf029619282f..41a6579092dfffff4aff7d8200ba4106b320129a 100644 --- a/docker/Dockerfile_cpu +++ b/docker/Dockerfile_cpu @@ -2,7 +2,7 @@ # docker build -t registry.ethz.ch/jweine/cmrsim/cpu -f Dockerfile_cpu . # Use the tensorflow image as base -FROM tensorflow/tensorflow:2.9.1 +FROM tensorflow/tensorflow:2.11.0 # maintainter from IBT LABEL maintainers="Jonathan Weine (weine@biomed.ee.ethz.ch)" @@ -16,13 +16,10 @@ RUN pip3 install h5py scipy tqdm deprecated RUN pip3 install pyvista ffmpeg ipyvtklink panel # Install cmrseq for sequence definitions. Latest version whl needs to be provided in context -#COPY $cmrseq_file $cmrseq_file -#RUN pip install $cmrseq_file -# When cmrseq will be public, this is not necessary anymore +RUN pip install cmrseq # Install cmrsim -COPY cmrsim-latest-py3-none-any.whl cmrsim-latest-py3-none-any.whl -RUN pip install cmrsim-latest-py3-none-any.whl +RUN pip install cmrsim ## # Environment diff --git a/docker/Dockerfile_gpu b/docker/Dockerfile_gpu index 4ce0bc75c714c19278955025f19111541728ea2f..d14f4e2cf15f35e687b8e5c3fa809c26d0344715 100644 --- a/docker/Dockerfile_gpu +++ b/docker/Dockerfile_gpu @@ -1,5 +1,7 @@ +## docker build -t registry.ethz.ch/jweine/cmrsim/gpu:latest -f ./Dockerfile_gpu . + # take the GPU python3 tensorflow image as base -FROM tensorflow/tensorflow:2.9.1-gpu +FROM tensorflow/tensorflow:2.11.0-gpu # maintainter from IBT LABEL maintainers="Jonathan Weine (weine@biomed.ee.ethz.ch)" @@ -12,14 +14,11 @@ RUN apt-get install -y libx11-6 libgl1-mesa-dev libxrender-dev RUN pip3 install h5py scipy tqdm deprecated RUN pip3 install pyvista ffmpeg ipyvtklink panel -# Install cmrseq for sequence definitions. Latest version whl needs to be provided in context -#COPY $cmrseq_file $cmrseq_file -#RUN pip install $cmrseq_file -# When cmrseq will be public, this is not necessary anymore +# Install cmrseq for sequence definitions. +RUN pip install cmrseq # Install cmrsim -COPY cmrsim-latest-py3-none-any.whl cmrsim-latest-py3-none-any.whl -RUN pip install cmrsim-latest-py3-none-any.whl +RUN pip install cmrsim ## # Environment diff --git a/docker/Dockerfile_jupyter b/docker/Dockerfile_jupyter index 7c214d6fdef5ffad9ec3a0658222548d04685cc7..599357abcf7c9a494fb93698d2875f944ec0d1b5 100644 --- a/docker/Dockerfile_jupyter +++ b/docker/Dockerfile_jupyter @@ -1,3 +1,4 @@ +## docker build -t registry.ethz.ch/jweine/cmrsim/jupyter:latest -f ./Dockerfile_jupyter . FROM registry.ethz.ch/jweine/cmrsim/gpu:latest # maintainter from IBT @@ -7,7 +8,7 @@ LABEL Description="This Docker image is meant to start MR simulations with the C # Set up jupyter lab with interactive plots and variable inspector RUN apt-get -y update RUN pip3 install jupyterlab -RUN pip3 install ipympl +RUN pip3 install ipympl==0.8.8 matplotlib==3.5 RUN apt-get install -y curl RUN curl -sL https://deb.nodesource.com/setup_14.x | bash -E RUN apt-get install -y nodejs @@ -16,8 +17,8 @@ RUN jupyter labextension install jupyter-matplotlib RUN pip install lckr-jupyterlab-variableinspector # Install pyvista and dependencies -RUN apt install -y libgl1-mesa-glx libxrender1 -RUN pip install vtk pyvista +RUN apt install -y libgl1-mesa-glx libxrender1 xvfb +RUN pip install vtk pyvista phantominator RUN apt-get install -y git RUN apt-get install apt-transport-https ca-certificates -y @@ -26,5 +27,4 @@ RUN update-ca-certificates ## Environment ## ENV LD_LIBRARY_PATH /usr/local/cuda/extras/CUPTI/lib64:$LD_LIBRARY_PATH RUN echo $LD_LIBRARY_PATH -ENV PATH /opt/conda/bin:$PATH diff --git a/docker/Dockerfile_unittest b/docker/Dockerfile_unittest index 7082a90c75a07ac46829764d54a33a8d6b4dd6af..9e752278e3484c66287c74cffa01fbac793435fe 100644 --- a/docker/Dockerfile_unittest +++ b/docker/Dockerfile_unittest @@ -11,8 +11,5 @@ RUN pip install coverage anybadge phantominator genbadge RUN pip install pylint pylint-exit RUN apt install -y libgl1-mesa-glx xvfb -COPY cmrseq-0.15.dev4-py3-none-any.whl cmrseq-0.15.dev4-py3-none-any.whl -RUN pip install cmrseq-0.15.dev4-py3-none-any.whl - # To push image: # docker push registry.ethz.ch/jweine/cmrsim:unittest \ No newline at end of file diff --git a/docker/cmrseq-0.15.dev4-py3-none-any.whl b/docker/cmrseq-0.15.dev4-py3-none-any.whl deleted file mode 100644 index ec4cc3610f303d909bb39d9614bd4743af47ad2b..0000000000000000000000000000000000000000 Binary files a/docker/cmrseq-0.15.dev4-py3-none-any.whl and /dev/null differ diff --git a/docs/source/_static/api/dataset_regular_grid_select_slice.png b/docs/source/_static/api/dataset_regular_grid_select_slice.png new file mode 100644 index 0000000000000000000000000000000000000000..6325cc53a6b56d583bc41466348c97aea5e293b9 --- /dev/null +++ b/docs/source/_static/api/dataset_regular_grid_select_slice.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f4d9883a1c80a8c2e212ee817d7184f77be2d51c5ea130b5296d7c4fa89e5120 +size 85981 diff --git a/docs/source/examples.rst b/docs/source/examples.rst index 4b3c27c6440a3efab5aee2b8dbaa87b107ab2ee4..58aede87acb1ce5c4040c2df0e3fdec7210c907d 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -15,6 +15,11 @@ Datasets .. nbgallery:: example_gallery/datasets/cardiac_mesh_dataset + .. card:: Volumes on regular Grids + + .. nbgallery:: example_gallery/datasets/regular_grid_dataset + + Trajectory Modules ==================== @@ -68,7 +73,7 @@ Static Phantoms .. card:: Off-resonance with BSSFP - .. nbgallery:: example_gallery/bloch_simulation/static_BSSFP_submodules.ipynb + .. nbgallery:: example_gallery/bloch_simulation/static_BSSFP_offresonance Flowing particles ^^^^^^^^^^^^^^^^^^^ @@ -92,7 +97,7 @@ Flowing particles .. nbgallery:: example_gallery/bloch_simulation/flow_bssfp_submodules_aorta -Mesh Motions +Mesh Motion ^^^^^^^^^^^^^ .. card-carousel:: 2 @@ -120,6 +125,10 @@ Analytic Simulation .. nbgallery:: example_gallery/analytic_simulation/3_multiple_coils + .. card:: BSSFP on spherical off-resonance Phantom + + .. nbgallery:: example_gallery/analytic_simulation/BSSFP_offresonance + .. card:: BSSFP on texturized XCAT phantom .. nbgallery:: example_gallery/analytic_simulation/simulate_mrxcat_cmrsim \ No newline at end of file diff --git a/docs/source/format_datasets.py b/docs/source/format_datasets.py index 86b05ba2fef53cdac1d28299ed4a6d92be8158f6..8e9d3b5c64281a8883af8d06a9d79414347dcc33 100644 --- a/docs/source/format_datasets.py +++ b/docs/source/format_datasets.py @@ -2,6 +2,7 @@ import cmrsim import glob import os + def format_dataset_index(): all_modules = glob.glob(f"api/datasets/cmrsim.datasets.*.rst") all_modules = [s.replace("api/", "") for s in all_modules] diff --git a/notebooks/analytic_simulation/BSSFP_offresonance.ipynb b/notebooks/analytic_simulation/BSSFP_offresonance.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..cd7b41c50bdaef93b493912fa805de6ca23c3226 --- /dev/null +++ b/notebooks/analytic_simulation/BSSFP_offresonance.ipynb @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bb017045678b9a3a2299fa7512712d6d045d67a3edd5158adc6ca6d7cbb57964 +size 539217 diff --git a/notebooks/bloch_simulation/static_BSSFP_offresonance.ipynb b/notebooks/bloch_simulation/static_BSSFP_offresonance.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..94eb5d2c75ae82a0bc06064a4480c15c70ccf88b --- /dev/null +++ b/notebooks/bloch_simulation/static_BSSFP_offresonance.ipynb @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2c0a3d4bc78511d90712facf3dd3de544e92fe84c721638e5fa6ded3f943b954 +size 919308 diff --git a/notebooks/bloch_simulation/static_BSSFP_submodules.ipynb b/notebooks/bloch_simulation/static_BSSFP_submodules.ipynb deleted file mode 100644 index 6dd7c66edc7a8ccd5be35517eecfc2f50872402d..0000000000000000000000000000000000000000 --- a/notebooks/bloch_simulation/static_BSSFP_submodules.ipynb +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:78ee646487c91de39dabe4fc4f65feb6271627c392654a05f21c099b546ca101 -size 15550702 diff --git a/notebooks/bloch_simulation/static_multi_coil_bssfp.ipynb b/notebooks/bloch_simulation/static_multi_coil_bssfp.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..fb272988219409e330da1b4f8c8d2051ff42c818 --- /dev/null +++ b/notebooks/bloch_simulation/static_multi_coil_bssfp.ipynb @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e2b13574fc91146c330f7e71c18a632e31e39f054e6960dc0446c60757205696 +size 2106036 diff --git a/notebooks/datasets/regular_grid_dataset.ipynb b/notebooks/datasets/regular_grid_dataset.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..1e4b62e1c44bfb96c7e3db4a5e15efad74f720ce --- /dev/null +++ b/notebooks/datasets/regular_grid_dataset.ipynb @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:389468d74ee83518a1045a2b2bc11189c8c333839623de4bc2d54c31125da0b3 +size 324270 diff --git a/notebooks/local_functions.py b/notebooks/local_functions.py index b71f1db5277a42671d7c76b1aec33c8b22dd4028..a4d93a0351ff94c91eb6b8c2074bfb3b0ee73505 100644 --- a/notebooks/local_functions.py +++ b/notebooks/local_functions.py @@ -1,22 +1,105 @@ +__all__ = ["transformed_box", "create_spherical_phantom", "add_custom_axes", + "add_refilling_box", "animate_trajectories", "get_custom_theme"] from typing import Sequence import pyvista import numpy as np +from pint import Quantity from tqdm.notebook import tqdm -def add_custom_axes(plotter): +def transformed_box(reference_basis: np.ndarray, slice_normal: np.ndarray, + readout_direction: np.ndarray, slice_position: Quantity, + extend: Quantity = Quantity([30, 30, 1], "cm")): + """ + + :param reference_basis: (3x3) matrix where row-index reference the x, y, z + basis vectors + + """ + from scipy.spatial.transform import Rotation + slice_normal /= np.linalg.norm(slice_normal) + + if readout_direction is None: + if np.allclose(slice_normal, np.array([1., 0., 0.]), atol=1e-3): + readout_direction = np.array([0., 1., 0.]) + else: + readout_direction = np.array([1., 0., 0.]) + else: + + if np.isclose( + np.dot(slice_normal, readout_direction / np.linalg.norm(readout_direction)), + 1., atol=1e-3): + raise ValueError("Slice normal and readout-direction are parallel") + + readout_direction = readout_direction - np.dot(slice_normal, + readout_direction) * slice_normal + readout_direction /= np.linalg.norm(readout_direction) + pe_direction = np.cross(readout_direction, slice_normal) + pe_direction /= np.linalg.norm(pe_direction) + new_basis = np.stack([readout_direction, pe_direction, slice_normal]) + rot, _ = Rotation.align_vectors(reference_basis, new_basis) + + total_transform = np.eye(4, 4) + total_transform[:3, :3] = rot.as_matrix() + total_transform[:3, 3] = slice_position.m_as("m") + bounds = np.stack([-extend.m_as("m")/2, extend.m_as("m")/2], axis=1).reshape(-1) + + box = pyvista.Box(bounds=bounds) + box.transform(total_transform, inplace=True) + return box + + +def create_spherical_phantom(dimensions=(120, 120, 120), spacing=(0.002, 0.002, 0.002), + big_radius=0.1) -> pyvista.UnstructuredGrid: + """ Creates a spherical phantom with 6 spherical holes on the coordinate axes and + a rectangular hole in the center. + + :param dimensions: (int, int, int) + :param spacing: (float, float, float) in meter + :param big_radius: float in meter + :return: pyvista.UnstructuredGrid + """ + raster = pyvista.UniformGrid(dimensions=dimensions, spacing=spacing, + origin=-np.array(dimensions) * np.array(spacing) / 2) + sphere_big = pyvista.Sphere(radius=big_radius).extract_surface() + box = pyvista.Box(bounds=(-big_radius * 0.3, big_radius * 0.3, -big_radius/5, + big_radius/5, -big_radius/10, big_radius/10)) + + o = big_radius * 0.7 + spheres_small = [pyvista.Sphere(radius=big_radius*0.05, center=c).extract_surface() for c in + ((-o, 0., 0.), (o, 0., 0.), (0., -o, 0.), (0., o, 0.), (0., 0., -o), + (0., 0., o))] + del o + + filled_sphere = raster.clip_surface(sphere_big, progress_bar=True) + filled_sphere = filled_sphere.clip_surface(box, progress_bar=True, invert=False) + for sms in spheres_small: + filled_sphere = filled_sphere.clip_surface(sms, progress_bar=True, invert=False) + return filled_sphere + + +def add_custom_axes(plotter: pyvista.Plotter, scale: float = 0.1): """ Add scaled arrows to the plotter at origin""" - origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, tip_radius=0.025, - tip_resolution=10, shaft_radius=0.005, shaft_resolution=10, scale=0.1) + origin_arrows = [pyvista.Arrow(start=(0., 0., 0.), direction=direction, tip_length=0.1, + tip_radius=0.025, tip_resolution=10, shaft_radius=0.005, + shaft_resolution=10, scale=scale) for direction in [(1., 0., 0.), (0., 1., 0.), (0., 0., 1.)]] plotter.add_mesh(origin_arrows[0], color="blue", label="x") plotter.add_mesh(origin_arrows[1], color="green", label="y") plotter.add_mesh(origin_arrows[2], color="red", label="z") -def add_refilling_box(plotter, bounding_box, thickness, position): - """ Adds an orange box oriented in (x,y,z) to the plotter, to illustrate a reseeding volume""" +def add_refilling_box(plotter: pyvista.Plotter, bounding_box: Quantity, + thickness: Quantity, position: Quantity): + """ Adds an orange box oriented in (x,y,z) to the plotter, to illustrate a reseeding volume + + :param plotter: - + :param bounding_box: (2, ) + :param thickness: ( ) + :param position: (3, ) + :return: + """ slice_mesh = pyvista.Box(bounds=(-bounding_box[0].m_as("m")/2, bounding_box[0].m_as("m")/2, -bounding_box[1].m_as("m")/2, bounding_box[1].m_as("m") /2, -thickness.m_as("m")/2, thickness.m_as("m")/2)) @@ -24,8 +107,8 @@ def add_refilling_box(plotter, bounding_box, thickness, position): plotter.add_mesh(slice_mesh, color="orange", opacity=0.5, show_edges=True) -def animate_trajectories(plotter, trajectories: Sequence, timing: np.ndarray, - filename: str, scalars: Sequence = None, +def animate_trajectories(plotter: pyvista.Plotter, trajectories: Sequence, + timing: np.ndarray, filename: str, scalars: Sequence = None, vectors: Sequence = None, multi_vecs: bool = False, vector_pos: Sequence = None, @@ -33,66 +116,83 @@ def animate_trajectories(plotter, trajectories: Sequence, timing: np.ndarray, text_kwargs: dict = None, vector_kwargs: dict = None, visualize_trajectories: bool = False, - trajectory_kwargs: dict = None - ): - - if mesh_kwargs is None: - mesh_kwargs = dict() - - if text_kwargs is None: - text_kwargs = dict(color='white', shadow=False, font_size=26, - position='upper_right') - if vector_kwargs is None: - vector_kwargs = dict(mag=1e-3) - - if vector_pos is None: - vector_pos = trajectories - - if trajectory_kwargs is None: - trajectory_kwargs = dict(color="red", show_scalar_bar=False) - - plotter.open_gif(f"{filename}") + trajectory_kwargs: dict = None): + """Animates Sequence of trajectories and saves the animation as gif + + :param plotter: + :param trajectories: sequence containing positions of arbitrary points. Each entry corresponds + to a single snapshot + :param timing: sequence of time-values corresponding to the snapshots in trajectories + :param filename: full path of the saved file (must include .gif ending) + :param scalars: scalar values corresponding to each entry of trajectories (shapes must match) + :param vectors: vectors corresponding to each entry of vector pos (shapes must match) + :param multi_vecs: bool - if true it is assumed that each entry of 'vectors' argument contains + multiple vectors per time-step / per point + :param vector_pos: positions of vectors specified in argument 'vectors' + :param mesh_kwargs: dictionary containing arguments to plotter.add_mesh(.., **mesh_kwargs) + :param text_kwargs: dictionary containing arguments to plotter.add_text(..., **text_kwargs) + :param vector_kwargs: dictionary containing arguments to plotter.add_arrows(..., **vector_kwargs) + :param visualize_trajectories: if true, particle trajectories are traced with a line + :param trajectory_kwargs: dictionary containing argument for lines tracing the trajectories + """ + + if mesh_kwargs is None: + mesh_kwargs = dict() + + if text_kwargs is None: + text_kwargs = dict(color='white', shadow=False, font_size=26, + position='upper_right') + if vector_kwargs is None: + vector_kwargs = dict(mag=1e-3) + + if vector_pos is None: + vector_pos = trajectories + + if trajectory_kwargs is None: + trajectory_kwargs = dict(color="red", show_scalar_bar=False) + + plotter.open_gif(f"{filename}") + if timing is not None: + desc = f"Rendering Mesh-Animation from {timing[0]:2.3f}ms to {timing[-1]:2.3f} ms" + + pbar = tqdm(range(len(trajectories)), desc=desc) + for step in pbar: + plot_mesh = pyvista.PolyData(trajectories[step]) + if timing is not None: - desc = f"Rendering Mesh-Animation from {timing[0]:2.3f}ms to {timing[-1]:2.3f} ms" - - pbar = tqdm(range(len(trajectories)), desc=desc) - for step in pbar: - plot_mesh = pyvista.PolyData(trajectories[step]) - - if timing is not None: - text_actor = plotter.add_text(f"{timing[step]:>2.2f} ms", **text_kwargs) - if scalars is not None: - point_actor = plotter.add_mesh(plot_mesh, scalars=scalars[step], **mesh_kwargs) + text_actor = plotter.add_text(f"{timing[step]:>2.2f} ms", **text_kwargs) + if scalars is not None: + point_actor = plotter.add_mesh(plot_mesh, scalars=scalars[step], **mesh_kwargs) + else: + point_actor = plotter.add_mesh(plot_mesh, **mesh_kwargs) + + if vectors is not None: + if multi_vecs: + arrow_actors = [plotter.add_arrows(cent=vector_pos[step], direction=vectors[i][step], + **vector_kwargs[i]) for i in range(len(vectors))] else: - point_actor = plotter.add_mesh(plot_mesh, **mesh_kwargs) - - if vectors is not None: - if multi_vecs: - arrow_actors = [plotter.add_arrows(cent=vector_pos[step], direction=vectors[i][step], - **vector_kwargs[i]) for i in range(len(vectors))] - else: - arrow_actor = plotter.add_arrows(cent=vector_pos[step], - direction=vectors[step], - **vector_kwargs) - - if visualize_trajectories: - if step >= 1: - plot_mesh["disp"] = trajectories[step]-trajectories[step-1] - glyphs = plot_mesh.glyph(orient="disp", geom=pyvista.Line()) - plotter.add_mesh(glyphs, **trajectory_kwargs) - - plotter.render() - plotter.write_frame() - if timing is not None: - plotter.remove_actor(text_actor) - if vectors is not None: - if multi_vecs: - for act in arrow_actors: - plotter.remove_actor(act) - else: - plotter.remove_actor(arrow_actor) - - plotter.remove_actor(point_actor) + arrow_actor = plotter.add_arrows(cent=vector_pos[step], + direction=vectors[step], + **vector_kwargs) + + if visualize_trajectories: + if step >= 1: + plot_mesh["disp"] = trajectories[step]-trajectories[step-1] + glyphs = plot_mesh.glyph(orient="disp", geom=pyvista.Line()) + plotter.add_mesh(glyphs, **trajectory_kwargs) + + plotter.render() + plotter.write_frame() + if timing is not None: + plotter.remove_actor(text_actor) + if vectors is not None: + if multi_vecs: + for act in arrow_actors: + plotter.remove_actor(act) + else: + plotter.remove_actor(arrow_actor) + + plotter.remove_actor(point_actor) def get_custom_theme(): diff --git a/requirements.txt b/requirements.txt index b6a6a5ecfd7e66d8b1929e41e6aef3862ed1fe26..c2efbd39f879cb8c33d6de4b3d77fd979f10ad46 100644 --- a/requirements.txt +++ b/requirements.txt @@ -8,5 +8,5 @@ tqdm scipy setuptools deprecated -cmrseq>=0.14 +cmrseq>=0.16 diff --git a/setup.py b/setup.py index 716b991a20173fe7388c36d519147b0e28cb9ed4..79ca2eccb4e10b5cf51b3c52c220971242726fe7 100644 --- a/setup.py +++ b/setup.py @@ -16,6 +16,12 @@ author_list = "Jonathan Weine, Charles McGrath" author_email_list = "weine@biomed.ee.ethz.ch, mcgrath@biomed.ee.ethz.ch" url = "https://gitlab.ethz.ch/jweine/cmrsim" +project_urls = { + 'Documentation': 'https://people.ee.ethz.ch/~jweine/cmrsim/latest/index.html', + 'Source': 'https://gitlab.ethz.ch/jweine/cmrsim', + 'Institute': 'https://cmr.ethz.ch/' +} + # Get version tag if '--version' in sys.argv: tag_index = sys.argv.index('--version') + 1 @@ -33,13 +39,14 @@ with open(f'{project_name}/__init__.py', 'r+') as module_header_file: setuptools.setup( name=project_name, url=url, + project_urls=project_urls, version=current_version, author=author_list, author_email_=author_email_list, long_description=long_description, packages=setuptools.find_packages(exclude=["*.tests", "*.tests.*", "tests.*", "tests"]), include_package_data=True, - install_requires=["tensorflow>=2.9", "pyvista>=0.34", "numpy>=0.19", "pint>=0.18", - "Deprecated>=1.2.13"], + install_requires=["tensorflow>=2.11", "pyvista>=0.37", "numpy>=1.19", "pint>=0.18", + "Deprecated>=1.2.13", "cmrseq>=0.17"], python_requires=">=3.6" ) diff --git a/tests/test_bloch_generic.py b/tests/test_bloch_generic.py index ecfc3285751b750b4126c318f424ba1b8a9ff6da..e20cf0ee30730caaa7508b4078e2974f153f1f81 100644 --- a/tests/test_bloch_generic.py +++ b/tests/test_bloch_generic.py @@ -92,47 +92,46 @@ class TestGenericBlochModule(unittest.TestCase): max_grad=Quantity(80., "mT/m"), grad_raster_time=Quantity(0.005, "ms"), rf_raster_time=Quantity(0.005, "ms")) - seq = cmrseq.seqdefs.diffusion.bipolar(system_specs=system_specs, - dt=Quantity(0.1, "ms"), - Dt=Quantity(10., "ms"), - direction=np.array([1., 0., 0.]), - amplitude=Quantity(1., "mT/m"), - rise_time=Quantity(0.03, "ms")) + seq = cmrseq.seqdefs.velocity.bipolar(system_specs=system_specs, + venc=Quantity(200, "cm/s"), + duration=Quantity(2.5, "ms"), + direction=np.array([-1., 0., 0.])) t, gridded_seq = seq.gradients_to_grid() wf = gridded_seq.T dt = system_specs.grad_raster_time.m_as("ms") module = cmrsim.bloch.GeneralBlochOperator("test_grad", - system_specs.gamma_rad.m_as("rad/ms/mT"), + gamma=system_specs.gamma_rad.m_as("rad/ms/mT"), time_grid=t, gradient_waveforms=wf[np.newaxis]) class TrajMod(tf.Module): - def increment_particles(self, r, dt): - v = Quantity([25, 0., 0.], "cm/s") - v = tf.constant(v.m_as("m/ms"), dtype=tf.float32, shape=(1, 3)) - return r + v * dt, dict() + def increment_particles(self, particle_positions: tf.Tensor, dt: tf.Tensor): + v = tf.constant(Quantity([25, 0., 0.], "cm/s").m_as("m/ms"), + dtype=tf.float32, shape=(1, 3)) + return particle_positions + v * dt, dict() + + with self.subTest("Mean Phase increase due to motion in direction of bipolar gradient"): # Set up phantom n_particles = 100 - r_init = tf.random.normal((n_particles, 3), 0., 0.05) + r_init = tf.random.normal((n_particles, 3), 0., 0.005) m = tf.tile(tf.constant([[1., 0., 0.], ], dtype=tf.complex64), [n_particles, 1]) properties = OrderedDict( magnetization=m, T1=tf.ones((n_particles,), dtype=tf.float32) * Quantity(1000, "ms").m, - T2=tf.ones((n_particles,), dtype=tf.float32) * Quantity(300, "ms").m, + T2=tf.ones((n_particles,), dtype=tf.float32) * Quantity(500, "ms").m, M0=tf.ones((n_particles,), dtype=tf.float32) * 10) # call module - m_new, _ = module(trajectory_module=TrajMod(), + m_new, r = module(trajectory_module=TrajMod(), initial_position=r_init, repetition_index=0, run_parallel=True, **properties) - + mean_phase = tf.reduce_mean(tf.math.angle(m_new[0, :, 0])) relaxation_factor = 1. - np.exp(-(wf.shape[0]-1)*dt/1000) # must match T1 definition self.assertGreater(tf.abs(mean_phase), tf.reduce_mean(tf.math.angle(m[:, 0]))) - self.assertTrue(np.allclose(tf.math.angle(m_new[0, :, 0]).numpy(), mean_phase.numpy(), - atol=1e-4)) + self.assertTrue(np.allclose(tf.math.angle(m_new[0, :, 0]).numpy(), mean_phase.numpy(), atol=1e-3)) self.assertTrue(np.allclose(tf.abs(m_new[0, :, 2]), relaxation_factor, atol=1e-4)) diff --git a/tests/test_bloch_multi_coil.py b/tests/test_bloch_multi_coil.py new file mode 100644 index 0000000000000000000000000000000000000000..3f284f9a28b48208d5e6ab31303b73c87e5dd92a --- /dev/null +++ b/tests/test_bloch_multi_coil.py @@ -0,0 +1,95 @@ +import unittest +from collections import OrderedDict +from copy import deepcopy + +import numpy as np +import tensorflow as tf +tf.config.run_functions_eagerly(True) + +from pint import Quantity +import cmrseq + +import cmrsim + + +class TestParallelReceiveBlochSimulation(unittest.TestCase): + def setUp(self) -> None: + self.n_reps = 10 + self.n_steps = 100 + self.time = tf.range(0., self.n_steps, dtype=tf.float32) * 0.01 + self.grads = tf.ones([self.n_reps, self.n_steps, 3], dtype=tf.float32) + self.rf = tf.zeros((self.n_reps, self.n_steps), dtype=tf.complex64) + adc_on = tf.roll(tf.eye(self.n_reps, self.n_steps, + dtype=tf.float32).numpy(), 1, axis=1) + adc_phase = tf.zeros_like(adc_on) + self.adc = tf.stack([adc_on, adc_phase], axis=-1) + + n_particles = 1000 + self.simulation_input = dict( + magnetization=tf.zeros([n_particles, 3], dtype=tf.complex64), + T1=tf.ones([n_particles, ], dtype=tf.float32), + T2=tf.ones([n_particles, ], dtype=tf.float32), + M0=tf.ones([n_particles, ], dtype=tf.float32), + initial_position=tf.random.normal([n_particles, 3], 0., 0.005), + off_res=tf.random.normal([n_particles, 1], 0, 1, dtype=tf.float32), + omega_t2s=tf.random.normal([n_particles, ], 0., 0.1, dtype=tf.float32) + ) + self.n_particles = n_particles + + self.coil_module = cmrsim.analytic.contrast.CoilSensitivity( + coil_sensitivities=np.ones([1, 20, 20, 20, 3], dtype=np.complex64), + map_dimensions=np.array([[-0.1, 0.1], [-0.1, 0.1], [-0.1, 0.1]])) + + def test_input_shapes(self): + + def _test_sequential_call(module: cmrsim.bloch.ParallelReceiveBlochOperator): + m, r = module(repetition_index=tf.constant(0), + run_parallel=False, **self.simulation_input) + self.assertTrue(all([i == j for i, j in zip(m.shape, [self.n_particles, 3])])) + self.assertTrue(m.dtype is tf.complex64) + self.assertTrue(r.dtype is tf.float32) + + def _test_parallel_call(module: cmrsim.bloch.ParallelReceiveBlochOperator): + m, r = module(run_parallel=True, **self.simulation_input) + self.assertTrue(all([i == j for i, j in + zip(m.shape, [self.n_reps, self.n_particles, 3])])) + self.assertTrue(m.dtype is tf.complex64) + self.assertTrue(r.dtype is tf.float32) + + second_iteration_input = deepcopy(self.simulation_input) + [second_iteration_input.pop(k) for k in ("magnetization", "initial_position")] + m, r = module(run_parallel=True, magnetization=m, initial_position=r, + **second_iteration_input) + self.assertTrue(all([i == j for i, j in + zip(m.shape, [self.n_reps, self.n_particles, 3])])) + self.assertTrue(m.dtype is tf.complex64) + self.assertTrue(r.dtype is tf.float32) + + # test __call__ without submodules + module_ = cmrsim.bloch.ParallelReceiveBlochOperator("test", 1., self.coil_module, + time_grid=self.time, + gradient_waveforms=self.grads, + adc_events=self.adc, + rf_waveforms=self.rf) + with self.subTest("No submodules, sequential simulation"): + _test_sequential_call(module_) + + with self.subTest("No submodules, parallel simulation"): + _test_parallel_call(module_) + + # test __call__ with submodules + submodules = (cmrsim.bloch.submodules.OffResonance(1), + cmrsim.bloch.submodules.T2starDistributionModel(), + cmrsim.bloch.submodules.ConcomitantFields(1, 1), + ) + module_ = cmrsim.bloch.ParallelReceiveBlochOperator("test", 1., self.coil_module, + time_grid=self.time, + gradient_waveforms=self.grads, + adc_events=self.adc, + rf_waveforms=self.rf, + submodules=submodules) + with self.subTest("With all submodules sequential call"): + _test_sequential_call(module_) + + with self.subTest("With all submodule parallel call"): + _test_parallel_call(module_) diff --git a/tests/test_datasets.py b/tests/test_datasets.py index e252a8ff3d72b04c6914f5117788d7e7da7e8d0a..1b1ad2036ce3ba622baf2734396359ed0371131f 100644 --- a/tests/test_datasets.py +++ b/tests/test_datasets.py @@ -203,7 +203,6 @@ class TestMeshDatasets(unittest.TestCase): slice_mesh.add_data_per_time(np.random.normal(0., 1., size=(3, npoints, 3)), name="random_vector") - from sys import platform if platform == "linux" or platform == "linux2": pyvista.start_xvfb() @@ -220,3 +219,64 @@ class TestMeshDatasets(unittest.TestCase): mesh_dataset.evaluate_cylinder_coordinates(mesh_dataset.mesh.points, 4, 40*30+1) mesh_dataset.compute_local_basis(time_stamp=Quantity(0., "ms")) mesh_dataset.compute_cylinder_coords(time_stamp=Quantity(0., "ms")) + + +class TestRegularGridDataset(unittest.TestCase): + def test_shepp_logan(self): + dataset = cmrsim.datasets.RegularGridDataset.from_shepp_logan( + map_size=(100, 100, 5), extent=Quantity([10, 10, 5], "cm")) + + def test_from_unstructured_grid(self): + sphere = pyvista.Sphere(radius=0.1, center=(0., 0., 0.)) + ds = cmrsim.datasets.RegularGridDataset.from_unstructured_grid( + sphere, pixel_spacing=Quantity([4., 4., 4.], "mm"), + padding=Quantity([1., 1., 1.], "cm")) + + def test_general_functionality(self): + sphere = pyvista.Sphere(radius=0.1, center=(0., 0., 0.)) + ds = cmrsim.datasets.RegularGridDataset.from_unstructured_grid( + sphere, pixel_spacing=Quantity([4., 4., 4.], "mm"), + padding=Quantity([1., 1., 1.], "cm")) + + ds.mesh["M0"] = np.where(ds.mesh["vtkValidPointMask"], + np.ones_like(ds.mesh["vtkValidPointMask"]), + np.zeros_like(ds.mesh["vtkValidPointMask"])) + + with self.subTest("compute_offres"): + # susceptibilty in ppm + ds.mesh["chi"] = np.where(ds.mesh["vtkValidPointMask"], + np.ones_like(ds.mesh["vtkValidPointMask"]) * (-9.05), + np.ones_like(ds.mesh["vtkValidPointMask"]) * 0.36) * 1e-6 + ds.compute_offresonance(b0=Quantity(3, "T"), susceptibility_key="chi") + + with self.subTest("select_slice"): + slice_normal = np.array([0., 1., 1.]) + readout_direction = np.array([1., 1., 0.]) + slice_position = Quantity([0., 0., 1], "cm") + field_of_view = Quantity([30, 30, 1], "cm") + spacing = Quantity([1, 1, 1], "mm") + + global_slice = ds.select_slice(slice_normal, spacing, slice_position, + field_of_view, readout_direction, in_mps=False) + mps_slice = ds.select_slice(slice_normal, spacing, slice_position, field_of_view, + readout_direction, in_mps=True) + global_slice = ds.select_slice(slice_normal, spacing, slice_position, + field_of_view, in_mps=False) + self.assertRaises(ValueError, lambda: ds.select_slice(slice_normal, spacing, + slice_position, field_of_view, slice_normal, in_mps=False)) + + with self.subTest("get_phantom_dict"): + ds.mesh["is_non_trivial"] = ds.mesh["M0"] > 0. + phantom_dict = ds.get_phantom_def(filter_by="is_non_trivial", keys=["M0", "offres"], + perturb_positions_std=Quantity(0.1, "mm").m_as("m")) + + with self.subTest("add_field"): + dat = np.ones(ds.mesh.dimensions) + ds.add_field(key="dummy", data=dat) + self.assertTrue("dummy" in ds.mesh.array_names) + + with self.subTest("resample_spacing"): + ds.resample_spacing(spacing=Quantity([3, 3, 3], "mm")) + + + diff --git a/tests/test_signal_so_phasetracking.py b/tests/test_signal_so_phasetracking.py index 6864d00cdd96dbe8efb8771d49a21fb328d17979..51968eb055c01524b0ab39380c6dcc2a2ee8722a 100644 --- a/tests/test_signal_so_phasetracking.py +++ b/tests/test_signal_so_phasetracking.py @@ -26,7 +26,9 @@ class TestDiffusionWaveFormTracking(unittest.TestCase): self.system_specs = cmrseq.SystemSpec(max_grad=Quantity(80., "mT/m"), max_slew=Quantity(80., "mT/m/ms"), grad_raster_time=Quantity(0.01, "ms")) - seq = cmrseq.seqdefs.diffusion.m012(self.system_specs, Quantity(100., "s/mm**2"), + seq = cmrseq.seqdefs.diffusion.m012(self.system_specs, lambda_=Quantity(5.6, "ms"), + zeta=Quantity(0.9, "ms"), + bvalue=Quantity(100., "s/mm**2"), direction=np.array([0., 1., 0.])) t, wf = seq.gradients_to_grid() wf = np.concatenate([t[np.newaxis], wf], axis=0).T @@ -53,7 +55,6 @@ class TestDiffusionWaveFormTracking(unittest.TestCase): signal_tensor = np.ones((rad_cords.shape[0], 1, 1), dtype=np.complex64) signal_tensor = self.mod(signal_tensor, trajectory) phase = np.angle(signal_tensor) - print(phase) self.assertTrue(np.allclose(phase, 0., atol=np.pi/1000), msg="Sum of all phases is greater than numerical precision " "tolerance for static phantom") diff --git a/tests/test_signal_so_sequences.py b/tests/test_signal_so_sequences.py index a7e6b3f212b4bbc2ffbcd99d47a7517b5e030c8b..1b59a3b4eb2fdf0d89dd8fcc49c78e56c0e1d7ea 100644 --- a/tests/test_signal_so_sequences.py +++ b/tests/test_signal_so_sequences.py @@ -86,18 +86,20 @@ class TestBSSFP(unittest.TestCase): T2star = tf.ones(shape=(10, 1, 1), dtype=tf.float32) offres = tf.zeros(shape=(10, 1, 1), dtype=tf.float32) - a = bssfp_module(m0, T1=T1, T2=T2, T2star=T2star, theta=offres) + a = bssfp_module(m0, T1=T1, T2=T2, T2star=T2star, off_res=offres) def test_banding(self): - theta_grid = tf.range(-2 * np.pi, 2 * np.pi, np.pi / 50, dtype=tf.float32) + TR = 6 # ms + theta_grid = tf.range(-2 * np.pi, 2 * np.pi, np.pi / 50, dtype=tf.float32) / TR + sequence_module = BSSFP(flip_angle=np.array([35, ]), echo_time=[3, ], - repetition_time=[6, ], expand_repetitions=True) - inp1 = dict(signal_tensor=tf.ones([1, 1, 1], dtype=tf.complex64), T1=tf.ones([1, 1, 1], dtype=tf.float32) * 600, + repetition_time=[TR, ], expand_repetitions=True) + inp1 = dict(signal_tensor=tf.ones([1, 1, 1], dtype=tf.complex64), T1=tf.ones([1, 1, 1] ,dtype=tf.float32) * 600, T2=tf.ones([1, 1, 1], dtype=tf.float32) * 100, T2star=tf.ones([1, 1, 1], dtype=tf.float32) * 1000, - theta=tf.reshape(theta_grid, [1, -1, 1])) + off_res=tf.reshape(theta_grid, [1, -1, 1])) inp2 = dict(signal_tensor=tf.ones([1, 1, 1], dtype=tf.complex64), T1=tf.ones([1, 1, 1], dtype=tf.float32) * 900, T2=tf.ones([1, 1, 1], dtype=tf.float32) * 50, T2star=tf.ones([1, 1, 1], dtype=tf.float32) * 1000, - theta=tf.reshape(theta_grid, [1, -1, 1])) + off_res=tf.reshape(theta_grid, [1, -1, 1])) res1 = sequence_module(**inp1) res2 = sequence_module(**inp2)