From 148ef91311d5deab0d5a28d1282c24302322a63e Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Fri, 28 Nov 2025 09:30:15 +0100 Subject: [PATCH 01/32] add EpiReadout (work in progress) --- .pre-commit-config.yaml | 1 + src/mrseq/utils/trajectory.py | 309 ++++++++++++++++++++++++++++++++++ 2 files changed, 310 insertions(+) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 87dd6b9..ee54a49 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -48,6 +48,7 @@ repos: - types-requests - typing-extensions - pytest + - matplotlib ci: autofix_commit_msg: | diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py index 5859a64..5836e3f 100644 --- a/src/mrseq/utils/trajectory.py +++ b/src/mrseq/utils/trajectory.py @@ -1,10 +1,14 @@ """Basic functionality for trajectory calculation.""" +import warnings from typing import Literal +import matplotlib.pyplot as plt import numpy as np import pypulseq as pp +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults from mrseq.utils import variable_density_spiral_trajectory @@ -302,3 +306,308 @@ def combine_gradients(*grad_objects, channel): time_to_echo = max_pre_duration + n_samples_to_echo * readout_oversampling * adc.dwell return gx_combined, gy_combined, adc, trajectory, time_to_echo + + +class EpiReadout: + """EPI readout module. + + Attributes + ---------- + system + System limits. + fov + Field of view in meters (square). + n_readout + Number of readout points. + n_phase_encoding + Number of phase encoding points. + bandwidth + Total receiver bandwidth in Hz (approx. 1/dwell_time_nyquist). + readout_duration will be approx matrix_size[0] / bandwidth. + oversampling + ADC oversampling factor. + ramp_sampling + If True, ADC is active during gradient ramps (optimized timing). + readout_type + Readout type ('flyback' or 'symmetric'). + partial_fourier_factor + Partial Fourier factor (0.5 to 1.0). + adc_freq_offset + Frequency offset for the ADC. + pe_enable + Enable phase encoding (useful for calibration if False). + spoiling_enable + Enable spoiling gradients (useful for calibration if False). + adc + ADC event. + gx_pre + Pre-phaser gradient in readout direction. + gy_pre + Pre-phaser gradient in phase encoding direction. + gx + Readout gradient. + gy + Phase encoding gradient. + gy_blip + Complete blip gradient used for flyback readout. + gy_blipup + Ramp-up part of blip gradient used for symmetric readout. + gy_blipdown + Ramp-down part of blip gradient used for symmetric readout. + gy_blipdownup + Combined gy_blipdown and gy_blipup gradient used for symmetric readout. + gz_spoil + Spoiler gradient in z-direction. + gx_flyback + Flyback gradient in x-direction. + """ + + def __init__( + self, + system: pp.Opts | None = None, + fov: float = 128e-3, + n_readout: int = 64, + n_phase_encoding: int = 64, + bandwidth: float = 50e3, + oversampling: int = 2, + ramp_sampling: bool = True, + readout_type: Literal['flyback', 'symmetric'] = 'symmetric', + partial_fourier_factor: float = 1.0, + adc_freq_offset: float = 0, + pe_enable: bool = True, + spoiling_enable: bool = True, + ): + """Initialize EPI Readout. + + Parameters + ---------- + system + System limits. + fov + Field of view in meters (square). + n_readout + Number of readout points. + n_phase_encoding + Number of phase encoding points. + bandwidth + Total receiver bandwidth in Hz (approx. 1/dwell_time_nyquist). + readout_duration will be approx matrix_size[0] / bandwidth. + oversampling + ADC oversampling factor. + ramp_sampling + If True, ADC is active during gradient ramps (optimized timing). + readout_type + Readout type ('flyback' or 'symmetric'). + partial_fourier_factor + Partial Fourier factor (0.5 to 1.0]. + adc_freq_offset + Frequency offset for the ADC. + pe_enable + Enable phase encoding (useful for calibration if False). + spoiling_enable + Enable spoiling gradients (useful for calibration if False). + """ + if not 0.5 < partial_fourier_factor <= 1.0: + raise ValueError('Desired partial Fourier factor must be larger than 0.5 and smaller or equal to 1.0.') + + # set system to default if not provided + if system is None: + system = sys_defaults + + self.system = system + self.fov = fov + self.n_readout = n_readout + self.n_phase_encoding = n_phase_encoding + self.bandwidth = bandwidth + self.oversampling = oversampling + self.ramp_sampling = ramp_sampling + self.readout_type = readout_type + self.partial_fourier_factor = partial_fourier_factor + self.adc_freq_offset = adc_freq_offset + self.pe_enable = pe_enable + self.spoiling_enable = spoiling_enable + + # Derived parameters + self.readout_time = self.n_readout / self.bandwidth + self.delta_k = 1 / self.fov + + # Initiate all optional gradients as None + self.gx_flyback = None + self.gy_blipdown = None + self.gy_blipup = None + self.gy_blipdownup = None + self.gz_spoil = None + + # Create blip gradient with shortest possible timing + gy_blip_duration = np.ceil(2 * np.sqrt(self.delta_k / system.max_slew) / 10e-6 / 2) * 10e-6 * 2 + gy_blip_half_dur = gy_blip_duration / 2 + self.gy_blip = pp.make_trapezoid(channel='y', system=self.system, area=-self.delta_k, duration=gy_blip_duration) + + # Create readout gradient + gx_encoding_area = self.n_readout * self.delta_k + if self.ramp_sampling: + # Calculate additional gradient area from gy_blip assuming maximum slew rate + gy_blip_area = np.power(gy_blip_half_dur, 2) * self.system.max_slew + + # Create gradient with additional area + gx = pp.make_trapezoid( + channel='x', + system=self.system, + area=gx_encoding_area + gy_blip_area, + duration=gy_blip_half_dur + self.readout_time + gy_blip_half_dur, + ) + + if not gx.fall_time == gx.rise_time: + raise ValueError('Gradient fall time must be equal to rise time for ramp sampling.') + + # Second, correct area taking actual slew rate into account + gx_slew = gx.amplitude / gx.rise_time + gx_area_reduced_slew = gx.area - np.power(gy_blip_half_dur, 2) * gx_slew + + gx.amplitude = float(gx.amplitude * gx_encoding_area / gx_area_reduced_slew) + gx.area = float(gx.amplitude * (gx.rise_time / 2 + gx.flat_time + gx.fall_time / 2)) + gx.flat_area = float(gx.amplitude * gx.flat_time) + self.gx = gx + else: + self.gx = pp.make_trapezoid( + channel='x', + system=self.system, + flat_area=gx_encoding_area, + flat_time=self.readout_time, + ) + + # Create ADC event + adc_dwell = self.delta_k / self.gx.amplitude / self.oversampling + adc_dwell = round_to_raster(adc_dwell, self.system.adc_raster_time, 'ceil') + adc_samples = int(np.ceil(self.readout_time / adc_dwell / 4) * 4) + + adc_time_to_center = adc_dwell * ((adc_samples - 1) / 2 + 0.5) + adc_delay = self.gx.rise_time + self.gx.flat_time / 2 - adc_time_to_center + adc_delay = round_to_raster(adc_delay, self.system.adc_raster_time) + + self.adc = pp.make_adc( + num_samples=adc_samples, + dwell=adc_dwell, + freq_offset=self.adc_freq_offset, + delay=adc_delay, + system=self.system, + ) + + # Create and align pre-phaser gradients considering partial fourier factor + # determine the number of pe lines after (and including) k-space center, which is independent of partial fourier + self.n_phase_enc_post_center = int(np.ceil(self.n_phase_encoding / 2 + 1)) + # find the closest number of lines to the desired partial fourier factor + valid_n_phase_total = int(np.ceil(self.partial_fourier_factor * self.n_phase_encoding)) + # ensure that at least one line before the center is acquired + self.n_phase_enc_pre_center = max(1, valid_n_phase_total - self.n_phase_enc_post_center) + self.n_phase_enc_total = self.n_phase_enc_pre_center + self.n_phase_enc_post_center + # recalculate the actual partial fourier factor + actual_pf_factor = self.n_phase_enc_total / self.n_phase_encoding + if actual_pf_factor != self.partial_fourier_factor: + warnings.warn( + f'Desired partial Fourier factor {self.partial_fourier_factor} adjusted to {actual_pf_factor}.', + UserWarning, + stacklevel=2, + ) + self.partial_fourier_factor = actual_pf_factor + + self.gx_pre = pp.make_trapezoid(channel='x', system=self.system, area=-self.gx.area / 2 - self.delta_k / 2) + self.gy_pre = pp.make_trapezoid( + channel='y', system=self.system, area=self.n_phase_enc_pre_center * self.delta_k + ) + self.gx_pre, self.gy_pre = pp.align(right=[self.gx_pre, self.gy_pre]) + + if self.readout_type == 'flyback': + # Create flyback gradient in case of flyback readout + self.gx_flyback = pp.make_trapezoid(channel='x', system=self.system, area=-self.gx.area) + elif self.readout_type == 'symmetric': + # Split and align blip gradient in case of symmetric readout + gy_parts = pp.split_gradient_at(grad=self.gy_blip, time_point=gy_blip_duration / 2, system=self.system) + self.gy_blipup, self.gy_blipdown, _ = pp.align(right=gy_parts[0], left=[gy_parts[1], self.gx]) + self.gy_blipdownup = pp.add_gradients((self.gy_blipdown, self.gy_blipup), system=self.system) + else: + raise NotImplementedError('Currently, only "symmetric" and "flyback" readout types are supported.') + + # Disable phase encoding if self.pe_enable is False + if not self.pe_enable: + self.gy_pre.amplitude = 0 + if ( + self.readout_type == 'symmetric' + and self.gy_blipup is not None + and self.gy_blipdown is not None + and self.gy_blipdownup is not None + ): + self.gy_blipup.waveform *= 0 + self.gy_blipdown.waveform *= 0 + self.gy_blipdownup.waveform *= 0 + + # Create spoiler gradient if spoiling is enabled + if self.spoiling_enable: + self.gz_spoil = pp.make_trapezoid(channel='z', system=self.system, area=4 * self.delta_k * n_readout) + + def add_to_seq(self, seq: pp.Sequence): + """Add EPI readout blocks to the sequence.""" + # (Re)set phase encoding (LIN) label + lin_label = pp.make_label(label='LIN', type='SET', value=0) + + # Add pre-phaser gradients + seq.add_block(self.gx_pre, self.gy_pre) + + for pe_idx in range(self.n_phase_enc_total): + rev_label = pp.make_label(type='SET', label='REV', value=self.gx.amplitude < 0) + if self.readout_type == 'symmetric': + # Select blip gradient based on phase encoding index + if pe_idx == 0: + gy_blip = self.gy_blipup + elif pe_idx == self.n_phase_enc_total - 1: + gy_blip = self.gy_blipdown + else: + gy_blip = self.gy_blipdownup + + # Add readout block and reverse polarity of readout gradient + seq.add_block(self.gx, gy_blip, self.adc, lin_label, rev_label) + self.gx.amplitude = -self.gx.amplitude + elif self.readout_type == 'flyback': + seq.add_block(self.gx, self.adc, lin_label, rev_label) + if pe_idx != self.n_phase_enc_total - 1: + seq.add_block(self.gx_flyback, self.gy_blip) + + lin_label = pp.make_label(label='LIN', type='INC', value=1) + + if self.spoiling_enable: + seq.add_block(self.gz_spoil) + + return seq + + def plot_sequence(self): + """Plot the sequence.""" + seq = pp.Sequence(self.system) + seq = self.add_to_seq(seq) + seq.plot(grad_disp='mT/m') + + def plot_trajectory(self): + """Plot k-space trajectory.""" + seq = pp.Sequence(self.system) + # add dummy excitation pulse for trajectory calculation + seq.add_block( + pp.make_block_pulse( + flip_angle=np.pi / 2, + duration=2e-3, + delay=self.system.rf_dead_time, + use='excitation', + system=self.system, + ) + ) + seq = self.add_to_seq(seq) + + # calculate k-space trajectory + k_traj_adc, k_traj, _, _, _ = seq.calculate_kspace() + + # plot trajectory + fig = plt.figure() + plt.plot(k_traj[0], k_traj[1], 'b') + plt.plot(k_traj_adc[0], k_traj_adc[1], 'x', color='red', markersize=4) + plt.grid() + plt.show() + + return fig From a7f092c4cb317013e58108690499dc0880f8a43c Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Mon, 15 Dec 2025 17:43:59 +0100 Subject: [PATCH 02/32] ensure round_to_raster returns floats --- src/mrseq/utils/sequence_helper.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mrseq/utils/sequence_helper.py b/src/mrseq/utils/sequence_helper.py index 2584bae..53458e5 100644 --- a/src/mrseq/utils/sequence_helper.py +++ b/src/mrseq/utils/sequence_helper.py @@ -23,11 +23,11 @@ def round_to_raster(value: float, raster_time: float, method: Literal['floor', ' Rounded value. """ if method == 'floor': - return raster_time * np.floor(value / raster_time) + return float(raster_time * np.floor(value / raster_time)) elif method == 'round': - return raster_time * np.round(value / raster_time) + return float(raster_time * np.round(value / raster_time)) elif method == 'ceil': - return raster_time * np.ceil(value / raster_time) + return float(raster_time * np.ceil(value / raster_time)) else: raise ValueError(f'Unknown rounding method: {method}. Expected: "floor", "round" or "ceil".') From 9b7a9fa84c09cc8c318ac9bc2096744492cf9dbd Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Mon, 15 Dec 2025 17:53:04 +0100 Subject: [PATCH 03/32] update EpiReadout and add epi script (WIP) --- src/mrseq/scripts/epi2d.py | 383 ++++++++++++++++++++++++++++++++++ src/mrseq/utils/trajectory.py | 166 ++++++++++----- 2 files changed, 494 insertions(+), 55 deletions(-) create mode 100644 src/mrseq/scripts/epi2d.py diff --git a/src/mrseq/scripts/epi2d.py b/src/mrseq/scripts/epi2d.py new file mode 100644 index 0000000..65b1a98 --- /dev/null +++ b/src/mrseq/scripts/epi2d.py @@ -0,0 +1,383 @@ +"""2D Echo Planar Imaging (EPI) sequence.""" + +from pathlib import Path +from typing import Literal + +import ismrmrd +import numpy as np +import pypulseq as pp + +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults +from mrseq.utils.ismrmrd import Fov +from mrseq.utils.ismrmrd import Limits +from mrseq.utils.ismrmrd import MatrixSize +from mrseq.utils.ismrmrd import create_header +from mrseq.utils.trajectory import EpiReadout + + +def epi2d_kernel( + system: pp.Opts, + te: float | None, + tr: float | None, + fov: float, + n_readout: int, + n_phase_encoding: int, + bandwidth: float, + slice_thickness: float, + n_slices: int, + rf_duration: float, + rf_flip_angle: float, + rf_bwt: float, + rf_apodization: float, + readout_type: Literal['symmetric', 'flyback'], + oversampling: Literal[1, 2, 4], + ramp_sampling: bool, + partial_fourier_factor: float, + pe_enable: bool, + spoiling_enable: bool, + mrd_header_file: str | Path | None, +) -> tuple[pp.Sequence, float, float]: + """Generate a 2D Echo Planar Imaging (EPI) sequence. + + Parameters + ---------- + system + PyPulseq system limits object. + te + Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. + tr + Desired repetition time (TR) (in seconds). + fov + Field of view in x and y direction (in meters). + n_readout + Number of frequency encoding steps. + n_phase_encoding + Number of phase encoding steps. + bandwidth + Total receiver bandwidth in Hz. + slice_thickness + Slice thickness of the 2D slice (in meters). + n_slices + Number of slices. + rf_duration + Duration of the rf excitation pulse (in seconds) + rf_flip_angle + Flip angle of rf excitation pulse (in degrees) + rf_bwt + Bandwidth-time product of rf excitation pulse (Hz * seconds) + rf_apodization + Apodization factor of rf excitation pulse + readout_type + Readout type ('symmetric' or 'flyback'). + oversampling + ADC oversampling factor. + ramp_sampling + If True, ADC is active during gradient ramps (optimized timing). + partial_fourier_factor + Partial Fourier factor (0.5 to 1.0). + pe_enable + Enable phase encoding (useful for calibration scans if False). + spoiling_enable + Enable spoiling gradients. + mrd_header_file + Filename of the ISMRMRD header file to be created. If None, no header file is created. + + Returns + ------- + seq + PyPulseq Sequence object + min_te + Shortest possible echo time. + min_tr + Shortest possible repetition time. + + """ + # create PyPulseq Sequence object and set system limits + seq = pp.Sequence(system=system) + + # create EpiReadout object + epi2d = EpiReadout( + system=system, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=bandwidth, + oversampling=oversampling, + readout_type=readout_type, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, + pe_enable=pe_enable, + spoiling_enable=spoiling_enable, + ) + + # create slice selective excitation pulse and gradients + rf, gz, gzr = pp.make_sinc_pulse( + flip_angle=np.deg2rad(rf_flip_angle), + duration=rf_duration, + slice_thickness=slice_thickness, + apodization=rf_apodization, + time_bw_product=rf_bwt, + delay=system.rf_dead_time, + system=system, + return_gz=True, + use='excitation', + ) + + # calculate minimum echo time + gzr_prephaser_dur = pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gx_pre) + min_te = rf.shape_dur / 2 # time from center to end of RF pulse + min_te += max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time + min_te += gzr_prephaser_dur # for minimum TE, gzr and pre-phasers are played out simultaneously + min_te += epi2d.time_to_center_without_prephaser + + # calculate echo time delay (te_delay) + if te is None: + te_delay = 0.0 + else: + te_delay = round_to_raster(te - min_te, system.block_duration_raster) + if te_delay < 0: + raise ValueError(f'TE must be larger than {min_te * 1000:.3f} ms. Current value is {te * 1000:.3f} ms.') + + # calculate minimum repetition time depending on chosen echo time + min_tr = pp.calc_duration(rf, gz) + min_tr += gzr_prephaser_dur # if TE is chosen minimal, gzr and pre-phasers are played out simultaneously + min_tr += epi2d.total_duration_without_prephaser + + # calculate repetition time delay (tr_delay) for current TE settings + current_min_tr = min_tr + te_delay + if tr is None: + tr_delay = 0.0 + else: + tr_delay = round_to_raster(tr - current_min_tr, system.block_duration_raster) + if tr_delay < 0: + raise ValueError( + f'TR must be larger than {current_min_tr * 1000:.2f} ms. Current value is {tr * 1000:.3f} ms.' + ) + + print(f'\nCurrent echo time = {(min_te + te_delay) * 1000:.3f} ms') + print(f'Current repetition time = {(current_min_tr + tr_delay) * 1000:.3f} ms') + + # create header + if mrd_header_file: + hdr = create_header( + traj_type='other', + encoding_fov=Fov(x=fov * oversampling, y=fov, z=slice_thickness), + recon_fov=Fov(x=fov, y=fov, z=slice_thickness), + encoding_matrix=MatrixSize(n_x=n_readout, n_y=epi2d.n_phase_encoding, n_z=1), + recon_matrix=MatrixSize(n_x=n_readout, n_y=n_phase_encoding, n_z=1), + dwell_time=epi2d.adc.dwell, + slice_limits=Limits(min=0, max=n_slices, center=n_slices // 2), + k1_limits=Limits(min=0, max=epi2d.n_phase_enc_total, center=epi2d.n_phase_enc_pre_center + 1), + k2_limits=Limits(min=0, max=1, center=0), + ) + + # write header to file + prot = ismrmrd.Dataset(mrd_header_file, 'w') + prot.write_xml_header(hdr.toXML('utf-8')) + + # obtain noise samples + seq.add_block( + pp.make_label(label='LIN', type='SET', value=0), + pp.make_label(label='SLC', type='SET', value=0), + pp.make_label(label='NOISE', type='SET', value=True), + ) + seq.add_block( + epi2d.adc, pp.make_delay(round_to_raster(pp.calc_duration(epi2d.adc), system.block_duration_raster, 'ceil')) + ) + seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) + seq.add_block(pp.make_delay(system.rf_dead_time)) + + if mrd_header_file: + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) + prot.append_acquisition(acq) + + for slice_ in range(n_slices): + # define label(s) + slice_label = pp.make_label(label='SLC', type='SET', value=slice_) + + # set frequency offset for current slice + rf.freq_offset = gz.amplitude * slice_thickness * (slice_ - (n_slices - 1) / 2) + + # add slice selective excitation pulse and set slice label + seq.add_block(rf, gz, slice_label) + + # add echo time delay + if te_delay > 0: + seq.add_block(pp.make_delay(te_delay)) + + # add slice selection rephaser and readout prephaser gradients + gzr, gx_pre, gy_pre = pp.align(left=[gzr], right=[epi2d.gx_pre, epi2d.gy_pre]) + seq.add_block(gzr, gx_pre, gy_pre) + + # add EPI readout block without pre-phaser gradients + seq, prot = epi2d.add_to_seq(seq, add_prephaser=False, mrd_dataset=prot) + + # add repetition time delay + if tr_delay > 0: + seq.add_block(pp.make_delay(tr_delay)) + + # calculate k-space trajectory from Sequence to add this info to mrd file + k_traj_adc, _, _, _, _ = seq.calculate_kspace() + samples_per_acq = epi2d.adc.num_samples + number_of_total_acq = k_traj_adc.shape[-1] // samples_per_acq + number_of_epi_acq = epi2d.n_phase_enc_total + if 'data' in prot._dataset: + number_of_noise_acq = prot.number_of_acquisitions() + else: + number_of_noise_acq = 0 + + if not number_of_epi_acq + number_of_noise_acq == number_of_total_acq: + raise (ValueError('Number of calculated acquisitions does not match expected number.')) + + k_traj_adc_readout = k_traj_adc[:, number_of_noise_acq * samples_per_acq :] + + # create mrd acquisition and add trajectory info for each EPI readout + for n in range(number_of_epi_acq): + start = n * samples_per_acq + end = start + samples_per_acq + + traj = np.zeros((samples_per_acq, 2), dtype=np.float32) + traj[:, 0] = np.round(k_traj_adc_readout[0, start:end] * fov, 3) + traj[:, 1] = np.round(k_traj_adc_readout[1, start:end] * fov, 3) + + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=samples_per_acq) + acq.traj[:] = traj + + prot.append_acquisition(acq) + + # close ISMRMRD file + if mrd_header_file: + prot.close() + + # set gridding definitions extracted from EpiReadout + if ramp_sampling: + gridding_params = [ + epi2d.gx.rise_time, + epi2d.gx.flat_time, + epi2d.gx.fall_time, + epi2d.adc.delay - epi2d.gx.delay, + epi2d.adc.num_samples * epi2d.adc.dwell, + ] + + seq.set_definition(key='TargetGriddedSamples', value=n_readout * oversampling) + seq.set_definition(key='TrapezoidGriddingParameters', value=gridding_params) + + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov, fov, slice_thickness * n_slices]) + seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, n_slices)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr or min_tr) + seq.set_definition('ReadoutOversamplingFactor', oversampling) + + return seq, min_te, min_tr + + +def main( + system: pp.Opts | None = None, + te: float | None = None, + tr: float | None = None, + fov: float = 256e-3, + n_readout: int = 64, + n_phase_encoding: int = 64, + n_slices: int = 1, + slice_thickness: float = 8e-3, + bandwidth: float = 64e3, + readout_type: Literal['symmetric', 'flyback'] = 'flyback', + oversampling: Literal[1, 2, 4] = 1, + ramp_sampling: bool = False, + partial_fourier_factor: float = 1, + show_plots: bool = True, + test_report: bool = True, + timing_check: bool = True, +) -> tuple[pp.Sequence, Path]: + """Generate an Echo-Planar Imaging (EPI) sequence. + + Returns + ------- + seq + Sequence object of radial FLASH sequence. + file_path + Path to the sequence file. + """ + if system is None: + system = sys_defaults + + # define settings of rf excitation pulse + rf_flip_angle: float = 90 # flip angle of the rf excitation pulse [°] + rf_duration = 2.56e-3 # duration of the rf excitation pulse [s] + rf_bwt = 4 # bandwidth-time product of rf excitation pulse [Hz*s] + rf_apodization = 0.5 # apodization factor of rf excitation pulse + + # define EPI settings + enable_phase_encoding = True + enable_gradient_spoiling = True + + # define sequence filename + rs_string = 'rs' if ramp_sampling else 'nors' + pf_string = f'{partial_fourier_factor}'.replace('.', 'p') + + filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny' + filename += f'_{readout_type}_{oversampling}ro_{rs_string}_{pf_string}' + + output_path = Path.cwd() / 'output' + output_path.mkdir(parents=True, exist_ok=True) + + # delete existing header file + if (output_path / Path(filename + '_header.h5')).exists(): + (output_path / Path(filename + '_header.h5')).unlink() + + mrd_file = output_path / Path(filename + '_header.h5') + + seq, _min_te, min_tr = epi2d_kernel( + system=system, + te=te, + tr=tr, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=bandwidth, + oversampling=oversampling, + slice_thickness=slice_thickness, + n_slices=n_slices, + rf_duration=rf_duration, + rf_flip_angle=rf_flip_angle, + rf_bwt=rf_bwt, + rf_apodization=rf_apodization, + readout_type=readout_type, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, + pe_enable=enable_phase_encoding, + spoiling_enable=enable_gradient_spoiling, + mrd_header_file=mrd_file, + ) + + # check timing of the sequence + if timing_check and not test_report: + ok, error_report = seq.check_timing() + if ok: + print('\nTiming check passed successfully') + else: + print('\nTiming check failed! Error listing follows\n') + print(error_report) + + # show advanced rest report + if test_report: + print('\nCreating advanced test report...') + print(seq.test_report()) + + # save seq-file to disk + print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") + seq.write(str(output_path / filename), create_signature=True) + + if show_plots: + seq.plot(time_range=(0, 10 * (tr or min_tr))) + + return seq, output_path / filename + + +if __name__ == '__main__': + main() diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py index 5836e3f..63c33ab 100644 --- a/src/mrseq/utils/trajectory.py +++ b/src/mrseq/utils/trajectory.py @@ -3,6 +3,7 @@ import warnings from typing import Literal +import ismrmrd import matplotlib.pyplot as plt import numpy as np import pypulseq as pp @@ -309,33 +310,16 @@ def combine_gradients(*grad_objects, channel): class EpiReadout: - """EPI readout module. + """EPI readout module supporting flyback and symmetric readout, ramp sampling, oversampling, and partial Fourier. Attributes ---------- system System limits. - fov - Field of view in meters (square). - n_readout - Number of readout points. - n_phase_encoding - Number of phase encoding points. - bandwidth - Total receiver bandwidth in Hz (approx. 1/dwell_time_nyquist). - readout_duration will be approx matrix_size[0] / bandwidth. - oversampling - ADC oversampling factor. ramp_sampling If True, ADC is active during gradient ramps (optimized timing). readout_type Readout type ('flyback' or 'symmetric'). - partial_fourier_factor - Partial Fourier factor (0.5 to 1.0). - adc_freq_offset - Frequency offset for the ADC. - pe_enable - Enable phase encoding (useful for calibration if False). spoiling_enable Enable spoiling gradients (useful for calibration if False). adc @@ -390,8 +374,7 @@ def __init__( n_phase_encoding Number of phase encoding points. bandwidth - Total receiver bandwidth in Hz (approx. 1/dwell_time_nyquist). - readout_duration will be approx matrix_size[0] / bandwidth. + Total receiver bandwidth in Hz. oversampling ADC oversampling factor. ramp_sampling @@ -399,13 +382,13 @@ def __init__( readout_type Readout type ('flyback' or 'symmetric'). partial_fourier_factor - Partial Fourier factor (0.5 to 1.0]. + Partial Fourier factor (0.5 to 1.0). adc_freq_offset Frequency offset for the ADC. pe_enable - Enable phase encoding (useful for calibration if False). + Enable phase encoding (useful for calibration scans if False). spoiling_enable - Enable spoiling gradients (useful for calibration if False). + Enable spoiling gradients. """ if not 0.5 < partial_fourier_factor <= 1.0: raise ValueError('Desired partial Fourier factor must be larger than 0.5 and smaller or equal to 1.0.') @@ -428,8 +411,9 @@ def __init__( self.spoiling_enable = spoiling_enable # Derived parameters - self.readout_time = self.n_readout / self.bandwidth - self.delta_k = 1 / self.fov + readout_time = n_readout / bandwidth + delta_kx = 1 / (fov * oversampling) + delta_ky = 1 / fov # Initiate all optional gradients as None self.gx_flyback = None @@ -439,22 +423,22 @@ def __init__( self.gz_spoil = None # Create blip gradient with shortest possible timing - gy_blip_duration = np.ceil(2 * np.sqrt(self.delta_k / system.max_slew) / 10e-6 / 2) * 10e-6 * 2 + gy_blip_duration = np.ceil(2 * np.sqrt(delta_ky / system.max_slew) / 10e-6 / 2) * 10e-6 * 2 gy_blip_half_dur = gy_blip_duration / 2 - self.gy_blip = pp.make_trapezoid(channel='y', system=self.system, area=-self.delta_k, duration=gy_blip_duration) + self.gy_blip = pp.make_trapezoid(channel='y', system=self.system, area=-delta_ky, duration=gy_blip_duration) # Create readout gradient - gx_encoding_area = self.n_readout * self.delta_k + gx_encoding_area = n_readout * delta_kx * oversampling if self.ramp_sampling: # Calculate additional gradient area from gy_blip assuming maximum slew rate - gy_blip_area = np.power(gy_blip_half_dur, 2) * self.system.max_slew + extra_area = np.power(gy_blip_half_dur, 2) * self.system.max_slew # Create gradient with additional area gx = pp.make_trapezoid( channel='x', system=self.system, - area=gx_encoding_area + gy_blip_area, - duration=gy_blip_half_dur + self.readout_time + gy_blip_half_dur, + area=gx_encoding_area + extra_area, + duration=gy_blip_half_dur + readout_time + gy_blip_half_dur, ) if not gx.fall_time == gx.rise_time: @@ -462,7 +446,7 @@ def __init__( # Second, correct area taking actual slew rate into account gx_slew = gx.amplitude / gx.rise_time - gx_area_reduced_slew = gx.area - np.power(gy_blip_half_dur, 2) * gx_slew + gx_area_reduced_slew = gx.area - gx_slew * np.power(gy_blip_half_dur, 2) gx.amplitude = float(gx.amplitude * gx_encoding_area / gx_area_reduced_slew) gx.area = float(gx.amplitude * (gx.rise_time / 2 + gx.flat_time + gx.fall_time / 2)) @@ -473,52 +457,62 @@ def __init__( channel='x', system=self.system, flat_area=gx_encoding_area, - flat_time=self.readout_time, + flat_time=readout_time, ) # Create ADC event - adc_dwell = self.delta_k / self.gx.amplitude / self.oversampling - adc_dwell = round_to_raster(adc_dwell, self.system.adc_raster_time, 'ceil') - adc_samples = int(np.ceil(self.readout_time / adc_dwell / 4) * 4) + adc_dwell = delta_kx / self.gx.amplitude + adc_dwell = round_to_raster(adc_dwell, self.system.adc_raster_time) + + adc_samples = int(round(readout_time / adc_dwell / 4) * 4) - adc_time_to_center = adc_dwell * ((adc_samples - 1) / 2 + 0.5) - adc_delay = self.gx.rise_time + self.gx.flat_time / 2 - adc_time_to_center - adc_delay = round_to_raster(adc_delay, self.system.adc_raster_time) + adc_time_to_center = adc_dwell * (adc_samples / 2 + 0.5) + adc_delay = self.gx.rise_time + self.gx.flat_time / 2 - adc_time_to_center + adc_dwell / 2 + # the adc delay has to be rounded to the rf raster time (not the adc raster time) + adc_delay = round_to_raster(adc_delay, self.system.adc_raster_time) # TODO: CHANGE BACK TO RF_RASTER_TIME self.adc = pp.make_adc( num_samples=adc_samples, dwell=adc_dwell, - freq_offset=self.adc_freq_offset, + freq_offset=adc_freq_offset, delay=adc_delay, system=self.system, ) # Create and align pre-phaser gradients considering partial fourier factor - # determine the number of pe lines after (and including) k-space center, which is independent of partial fourier + # determine the number of "PE" lines after (and including) k-space center (independent of partial fourier) self.n_phase_enc_post_center = int(np.ceil(self.n_phase_encoding / 2 + 1)) # find the closest number of lines to the desired partial fourier factor - valid_n_phase_total = int(np.ceil(self.partial_fourier_factor * self.n_phase_encoding)) + valid_n_phase_total = int(np.ceil(partial_fourier_factor * self.n_phase_encoding)) # ensure that at least one line before the center is acquired self.n_phase_enc_pre_center = max(1, valid_n_phase_total - self.n_phase_enc_post_center) + # update the total number of "PE" lines self.n_phase_enc_total = self.n_phase_enc_pre_center + self.n_phase_enc_post_center # recalculate the actual partial fourier factor actual_pf_factor = self.n_phase_enc_total / self.n_phase_encoding - if actual_pf_factor != self.partial_fourier_factor: + if actual_pf_factor != partial_fourier_factor: warnings.warn( - f'Desired partial Fourier factor {self.partial_fourier_factor} adjusted to {actual_pf_factor}.', + f'Desired partial Fourier factor {partial_fourier_factor} adjusted to {actual_pf_factor}.', UserWarning, stacklevel=2, ) self.partial_fourier_factor = actual_pf_factor - self.gx_pre = pp.make_trapezoid(channel='x', system=self.system, area=-self.gx.area / 2 - self.delta_k / 2) + # Create pre-phaser gradients + self.gx_pre = pp.make_trapezoid( + channel='x', + system=self.system, + area=-self.gx.area / 2 - delta_kx / 2, + ) self.gy_pre = pp.make_trapezoid( - channel='y', system=self.system, area=self.n_phase_enc_pre_center * self.delta_k + channel='y', + system=self.system, + area=self.n_phase_enc_pre_center * delta_ky, ) self.gx_pre, self.gy_pre = pp.align(right=[self.gx_pre, self.gy_pre]) + # Create and align "phase encoding" gradients if self.readout_type == 'flyback': - # Create flyback gradient in case of flyback readout self.gx_flyback = pp.make_trapezoid(channel='x', system=self.system, area=-self.gx.area) elif self.readout_type == 'symmetric': # Split and align blip gradient in case of symmetric readout @@ -540,18 +534,75 @@ def __init__( self.gy_blipup.waveform *= 0 self.gy_blipdown.waveform *= 0 self.gy_blipdownup.waveform *= 0 + elif self.readout_type == 'flyback': + self.gy_blip.waveform *= 0 # Create spoiler gradient if spoiling is enabled if self.spoiling_enable: - self.gz_spoil = pp.make_trapezoid(channel='z', system=self.system, area=4 * self.delta_k * n_readout) + self.gz_spoil = pp.make_trapezoid(channel='z', system=self.system, area=4 * delta_ky * n_readout) + + @property + def time_to_center(self) -> float: + """Return time from beginning of readout to center of k-space (needed for TE calculations).""" + # i) add time for pre phaser + time_to_center = pp.calc_duration(self.gx_pre, self.gy_pre) + # ii) add time for completed k-space lines before central (ky = 0) line + if self.readout_type == 'flyback': + time_to_center += self.n_phase_enc_pre_center * ( + pp.calc_duration(self.gx) + pp.calc_duration(self.gx_flyback, self.gy_blip) + ) + elif self.readout_type == 'symmetric': + time_to_center += self.n_phase_enc_pre_center * pp.calc_duration(self.gx) + # iii) add time before start of ADC of central k-space line + if self.ramp_sampling: + time_to_center += self.adc.delay + else: + time_to_center += self.gx.rise_time + # iv) add time from start of ADC (for ky = 0) to timepoint when kx = 0 as well + time_to_center += self.adc.dwell * (self.adc.num_samples / 2 + 0.5) - def add_to_seq(self, seq: pp.Sequence): + return float(time_to_center) + + @property + def time_to_center_without_prephaser(self) -> float: + """Return time from after pre-phasers to center of k-space (needed for TE calculations).""" + return self.time_to_center - pp.calc_duration(self.gx_pre, self.gy_pre) + + @property + def total_duration(self) -> float: + """Return total duration of readout including pre-phaser and optional spoiler.""" + total_duration = pp.calc_duration(self.gx_pre, self.gy_pre) + if self.readout_type == 'flyback': + total_duration += self.n_phase_enc_total * ( + pp.calc_duration(self.gx) + pp.calc_duration(self.gx_flyback, self.gy_blip) + ) + total_duration -= pp.calc_duration(self.gx_flyback, self.gy_blip) + elif self.readout_type == 'symmetric': + total_duration += self.n_phase_enc_total * pp.calc_duration(self.gx) + # add time for spoiler if enabled + if self.spoiling_enable: + total_duration += pp.calc_duration(self.gz_spoil) + + return float(total_duration) + + @property + def total_duration_without_prephaser(self) -> float: + """Return total duration of readout excluding pre-phasers but including optional spoiler.""" + return self.total_duration - pp.calc_duration(self.gx_pre, self.gy_pre) + + def add_to_seq( + self, + seq: pp.Sequence, + add_prephaser: bool = True, + mrd_dataset: ismrmrd.Dataset | None = None, + ) -> tuple[pp.Sequence, ismrmrd.Dataset | None]: """Add EPI readout blocks to the sequence.""" # (Re)set phase encoding (LIN) label lin_label = pp.make_label(label='LIN', type='SET', value=0) - # Add pre-phaser gradients - seq.add_block(self.gx_pre, self.gy_pre) + # Add pre-phaser gradients if enabled + if add_prephaser: + seq.add_block(self.gx_pre, self.gy_pre) for pe_idx in range(self.n_phase_enc_total): rev_label = pp.make_label(type='SET', label='REV', value=self.gx.amplitude < 0) @@ -567,6 +618,7 @@ def add_to_seq(self, seq: pp.Sequence): # Add readout block and reverse polarity of readout gradient seq.add_block(self.gx, gy_blip, self.adc, lin_label, rev_label) self.gx.amplitude = -self.gx.amplitude + elif self.readout_type == 'flyback': seq.add_block(self.gx, self.adc, lin_label, rev_label) if pe_idx != self.n_phase_enc_total - 1: @@ -577,13 +629,17 @@ def add_to_seq(self, seq: pp.Sequence): if self.spoiling_enable: seq.add_block(self.gz_spoil) - return seq + return seq, mrd_dataset def plot_sequence(self): """Plot the sequence.""" seq = pp.Sequence(self.system) - seq = self.add_to_seq(seq) - seq.plot(grad_disp='mT/m') + seq, _ = self.add_to_seq(seq) + _, axs1, _, axs2 = seq.plot(grad_disp='mT/m', plot_now=False) + # add vertical line at time to center to all plots in both figures + for ax in axs1 + axs2: + ax.axvline(x=self.time_to_center, color='r', linestyle='--') + plt.show() def plot_trajectory(self): """Plot k-space trajectory.""" @@ -598,7 +654,7 @@ def plot_trajectory(self): system=self.system, ) ) - seq = self.add_to_seq(seq) + seq, _ = self.add_to_seq(seq) # calculate k-space trajectory k_traj_adc, k_traj, _, _, _ = seq.calculate_kspace() From ec7e3dc4b094a554a9c8aa387f327738d0f116c4 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Mon, 22 Dec 2025 11:27:56 +0100 Subject: [PATCH 04/32] add segmentation (SEG) label for reverse acquisitions --- src/mrseq/utils/trajectory.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py index 63c33ab..1ba3959 100644 --- a/src/mrseq/utils/trajectory.py +++ b/src/mrseq/utils/trajectory.py @@ -469,7 +469,7 @@ def __init__( adc_time_to_center = adc_dwell * (adc_samples / 2 + 0.5) adc_delay = self.gx.rise_time + self.gx.flat_time / 2 - adc_time_to_center + adc_dwell / 2 # the adc delay has to be rounded to the rf raster time (not the adc raster time) - adc_delay = round_to_raster(adc_delay, self.system.adc_raster_time) # TODO: CHANGE BACK TO RF_RASTER_TIME + adc_delay = round_to_raster(adc_delay, self.system.rf_raster_time) self.adc = pp.make_adc( num_samples=adc_samples, @@ -606,6 +606,7 @@ def add_to_seq( for pe_idx in range(self.n_phase_enc_total): rev_label = pp.make_label(type='SET', label='REV', value=self.gx.amplitude < 0) + seg_label = pp.make_label(type='SET', label='SEG', value=self.gx.amplitude < 0) if self.readout_type == 'symmetric': # Select blip gradient based on phase encoding index if pe_idx == 0: @@ -616,7 +617,7 @@ def add_to_seq( gy_blip = self.gy_blipdownup # Add readout block and reverse polarity of readout gradient - seq.add_block(self.gx, gy_blip, self.adc, lin_label, rev_label) + seq.add_block(self.gx, gy_blip, self.adc, lin_label, rev_label, seg_label) self.gx.amplitude = -self.gx.amplitude elif self.readout_type == 'flyback': From 6ae5063c6170561966517b5b0d8d60798288d4b0 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 6 Jan 2026 11:31:45 +0100 Subject: [PATCH 05/32] update epi2d script (deprecated) --- src/mrseq/scripts/epi2d.py | 288 +++++++++++++++++++++++++++++-------- 1 file changed, 229 insertions(+), 59 deletions(-) diff --git a/src/mrseq/scripts/epi2d.py b/src/mrseq/scripts/epi2d.py index 65b1a98..6edae26 100644 --- a/src/mrseq/scripts/epi2d.py +++ b/src/mrseq/scripts/epi2d.py @@ -1,9 +1,11 @@ """2D Echo Planar Imaging (EPI) sequence.""" +from math import floor from pathlib import Path from typing import Literal import ismrmrd +import matplotlib.pyplot as plt import numpy as np import pypulseq as pp @@ -31,11 +33,14 @@ def epi2d_kernel( rf_bwt: float, rf_apodization: float, readout_type: Literal['symmetric', 'flyback'], + echo_type: Literal['FID', 'SE'], oversampling: Literal[1, 2, 4], ramp_sampling: bool, partial_fourier_factor: float, pe_enable: bool, spoiling_enable: bool, + add_noise_acq: bool, + add_navigator_acq: bool, mrd_header_file: str | Path | None, ) -> tuple[pp.Sequence, float, float]: """Generate a 2D Echo Planar Imaging (EPI) sequence. @@ -70,6 +75,8 @@ def epi2d_kernel( Apodization factor of rf excitation pulse readout_type Readout type ('symmetric' or 'flyback'). + echo_type + Echo type ('FID' or 'SE'). oversampling ADC oversampling factor. ramp_sampling @@ -124,39 +131,131 @@ def epi2d_kernel( use='excitation', ) + # create refocussing pulse and gradients if echo type is 'SE (spin echo) + if echo_type == 'SE': + rf180, gz180, _ = pp.make_sinc_pulse( + flip_angle=np.pi, + system=system, + duration=rf_duration * 2, + slice_thickness=slice_thickness, + apodization=0.5, + time_bw_product=4, + phase_offset=np.pi / 2, + use='refocusing', + return_gz=True, + delay=system.rf_dead_time, + ) + _, gzr1_t, gzr1_a = pp.make_extended_trapezoid_area( + channel='z', + grad_start=0, + grad_end=gz180.amplitude, + area=1.5 * gz.area, + system=system, + ) + _, gzr2_t, gzr2_a = pp.make_extended_trapezoid_area( + channel='z', + grad_start=gz180.amplitude, + grad_end=0, + area=-gzr.area + 1.5 * gz.area, + system=system, + ) + + # create combined gradient including pre/post spoiler + gz180n = pp.make_extended_trapezoid( + channel='z', + system=system, + times=np.array([*gzr1_t, *gzr2_t + gzr1_t[3] + gz180.flat_time]), + amplitudes=np.array([*gzr1_a, *gzr2_a]), + ) + + # update rf delay of refocussing pulse to ensure it's centered in plateau of combined gradient + rf180.delay = gzr1_t[-1] + # calculate minimum echo time gzr_prephaser_dur = pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gx_pre) - min_te = rf.shape_dur / 2 # time from center to end of RF pulse - min_te += max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time - min_te += gzr_prephaser_dur # for minimum TE, gzr and pre-phasers are played out simultaneously - min_te += epi2d.time_to_center_without_prephaser - - # calculate echo time delay (te_delay) - if te is None: - te_delay = 0.0 - else: - te_delay = round_to_raster(te - min_te, system.block_duration_raster) - if te_delay < 0: - raise ValueError(f'TE must be larger than {min_te * 1000:.3f} ms. Current value is {te * 1000:.3f} ms.') - # calculate minimum repetition time depending on chosen echo time - min_tr = pp.calc_duration(rf, gz) - min_tr += gzr_prephaser_dur # if TE is chosen minimal, gzr and pre-phasers are played out simultaneously - min_tr += epi2d.total_duration_without_prephaser + # calculate echo time delay(s) + if echo_type == 'FID': + min_te = rf.shape_dur / 2 # time from center to end of RF pulse + min_te += max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time + if add_navigator_acq: + min_te += pp.calc_duration(gzr, epi2d.gx) + min_te += 3 * pp.calc_duration(epi2d.gx) + min_te += pp.calc_duration(epi2d.gy_pre) + else: + min_te += gzr_prephaser_dur # for minimum TE, gzr and pre-phasers are played out simultaneously + min_te += epi2d.time_to_center_without_prephaser + + if te is None: + te_delay = 0.0 + else: + te_delay = round_to_raster(te - min_te, system.block_duration_raster) + if te_delay < 0: + raise ValueError(f'TE must be larger than {min_te * 1000:.3f} ms. Current value is {te * 1000:.3f} ms.') + elif echo_type == 'SE': + t_mid_ref_to_center = rf180.shape_dur / 2 + t_mid_ref_to_center += gzr2_t[-1] + if add_navigator_acq: + t_mid_ref_to_center += pp.calc_duration(epi2d.gx) + t_mid_ref_to_center += 3 * pp.calc_duration(epi2d.gx) + t_mid_ref_to_center += pp.calc_duration(epi2d.gy_pre) + else: + t_mid_ref_to_center += gzr_prephaser_dur + t_mid_ref_to_center += epi2d.time_to_center_without_prephaser + + min_te = 2 * t_mid_ref_to_center + + te_delay2 = 0.0 + te_delay1 = t_mid_ref_to_center + te_delay1 -= rf.shape_dur / 2 + te_delay1 -= max(rf.ringdown_time, gz.fall_time) + te_delay1 -= gzr1_t[-1] + te_delay1 -= rf180.shape_dur / 2 + te_delay1 = round_to_raster(te_delay1, system.block_duration_raster) + + if te is not None: + if te > min_te: + additional_te_delay_half = round_to_raster((te - min_te) / 2, system.block_duration_raster) + te_delay1 += additional_te_delay_half + te_delay2 += additional_te_delay_half + else: + raise ValueError( + f'Desired TE ({te * 1000:.3f} ms) is smaller than minimum TE ({min_te * 1000:.3f} ms).' + ) # calculate repetition time delay (tr_delay) for current TE settings - current_min_tr = min_tr + te_delay + if echo_type == 'FID': + min_tr = pp.calc_duration(rf, gz) + if add_navigator_acq: + min_tr += pp.calc_duration(gzr, epi2d.gx) + min_tr += 3 * pp.calc_duration(epi2d.gx) + min_tr += pp.calc_duration(epi2d.gy_pre) + else: + min_tr += gzr_prephaser_dur + min_tr += epi2d.total_duration_without_prephaser + min_tr += te_delay + else: + min_tr = pp.calc_duration(rf, gz) + min_tr += te_delay1 + min_tr += pp.calc_duration(rf180, gz180n) + min_tr += te_delay2 + if add_navigator_acq: + min_tr += pp.calc_duration(gzr, epi2d.gx) + min_tr += 3 * pp.calc_duration(epi2d.gx) + min_tr += pp.calc_duration(epi2d.gy_pre) + else: + min_tr += gzr_prephaser_dur + min_tr += epi2d.total_duration_without_prephaser + if tr is None: tr_delay = 0.0 else: - tr_delay = round_to_raster(tr - current_min_tr, system.block_duration_raster) + tr_delay = round_to_raster(tr - min_tr, system.block_duration_raster) if tr_delay < 0: - raise ValueError( - f'TR must be larger than {current_min_tr * 1000:.2f} ms. Current value is {tr * 1000:.3f} ms.' - ) + raise ValueError(f'TR must be larger than {min_tr * 1000:.2f} ms. Current value is {tr * 1000:.3f} ms.') - print(f'\nCurrent echo time = {(min_te + te_delay) * 1000:.3f} ms') - print(f'Current repetition time = {(current_min_tr + tr_delay) * 1000:.3f} ms') + # print(f'\nCurrent echo time = {(min_te + te_delay) * 1000:.3f} ms') + # print(f'Current repetition time = {(min_tr + tr_delay) * 1000:.3f} ms') # create header if mrd_header_file: @@ -176,22 +275,25 @@ def epi2d_kernel( prot = ismrmrd.Dataset(mrd_header_file, 'w') prot.write_xml_header(hdr.toXML('utf-8')) - # obtain noise samples - seq.add_block( - pp.make_label(label='LIN', type='SET', value=0), - pp.make_label(label='SLC', type='SET', value=0), - pp.make_label(label='NOISE', type='SET', value=True), - ) - seq.add_block( - epi2d.adc, pp.make_delay(round_to_raster(pp.calc_duration(epi2d.adc), system.block_duration_raster, 'ceil')) - ) - seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) - seq.add_block(pp.make_delay(system.rf_dead_time)) + # obtain noise samples if selected + if add_noise_acq: + seq.add_block( + pp.make_label(label='LIN', type='SET', value=0), + pp.make_label(label='SLC', type='SET', value=0), + pp.make_label(label='NOISE', type='SET', value=True), + ) + seq.add_block( + epi2d.adc, pp.make_delay(round_to_raster(pp.calc_duration(epi2d.adc), system.block_duration_raster, 'ceil')) + ) + seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) + seq.add_block(pp.make_delay(system.rf_dead_time)) - if mrd_header_file: - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) - prot.append_acquisition(acq) + if mrd_header_file: + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) + prot.append_acquisition(acq) + + t_after_noise = sum(seq.block_durations.values()) for slice_ in range(n_slices): # define label(s) @@ -203,14 +305,60 @@ def epi2d_kernel( # add slice selective excitation pulse and set slice label seq.add_block(rf, gz, slice_label) - # add echo time delay - if te_delay > 0: + if echo_type == 'FID' and te_delay > 0: seq.add_block(pp.make_delay(te_delay)) - - # add slice selection rephaser and readout prephaser gradients - gzr, gx_pre, gy_pre = pp.align(left=[gzr], right=[epi2d.gx_pre, epi2d.gy_pre]) - seq.add_block(gzr, gx_pre, gy_pre) - + elif echo_type == 'SE': + seq.add_block(pp.make_delay(te_delay1)) + t_after_te_delay1 = sum(seq.block_durations.values()) + seq.add_block(rf180, gz180n) + if te_delay2 > 0: + seq.add_block(pp.make_delay(te_delay2)) + + # add navigator scans for ghost correction + if add_navigator_acq: + # reverse the readout gradient in advance for navigator + gx_pre = pp.scale_grad(epi2d.gx_pre, -1) + gx = pp.scale_grad(epi2d.gx, -1) + block_content = [ + gx_pre, + pp.make_label(label='NAV', type='SET', value=1), + pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), + ] + if echo_type == 'FID': + block_content.append(gzr) # gzr already included in gz180n for SE + seq.add_block(*block_content) + + # reverse gx_pre back after addBlock + gx_pre = pp.scale_grad(gx_pre, -1) + for n in range(3): + seq.add_block( + gx, + epi2d.adc, + pp.make_label(label='REV', type='SET', value=gx.amplitude < 0), + pp.make_label(label='SEG', type='SET', value=gx.amplitude < 0), + pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), + ) + gx = pp.scale_grad(gx, -1) + if mrd_header_file: + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) + prot.append_acquisition(acq) + + # add gy_pre and reset labels + seq.add_block( + epi2d.gy_pre, + pp.make_label(label='NAV', type='SET', value=0), + pp.make_label(label='AVG', type='SET', value=0), + ) + else: + if echo_type == 'FID': + gzr, gx_pre, gy_pre = pp.align(left=[gzr], right=[epi2d.gx_pre, epi2d.gy_pre]) + seq.add_block(gzr, gx_pre, gy_pre) + else: + gx_pre, gy_pre = pp.align(right=[epi2d.gx_pre, epi2d.gy_pre]) + seq.add_block(gx_pre, gy_pre) + + t_before_readout = sum(seq.block_durations.values()) # add EPI readout block without pre-phaser gradients seq, prot = epi2d.add_to_seq(seq, add_prephaser=False, mrd_dataset=prot) @@ -273,23 +421,33 @@ def epi2d_kernel( seq.set_definition('TR', tr or min_tr) seq.set_definition('ReadoutOversamplingFactor', oversampling) - return seq, min_te, min_tr + t_exc = t_after_noise + rf.delay + rf.shape_dur / 2 + if echo_type == 'SE': + t_ref = t_after_te_delay1 + rf180.delay + rf180.shape_dur / 2 + else: + t_ref = 0.0 + t_echo = t_before_readout + epi2d.time_to_center_without_prephaser + + return seq, min_te, min_tr, t_exc, t_ref, t_echo def main( system: pp.Opts | None = None, te: float | None = None, tr: float | None = None, - fov: float = 256e-3, - n_readout: int = 64, - n_phase_encoding: int = 64, + fov: float = 200e-3, + n_readout: int = 16, + n_phase_encoding: int = 16, n_slices: int = 1, - slice_thickness: float = 8e-3, + slice_thickness: float = 4e-3, bandwidth: float = 64e3, - readout_type: Literal['symmetric', 'flyback'] = 'flyback', - oversampling: Literal[1, 2, 4] = 1, + readout_type: Literal['symmetric', 'flyback'] = 'symmetric', + echo_type: Literal['FID', 'SE'] = 'FID', + oversampling: Literal[1, 2, 4] = 2, ramp_sampling: bool = False, partial_fourier_factor: float = 1, + add_navigator_acq: bool = True, + add_noise_acq: bool = True, show_plots: bool = True, test_report: bool = True, timing_check: bool = True, @@ -307,9 +465,9 @@ def main( system = sys_defaults # define settings of rf excitation pulse - rf_flip_angle: float = 90 # flip angle of the rf excitation pulse [°] - rf_duration = 2.56e-3 # duration of the rf excitation pulse [s] - rf_bwt = 4 # bandwidth-time product of rf excitation pulse [Hz*s] + rf_flip_angle = 90.0 # flip angle of the rf excitation pulse [°] + rf_duration = 1.28e-3 # duration of the rf excitation pulse [s] + rf_bwt = 4.0 # bandwidth-time product of rf excitation pulse [Hz*s] rf_apodization = 0.5 # apodization factor of rf excitation pulse # define EPI settings @@ -318,10 +476,10 @@ def main( # define sequence filename rs_string = 'rs' if ramp_sampling else 'nors' - pf_string = f'{partial_fourier_factor}'.replace('.', 'p') + pf_string = f'{partial_fourier_factor}pf'.replace('.', 'p') filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny' - filename += f'_{readout_type}_{oversampling}ro_{rs_string}_{pf_string}' + filename += f'_{readout_type}_{echo_type}_{oversampling}ro_{rs_string}_{pf_string}_with_nav' output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) @@ -332,7 +490,7 @@ def main( mrd_file = output_path / Path(filename + '_header.h5') - seq, _min_te, min_tr = epi2d_kernel( + seq, _min_te, min_tr, t_exc, t_ref, t_echo = epi2d_kernel( system=system, te=te, tr=tr, @@ -348,13 +506,20 @@ def main( rf_bwt=rf_bwt, rf_apodization=rf_apodization, readout_type=readout_type, + echo_type=echo_type, ramp_sampling=ramp_sampling, partial_fourier_factor=partial_fourier_factor, pe_enable=enable_phase_encoding, spoiling_enable=enable_gradient_spoiling, + add_noise_acq=add_noise_acq, + add_navigator_acq=add_navigator_acq, mrd_header_file=mrd_file, ) + print(f'Time from exc to ref: {(t_ref - t_exc) * 1000}') + print(f'Time from ref to echo: {(t_echo - t_ref) * 1000}') + print(f'Time from exc to echo: {(t_echo - t_exc) * 1000}') + # check timing of the sequence if timing_check and not test_report: ok, error_report = seq.check_timing() @@ -374,7 +539,12 @@ def main( seq.write(str(output_path / filename), create_signature=True) if show_plots: - seq.plot(time_range=(0, 10 * (tr or min_tr))) + fig1, axs1, fig2, axs2 = seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=False) + for ax in [*axs1, *axs2]: + ax.axvline(t_exc, color='red') + ax.axvline(t_ref, color='red') + ax.axvline(t_echo, color='red') + plt.show() return seq, output_path / filename From 005ec8dbdc74023997228aea9c8f00b435aac749 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 6 Jan 2026 11:37:30 +0100 Subject: [PATCH 06/32] split epi2d script into dedicated epi2d_fid and epi2d_se scripts (WIP) --- src/mrseq/scripts/epi2d_fid.py | 439 ++++++++++++++++++++ src/mrseq/scripts/{epi2d.py => epi2d_se.py} | 317 +++++++------- 2 files changed, 589 insertions(+), 167 deletions(-) create mode 100644 src/mrseq/scripts/epi2d_fid.py rename src/mrseq/scripts/{epi2d.py => epi2d_se.py} (64%) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py new file mode 100644 index 0000000..9038823 --- /dev/null +++ b/src/mrseq/scripts/epi2d_fid.py @@ -0,0 +1,439 @@ +"""2D Echo Planar Imaging (EPI) sequence.""" + +from math import floor +from pathlib import Path +from typing import Literal + +import ismrmrd +import numpy as np +import pypulseq as pp + +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults +from mrseq.utils.ismrmrd import Fov +from mrseq.utils.ismrmrd import Limits +from mrseq.utils.ismrmrd import MatrixSize +from mrseq.utils.ismrmrd import create_header +from mrseq.utils.trajectory import EpiReadout + + +def epi2d_fid_kernel( + system: pp.Opts, + te: float | None, + tr: float | None, + fov: float, + n_readout: int, + n_phase_encoding: int, + bandwidth: float, + slice_thickness: float, + n_slices: int, + rf_duration: float, + rf_flip_angle: float, + rf_bwt: float, + rf_apodization: float, + readout_type: Literal['symmetric', 'flyback'], + oversampling: Literal[1, 2, 4], + ramp_sampling: bool, + partial_fourier_factor: float, + add_spoiler: bool, + add_noise_acq: bool, + add_navigator_acq: bool, + mrd_header_file: str | Path | None, +) -> tuple[pp.Sequence, float, float]: + """Generate a 2D Echo Planar Imaging (EPI) sequence. + + Parameters + ---------- + system + PyPulseq system limits object. + te + Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. + tr + Desired repetition time (TR) (in seconds). + fov + Field of view in x and y direction (in meters). + n_readout + Number of frequency encoding steps. + n_phase_encoding + Number of phase encoding steps. + bandwidth + Total receiver bandwidth (in Hz). + slice_thickness + Slice thickness of the 2D slice (in meters). + n_slices + Number of slices. + rf_duration + Duration of the rf excitation pulse (in seconds) + rf_flip_angle + Flip angle of rf excitation pulse (in degrees) + rf_bwt + Bandwidth-time product of rf excitation pulse (in Hz * seconds) + rf_apodization + Apodization factor of rf excitation pulse. + readout_type + Readout type ('symmetric' or 'flyback'). + oversampling + Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. + ramp_sampling + If True, ADC is active during gradient ramps for optimized timing. + partial_fourier_factor + Desired partial Fourier factor in "phase encoding" direction. Must be larger than 0.5 and smaller or equal to 1. + The actual partial Fourier factor might slightly deviate from the desired value. + add_spoiler + If True, a spoiler gradients will be added to the sequence after the EPI readout. + add_noise_acq + If True, noise acquisitions will be added at the beginning of the sequence. + add_navigator_acq + If True, 3 navigator acquisitions will be added to the sequence to allow for ghost corrections. + The navigator acquisitions are added between the rf excitation pulse and the EPI readout. + Be aware that navigator acquisitions will increase the minimum echo and repetition times. + mrd_header_file + Filename of the ISMRMRD header file to be created. If None, no header file is created. + + Returns + ------- + seq + PyPulseq Sequence object + min_te + Shortest possible echo time. + min_tr + Shortest possible repetition time. + + """ + # create PyPulseq Sequence object and set system limits + seq = pp.Sequence(system=system) + + # create EpiReadout object + epi2d = EpiReadout( + system=system, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=bandwidth, + oversampling=oversampling, + readout_type=readout_type, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, + pe_enable=True, + spoiling_enable=add_spoiler, + ) + + # create slice selective excitation pulse and gradients + rf, gz, gzr = pp.make_sinc_pulse( + flip_angle=np.deg2rad(rf_flip_angle), + duration=rf_duration, + slice_thickness=slice_thickness, + apodization=rf_apodization, + time_bw_product=rf_bwt, + delay=system.rf_dead_time, + system=system, + return_gz=True, + use='excitation', + ) + + # calculate echo time delay + min_te = rf.shape_dur / 2 # time from center to end of RF pulse + min_te += max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time + if add_navigator_acq: + min_te += pp.calc_duration(gzr, epi2d.gx_pre) + min_te += 3 * pp.calc_duration(epi2d.gx) + min_te += pp.calc_duration(epi2d.gy_pre) + else: + min_te += pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gx_pre) + min_te += epi2d.time_to_center_without_prephaser + + if te is None: + te_delay = 0.0 + else: + te_delay = round_to_raster(te - min_te, system.block_duration_raster) + if te_delay < 0: + raise ValueError(f'TE must be larger than {min_te * 1000:.3f} ms. Current value is {te * 1000:.3f} ms.') + + # calculate repetition time delay (tr_delay) for current TE settings + min_tr = pp.calc_duration(rf, gz) + if add_navigator_acq: + min_tr += pp.calc_duration(gzr, epi2d.gx_pre) + min_tr += 3 * pp.calc_duration(epi2d.gx) + min_tr += pp.calc_duration(epi2d.gy_pre) + else: + min_tr += pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gx_pre) + min_tr += epi2d.total_duration_without_prephaser + min_tr += te_delay + + if tr is None: + tr_delay = 0.0 + else: + tr_delay = round_to_raster(tr - min_tr, system.block_duration_raster) + if tr_delay < 0: + raise ValueError(f'TR must be larger than {min_tr * 1000:.2f} ms. Current value is {tr * 1000:.3f} ms.') + + # create header + if mrd_header_file: + hdr = create_header( + traj_type='other', + encoding_fov=Fov(x=fov * oversampling, y=fov, z=slice_thickness), + recon_fov=Fov(x=fov, y=fov, z=slice_thickness), + encoding_matrix=MatrixSize(n_x=n_readout, n_y=epi2d.n_phase_encoding, n_z=1), + recon_matrix=MatrixSize(n_x=n_readout, n_y=n_phase_encoding, n_z=1), + dwell_time=epi2d.adc.dwell, + slice_limits=Limits(min=0, max=n_slices, center=n_slices // 2), + k1_limits=Limits(min=0, max=epi2d.n_phase_enc_total, center=epi2d.n_phase_enc_pre_center + 1), + k2_limits=Limits(min=0, max=1, center=0), + ) + + # write header to file + prot = ismrmrd.Dataset(mrd_header_file, 'w') + prot.write_xml_header(hdr.toXML('utf-8')) + + # obtain noise samples if selected + if add_noise_acq: + seq.add_block( + pp.make_label(label='LIN', type='SET', value=0), + pp.make_label(label='SLC', type='SET', value=0), + pp.make_label(label='NOISE', type='SET', value=True), + ) + seq.add_block( + epi2d.adc, pp.make_delay(round_to_raster(pp.calc_duration(epi2d.adc), system.block_duration_raster, 'ceil')) + ) + seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) + seq.add_block(pp.make_delay(system.rf_dead_time)) + + if mrd_header_file: + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) + prot.append_acquisition(acq) + + for slice_ in range(n_slices): + # define label(s) + slice_label = pp.make_label(label='SLC', type='SET', value=slice_) + + # set frequency offset for current slice + rf.freq_offset = gz.amplitude * slice_thickness * (slice_ - (n_slices - 1) / 2) + # rf.phase_offset = - 2 * np.pi * rf.freq_offset * pp.calc_rf_center(rf) + + # add slice-selective excitation pulse and set slice label + seq.add_block(rf, gz, slice_label) + + # add navigator scans for ghost correction + if add_navigator_acq: + # reverse the readout gradient and pre-winder in advance for navigator + gx_pre = pp.scale_grad(epi2d.gx_pre, -1) + gx = pp.scale_grad(epi2d.gx, -1) + # add slice selection rewinder and readout pre-winder in x direction (gy_pre will be added after navigators) + gzr, gx_pre = pp.align(left=[gzr], right=[gx_pre]) + seq.add_block( + gzr, + gx_pre, + pp.make_label(label='NAV', type='SET', value=1), + pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), + ) + gx_pre = pp.scale_grad(gx_pre, -1) + for n in range(3): + seq.add_block( + gx, + epi2d.adc, + pp.make_label(label='REV', type='SET', value=gx.amplitude < 0), + pp.make_label(label='SEG', type='SET', value=gx.amplitude < 0), + pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), + ) + gx = pp.scale_grad(gx, -1) + # add navigator acquisitions to ISMRMRD file + if mrd_header_file: + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) + prot.append_acquisition(acq) + + # add echo time delay + seq.add_block(pp.make_delay(te_delay)) + + # add gy_pre and reset labels + seq.add_block( + epi2d.gy_pre, + pp.make_label(label='NAV', type='SET', value=0), + pp.make_label(label='AVG', type='SET', value=0), + ) + else: + # add echo time delay + seq.add_block(pp.make_delay(te_delay)) + + # align and add slice-selection rewinder and readout pre-winder gradients + gzr, gx_pre, gy_pre = pp.align(left=[gzr], right=[epi2d.gx_pre, epi2d.gy_pre]) + seq.add_block(gzr, gx_pre, gy_pre) + + # add EPI readout block without pre-phaser gradients + seq, prot = epi2d.add_to_seq(seq, add_prephaser=False, mrd_dataset=prot) + + # add repetition time delay + if tr_delay > 0: + seq.add_block(pp.make_delay(tr_delay)) + + if mrd_header_file: + # calculate k-space trajectory from Sequence to add this info to mrd file + k_traj_adc, _, _, _, _ = seq.calculate_kspace() + samples_per_acq = epi2d.adc.num_samples + number_of_total_acq = k_traj_adc.shape[-1] // samples_per_acq + number_of_epi_acq = epi2d.n_phase_enc_total + if 'data' in prot._dataset: + number_of_noise_acq = prot.number_of_acquisitions() + else: + number_of_noise_acq = 0 + + if not number_of_epi_acq + number_of_noise_acq == number_of_total_acq: + raise (ValueError('Number of calculated acquisitions does not match expected number.')) + + k_traj_adc_readout = k_traj_adc[:, number_of_noise_acq * samples_per_acq :] + + # create mrd acquisition and add trajectory info for each EPI readout + for n in range(number_of_epi_acq): + start = n * samples_per_acq + end = start + samples_per_acq + + traj = np.zeros((samples_per_acq, 2), dtype=np.float32) + traj[:, 0] = np.round(k_traj_adc_readout[0, start:end] * fov, 3) + traj[:, 1] = np.round(k_traj_adc_readout[1, start:end] * fov, 3) + + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=samples_per_acq) + acq.traj[:] = traj + + prot.append_acquisition(acq) + + # close ISMRMRD file + if mrd_header_file: + prot.close() + + # set gridding definitions extracted from EpiReadout + if ramp_sampling: + gridding_params = [ + epi2d.gx.rise_time, + epi2d.gx.flat_time, + epi2d.gx.fall_time, + epi2d.adc.delay - epi2d.gx.delay, + epi2d.adc.num_samples * epi2d.adc.dwell, + ] + + seq.set_definition(key='TargetGriddedSamples', value=n_readout * oversampling) + seq.set_definition(key='TrapezoidGriddingParameters', value=gridding_params) + + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov, fov, slice_thickness * n_slices]) + seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, n_slices)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr or min_tr) + seq.set_definition('ReadoutOversamplingFactor', oversampling) + + return seq, min_te, min_tr + + +def main( + system: pp.Opts | None = None, + te: float | None = None, + tr: float | None = None, + fov: float = 200e-3, + n_readout: int = 64, + n_phase_encoding: int = 64, + n_slices: int = 1, + slice_thickness: float = 8e-3, + bandwidth: float = 100e3, + readout_type: Literal['symmetric', 'flyback'] = 'symmetric', + oversampling: Literal[1, 2, 4] = 2, + ramp_sampling: bool = True, + partial_fourier_factor: float = 0.7, + add_navigator_acq: bool = True, + add_noise_acq: bool = True, + show_plots: bool = True, + test_report: bool = True, + timing_check: bool = True, +) -> tuple[pp.Sequence, Path]: + """Generate an Echo-Planar Imaging (EPI) sequence. + + Returns + ------- + seq + Sequence object of radial FLASH sequence. + file_path + Path to the sequence file. + """ + if system is None: + system = sys_defaults + + # define settings of rf excitation pulse + rf_flip_angle = 90.0 # flip angle of the rf excitation pulse [°] + rf_duration = 1.28e-3 # duration of the rf excitation pulse [s] + rf_bwt = 4.0 # bandwidth-time product of rf excitation pulse [Hz*s] + rf_apodization = 0.5 # apodization factor of rf excitation pulse + + # define spoiling settings + enable_gradient_spoiling = True + + # define sequence filename + rs_string = 'rs' if ramp_sampling else 'nors' + pf_string = f'{partial_fourier_factor}pf'.replace('.', 'p') + readout_string = 'sym' if readout_type == 'symmetric' else 'flyb' + noise_string = 'withnoise' if add_noise_acq else 'nonoise' + nav_string = 'withnav' if add_navigator_acq else 'nonav' + + filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}px' + filename += f'_{readout_string}_fid_{oversampling}ro_{rs_string}_{pf_string}' + filename += f'_{noise_string}_{nav_string}' + + output_path = Path.cwd() / 'output' + output_path.mkdir(parents=True, exist_ok=True) + + # delete existing header file + if (output_path / Path(filename + '_header.h5')).exists(): + (output_path / Path(filename + '_header.h5')).unlink() + + mrd_file = output_path / Path(filename + '_header.h5') + + seq, _min_te, _min_tr = epi2d_fid_kernel( + system=system, + te=te, + tr=tr, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=bandwidth, + oversampling=oversampling, + slice_thickness=slice_thickness, + n_slices=n_slices, + rf_duration=rf_duration, + rf_flip_angle=rf_flip_angle, + rf_bwt=rf_bwt, + rf_apodization=rf_apodization, + readout_type=readout_type, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, + add_spoiler=enable_gradient_spoiling, + add_noise_acq=add_noise_acq, + add_navigator_acq=add_navigator_acq, + mrd_header_file=mrd_file, + ) + + # check timing of the sequence + if timing_check and not test_report: + ok, error_report = seq.check_timing() + if ok: + print('\nTiming check passed successfully') + else: + print('\nTiming check failed! Error listing follows\n') + print(error_report) + + # show advanced rest report + if test_report: + print('\nCreating advanced test report...') + print(seq.test_report()) + + # save seq-file to disk + print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") + seq.write(str(output_path / filename), create_signature=True) + + if show_plots: + seq.plot(time_range=(0, tr or _min_tr), plot_now=True) + + return seq, output_path / filename + + +if __name__ == '__main__': + main() diff --git a/src/mrseq/scripts/epi2d.py b/src/mrseq/scripts/epi2d_se.py similarity index 64% rename from src/mrseq/scripts/epi2d.py rename to src/mrseq/scripts/epi2d_se.py index 6edae26..c2db65a 100644 --- a/src/mrseq/scripts/epi2d.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -18,7 +18,7 @@ from mrseq.utils.trajectory import EpiReadout -def epi2d_kernel( +def epi2d_se_kernel( system: pp.Opts, te: float | None, tr: float | None, @@ -33,12 +33,10 @@ def epi2d_kernel( rf_bwt: float, rf_apodization: float, readout_type: Literal['symmetric', 'flyback'], - echo_type: Literal['FID', 'SE'], oversampling: Literal[1, 2, 4], ramp_sampling: bool, partial_fourier_factor: float, - pe_enable: bool, - spoiling_enable: bool, + add_spoiler: bool, add_noise_acq: bool, add_navigator_acq: bool, mrd_header_file: str | Path | None, @@ -83,10 +81,12 @@ def epi2d_kernel( If True, ADC is active during gradient ramps (optimized timing). partial_fourier_factor Partial Fourier factor (0.5 to 1.0). - pe_enable - Enable phase encoding (useful for calibration scans if False). - spoiling_enable + add_spoiler Enable spoiling gradients. + add_noise_acq + Enable noise acquisition. + add_navigator_acq + Enable navigator acquisition. mrd_header_file Filename of the ISMRMRD header file to be created. If None, no header file is created. @@ -114,8 +114,8 @@ def epi2d_kernel( readout_type=readout_type, ramp_sampling=ramp_sampling, partial_fourier_factor=partial_fourier_factor, - pe_enable=pe_enable, - spoiling_enable=spoiling_enable, + pe_enable=True, + spoiling_enable=add_spoiler, ) # create slice selective excitation pulse and gradients @@ -132,120 +132,93 @@ def epi2d_kernel( ) # create refocussing pulse and gradients if echo type is 'SE (spin echo) - if echo_type == 'SE': - rf180, gz180, _ = pp.make_sinc_pulse( - flip_angle=np.pi, - system=system, - duration=rf_duration * 2, - slice_thickness=slice_thickness, - apodization=0.5, - time_bw_product=4, - phase_offset=np.pi / 2, - use='refocusing', - return_gz=True, - delay=system.rf_dead_time, - ) - _, gzr1_t, gzr1_a = pp.make_extended_trapezoid_area( - channel='z', - grad_start=0, - grad_end=gz180.amplitude, - area=1.5 * gz.area, - system=system, - ) - _, gzr2_t, gzr2_a = pp.make_extended_trapezoid_area( - channel='z', - grad_start=gz180.amplitude, - grad_end=0, - area=-gzr.area + 1.5 * gz.area, - system=system, - ) + rf180, gz180, _ = pp.make_sinc_pulse( + flip_angle=np.pi, + system=system, + duration=rf_duration * 2, + slice_thickness=slice_thickness, + apodization=0.5, + time_bw_product=4, + phase_offset=np.pi / 2, + use='refocusing', + return_gz=True, + delay=system.rf_dead_time, + ) - # create combined gradient including pre/post spoiler - gz180n = pp.make_extended_trapezoid( - channel='z', - system=system, - times=np.array([*gzr1_t, *gzr2_t + gzr1_t[3] + gz180.flat_time]), - amplitudes=np.array([*gzr1_a, *gzr2_a]), - ) + _spoil_factor = 1.5 + _, gzr1_t, gzr1_a = pp.make_extended_trapezoid_area( + channel='z', + grad_start=0, + grad_end=gz180.amplitude, + area=_spoil_factor * gz.area, + system=system, + ) + _, gzr2_t, gzr2_a = pp.make_extended_trapezoid_area( + channel='z', + grad_start=gz180.amplitude, + grad_end=0, + area=_spoil_factor * gz.area, + system=system, + ) - # update rf delay of refocussing pulse to ensure it's centered in plateau of combined gradient - rf180.delay = gzr1_t[-1] + # create combined gradient including pre/post spoiler + gz180n = pp.make_extended_trapezoid( + channel='z', + system=system, + times=np.array([*gzr1_t, *gzr2_t + gzr1_t[3] + gz180.flat_time]), + amplitudes=np.array([*gzr1_a, *gzr2_a]), + ) - # calculate minimum echo time - gzr_prephaser_dur = pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gx_pre) + # update rf delay of refocussing pulse to ensure it's centered in plateau of combined gradient + rf180.delay = gzr1_t[-1] # calculate echo time delay(s) - if echo_type == 'FID': - min_te = rf.shape_dur / 2 # time from center to end of RF pulse - min_te += max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time - if add_navigator_acq: - min_te += pp.calc_duration(gzr, epi2d.gx) - min_te += 3 * pp.calc_duration(epi2d.gx) - min_te += pp.calc_duration(epi2d.gy_pre) - else: - min_te += gzr_prephaser_dur # for minimum TE, gzr and pre-phasers are played out simultaneously - min_te += epi2d.time_to_center_without_prephaser + t_exc_to_ref = rf.shape_dur / 2 + t_exc_to_ref += max(rf.ringdown_time, gz.fall_time) + if add_navigator_acq: + t_exc_to_ref += pp.calc_duration(gzr, epi2d.gx_pre) + t_exc_to_ref += 3 * pp.calc_duration(epi2d.gx) + t_exc_to_ref += pp.calc_duration(epi2d.gy_pre) + else: + t_exc_to_ref += pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gy_pre) + t_exc_to_ref += gzr1_t[-1] + t_exc_to_ref += rf180.shape_dur / 2 - if te is None: - te_delay = 0.0 - else: - te_delay = round_to_raster(te - min_te, system.block_duration_raster) - if te_delay < 0: - raise ValueError(f'TE must be larger than {min_te * 1000:.3f} ms. Current value is {te * 1000:.3f} ms.') - elif echo_type == 'SE': - t_mid_ref_to_center = rf180.shape_dur / 2 - t_mid_ref_to_center += gzr2_t[-1] - if add_navigator_acq: - t_mid_ref_to_center += pp.calc_duration(epi2d.gx) - t_mid_ref_to_center += 3 * pp.calc_duration(epi2d.gx) - t_mid_ref_to_center += pp.calc_duration(epi2d.gy_pre) + t_ref_to_kcenter = rf180.shape_dur / 2 + t_ref_to_kcenter += gzr2_t[-1] + t_ref_to_kcenter += epi2d.time_to_center_without_prephaser + + # calculate minimum echo time + min_te = 2 * round_to_raster(max(t_exc_to_ref, t_ref_to_kcenter), system.block_duration_raster) + + # calculate echo time delays for minimum echo time + te_delay1 = round_to_raster(min_te / 2 - t_exc_to_ref, system.block_duration_raster) + te_delay2 = round_to_raster(min_te / 2 - t_ref_to_kcenter, system.block_duration_raster) + + if te is not None: + if te > min_te: + # calculate echo time delays for desired echo time + additional_te_delay_half = round_to_raster((te - min_te) / 2, system.block_duration_raster) + te_delay1 += additional_te_delay_half + te_delay2 += additional_te_delay_half else: - t_mid_ref_to_center += gzr_prephaser_dur - t_mid_ref_to_center += epi2d.time_to_center_without_prephaser - - min_te = 2 * t_mid_ref_to_center - - te_delay2 = 0.0 - te_delay1 = t_mid_ref_to_center - te_delay1 -= rf.shape_dur / 2 - te_delay1 -= max(rf.ringdown_time, gz.fall_time) - te_delay1 -= gzr1_t[-1] - te_delay1 -= rf180.shape_dur / 2 - te_delay1 = round_to_raster(te_delay1, system.block_duration_raster) - - if te is not None: - if te > min_te: - additional_te_delay_half = round_to_raster((te - min_te) / 2, system.block_duration_raster) - te_delay1 += additional_te_delay_half - te_delay2 += additional_te_delay_half - else: - raise ValueError( - f'Desired TE ({te * 1000:.3f} ms) is smaller than minimum TE ({min_te * 1000:.3f} ms).' - ) + raise ValueError(f'Desired TE ({te * 1000:.3f} ms) is smaller than minimum TE ({min_te * 1000:.3f} ms).') + + t_exc_to_ref = t_exc_to_ref + te_delay1 # might NOT be on block_raster + t_ref_to_kcenter = t_ref_to_kcenter + te_delay2 # might NOT be on block_raster # calculate repetition time delay (tr_delay) for current TE settings - if echo_type == 'FID': - min_tr = pp.calc_duration(rf, gz) - if add_navigator_acq: - min_tr += pp.calc_duration(gzr, epi2d.gx) - min_tr += 3 * pp.calc_duration(epi2d.gx) - min_tr += pp.calc_duration(epi2d.gy_pre) - else: - min_tr += gzr_prephaser_dur - min_tr += epi2d.total_duration_without_prephaser - min_tr += te_delay + min_tr = pp.calc_duration(rf, gz) + if add_navigator_acq: + min_tr += pp.calc_duration(gzr, epi2d.gx_pre) + min_tr += 3 * pp.calc_duration(epi2d.gx) + min_tr += pp.calc_duration(epi2d.gy_pre) else: - min_tr = pp.calc_duration(rf, gz) - min_tr += te_delay1 - min_tr += pp.calc_duration(rf180, gz180n) - min_tr += te_delay2 - if add_navigator_acq: - min_tr += pp.calc_duration(gzr, epi2d.gx) - min_tr += 3 * pp.calc_duration(epi2d.gx) - min_tr += pp.calc_duration(epi2d.gy_pre) - else: - min_tr += gzr_prephaser_dur - min_tr += epi2d.total_duration_without_prephaser + min_tr += pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gy_pre) + min_tr += te_delay1 + min_tr += pp.calc_duration(rf180, gz180n) + min_tr += te_delay2 + min_tr += epi2d.total_duration_without_prephaser if tr is None: tr_delay = 0.0 @@ -254,8 +227,8 @@ def epi2d_kernel( if tr_delay < 0: raise ValueError(f'TR must be larger than {min_tr * 1000:.2f} ms. Current value is {tr * 1000:.3f} ms.') - # print(f'\nCurrent echo time = {(min_te + te_delay) * 1000:.3f} ms') - # print(f'Current repetition time = {(min_tr + tr_delay) * 1000:.3f} ms') + print(f'\nCurrent echo time = {(t_exc_to_ref + t_ref_to_kcenter) * 1000:.4f} ms') + print(f'Current repetition time = {(min_tr + tr_delay) * 1000:.4f} ms') # create header if mrd_header_file: @@ -301,32 +274,25 @@ def epi2d_kernel( # set frequency offset for current slice rf.freq_offset = gz.amplitude * slice_thickness * (slice_ - (n_slices - 1) / 2) + # rf.phase_offset = - 2 * np.pi * rf.freq_offset * pp.calc_rf_center(rf) # add slice selective excitation pulse and set slice label seq.add_block(rf, gz, slice_label) - if echo_type == 'FID' and te_delay > 0: - seq.add_block(pp.make_delay(te_delay)) - elif echo_type == 'SE': - seq.add_block(pp.make_delay(te_delay1)) - t_after_te_delay1 = sum(seq.block_durations.values()) - seq.add_block(rf180, gz180n) - if te_delay2 > 0: - seq.add_block(pp.make_delay(te_delay2)) - # add navigator scans for ghost correction if add_navigator_acq: - # reverse the readout gradient in advance for navigator - gx_pre = pp.scale_grad(epi2d.gx_pre, -1) - gx = pp.scale_grad(epi2d.gx, -1) - block_content = [ + # reverse the readout gradient and pre-winder in advance for navigator + gx_pre = pp.scale_grad(epi2d.gx_pre, 1) + gx = pp.scale_grad(epi2d.gx, 1) + # add slice selection rewinder and readout pre-winder in x direction (gy_pre will be added after navigators) + gzr, gx_pre = pp.align(left=[gzr], right=[gx_pre]) + seq.add_block( + gzr, gx_pre, pp.make_label(label='NAV', type='SET', value=1), pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), - ] - if echo_type == 'FID': - block_content.append(gzr) # gzr already included in gz180n for SE - seq.add_block(*block_content) + ) + gx_pre = pp.scale_grad(gx_pre, -1) # reverse gx_pre back after addBlock gx_pre = pp.scale_grad(gx_pre, -1) @@ -346,19 +312,22 @@ def epi2d_kernel( # add gy_pre and reset labels seq.add_block( - epi2d.gy_pre, + pp.scale_grad(epi2d.gy_pre, -1), pp.make_label(label='NAV', type='SET', value=0), pp.make_label(label='AVG', type='SET', value=0), ) else: - if echo_type == 'FID': - gzr, gx_pre, gy_pre = pp.align(left=[gzr], right=[epi2d.gx_pre, epi2d.gy_pre]) - seq.add_block(gzr, gx_pre, gy_pre) - else: - gx_pre, gy_pre = pp.align(right=[epi2d.gx_pre, epi2d.gy_pre]) - seq.add_block(gx_pre, gy_pre) - - t_before_readout = sum(seq.block_durations.values()) + seq.add_block(gzr, pp.scale_grad(epi2d.gx_pre, -1), pp.scale_grad(epi2d.gy_pre, -1)) + + if te_delay1 > 0: + seq.add_block(pp.make_delay(te_delay1)) + + # add refocusing pulse + combined gradient + seq.add_block(rf180, gz180n) + + if te_delay2 > 0: + seq.add_block(pp.make_delay(te_delay2)) + # add EPI readout block without pre-phaser gradients seq, prot = epi2d.add_to_seq(seq, add_prephaser=False, mrd_dataset=prot) @@ -417,16 +386,13 @@ def epi2d_kernel( seq.set_definition('FOV', [fov, fov, slice_thickness * n_slices]) seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, n_slices)) seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr or min_tr) + seq.set_definition('TE', float(t_exc_to_ref + t_ref_to_kcenter)) + seq.set_definition('TR', tr or float(min_tr)) seq.set_definition('ReadoutOversamplingFactor', oversampling) t_exc = t_after_noise + rf.delay + rf.shape_dur / 2 - if echo_type == 'SE': - t_ref = t_after_te_delay1 + rf180.delay + rf180.shape_dur / 2 - else: - t_ref = 0.0 - t_echo = t_before_readout + epi2d.time_to_center_without_prephaser + t_ref = t_exc + t_exc_to_ref + t_echo = t_exc + t_exc_to_ref + t_ref_to_kcenter return seq, min_te, min_tr, t_exc, t_ref, t_echo @@ -436,16 +402,15 @@ def main( te: float | None = None, tr: float | None = None, fov: float = 200e-3, - n_readout: int = 16, - n_phase_encoding: int = 16, + n_readout: int = 64, + n_phase_encoding: int = 64, n_slices: int = 1, - slice_thickness: float = 4e-3, - bandwidth: float = 64e3, + slice_thickness: float = 8e-3, + bandwidth: float = 100e3, readout_type: Literal['symmetric', 'flyback'] = 'symmetric', - echo_type: Literal['FID', 'SE'] = 'FID', oversampling: Literal[1, 2, 4] = 2, - ramp_sampling: bool = False, - partial_fourier_factor: float = 1, + ramp_sampling: bool = True, + partial_fourier_factor: float = 0.7, add_navigator_acq: bool = True, add_noise_acq: bool = True, show_plots: bool = True, @@ -470,16 +435,19 @@ def main( rf_bwt = 4.0 # bandwidth-time product of rf excitation pulse [Hz*s] rf_apodization = 0.5 # apodization factor of rf excitation pulse - # define EPI settings - enable_phase_encoding = True + # define spoiling settings enable_gradient_spoiling = True # define sequence filename rs_string = 'rs' if ramp_sampling else 'nors' pf_string = f'{partial_fourier_factor}pf'.replace('.', 'p') + readout_string = 'sym' if readout_type == 'symmetric' else 'flyb' + noise_string = 'withnoise' if add_noise_acq else 'nonoise' + nav_string = 'withnav' if add_navigator_acq else 'nonav' - filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny' - filename += f'_{readout_type}_{echo_type}_{oversampling}ro_{rs_string}_{pf_string}_with_nav' + filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}px' + filename += f'_{readout_string}_se_{oversampling}ro_{rs_string}_{pf_string}' + filename += f'_{noise_string}_{nav_string}' output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) @@ -490,7 +458,7 @@ def main( mrd_file = output_path / Path(filename + '_header.h5') - seq, _min_te, min_tr, t_exc, t_ref, t_echo = epi2d_kernel( + seq, _min_te, min_tr, t_exc, t_ref, t_echo = epi2d_se_kernel( system=system, te=te, tr=tr, @@ -506,19 +474,24 @@ def main( rf_bwt=rf_bwt, rf_apodization=rf_apodization, readout_type=readout_type, - echo_type=echo_type, ramp_sampling=ramp_sampling, partial_fourier_factor=partial_fourier_factor, - pe_enable=enable_phase_encoding, - spoiling_enable=enable_gradient_spoiling, + add_spoiler=enable_gradient_spoiling, add_noise_acq=add_noise_acq, add_navigator_acq=add_navigator_acq, mrd_header_file=mrd_file, ) - print(f'Time from exc to ref: {(t_ref - t_exc) * 1000}') - print(f'Time from ref to echo: {(t_echo - t_ref) * 1000}') - print(f'Time from exc to echo: {(t_echo - t_exc) * 1000}') + # print(f'min_te = {_min_te * 1000:.4f}') + # print(f'min_tr = {min_tr * 1000:.4f}') + # print(f't_exc = {t_exc * 1000:.4f}') + # print(f't_ref = {t_ref * 1000:.4f}') + # print(f't_echo = {t_echo * 1000:.4f}') + # print(f'Echo time (TE) = {(t_echo - t_exc) * 1000:.4f}') + + # print(f'Time from exc to ref: {(t_ref - t_exc) * 1000:.4f}') + # print(f'Time from ref to echo: {(t_echo - t_ref) * 1000:.4f}') + # print(f'added together = {((t_ref - t_exc) + (t_echo - t_ref)) * 1000:.4f}') # check timing of the sequence if timing_check and not test_report: @@ -538,8 +511,18 @@ def main( print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") seq.write(str(output_path / filename), create_signature=True) + # # calculate k-space trajectory + # k_traj_adc, k_traj, _, _, _ = seq.calculate_kspace() + + # # plot trajectory + # fig = plt.figure() + # plt.plot(k_traj[0], k_traj[1], 'b') + # plt.plot(k_traj_adc[0], k_traj_adc[1], 'x', color='red', markersize=4) + # plt.grid() + # plt.show() + if show_plots: - fig1, axs1, fig2, axs2 = seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=False) + _, axs1, _, axs2 = seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=False) for ax in [*axs1, *axs2]: ax.axvline(t_exc, color='red') ax.axvline(t_ref, color='red') From 406e7a848662145c2407e52064ce1d5bd7c55cac Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 3 Feb 2026 09:49:08 +0100 Subject: [PATCH 07/32] minor adjustments to epi2d_se --- src/mrseq/scripts/epi2d_se.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index c2db65a..ccd005b 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -165,7 +165,7 @@ def epi2d_se_kernel( gz180n = pp.make_extended_trapezoid( channel='z', system=system, - times=np.array([*gzr1_t, *gzr2_t + gzr1_t[3] + gz180.flat_time]), + times=np.array([*gzr1_t, *gzr2_t + gzr1_t[-1] + gz180.flat_time]), amplitudes=np.array([*gzr1_a, *gzr2_a]), ) @@ -522,12 +522,7 @@ def main( # plt.show() if show_plots: - _, axs1, _, axs2 = seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=False) - for ax in [*axs1, *axs2]: - ax.axvline(t_exc, color='red') - ax.axvline(t_ref, color='red') - ax.axvline(t_echo, color='red') - plt.show() + seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=False) return seq, output_path / filename From 7459abe6d556395bee9f898ff6d919f679671221 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 24 Feb 2026 08:52:46 +0100 Subject: [PATCH 08/32] cleanup docstrings --- src/mrseq/scripts/epi2d_fid.py | 62 ++++++++++++++--- src/mrseq/scripts/epi2d_se.py | 119 ++++++++++++++++++--------------- 2 files changed, 118 insertions(+), 63 deletions(-) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index 9038823..c476ca1 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -1,4 +1,4 @@ -"""2D Echo Planar Imaging (EPI) sequence.""" +"""2D Echo Planar Imaging (EPI) FID sequence.""" from math import floor from pathlib import Path @@ -204,7 +204,7 @@ def epi2d_fid_kernel( prot.append_acquisition(acq) for slice_ in range(n_slices): - # define label(s) + # define slice label slice_label = pp.make_label(label='SLC', type='SET', value=slice_) # set frequency offset for current slice @@ -227,7 +227,10 @@ def epi2d_fid_kernel( pp.make_label(label='NAV', type='SET', value=1), pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), ) + # reverse gx_pre back after adding to sequence gx_pre = pp.scale_grad(gx_pre, -1) + + # add 3 navigator acquisitions for n in range(3): seq.add_block( gx, @@ -346,12 +349,51 @@ def main( test_report: bool = True, timing_check: bool = True, ) -> tuple[pp.Sequence, Path]: - """Generate an Echo-Planar Imaging (EPI) sequence. + """Generate a 2D Echo Planar Imaging (EPI) FID sequence. + + Parameters + ---------- + system + PyPulseq system limits object. + te + Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. + tr + Desired repetition time (TR) (in seconds). Minimum repetition time is used if set to None. + fov + Field of view in x and y direction (in meters). + n_readout + Number of frequency encoding steps. + n_phase_encoding + Number of phase encoding steps. + n_slices + Number of slices. + slice_thickness + Slice thickness of the 2D slice (in meters). + bandwidth + Total receiver bandwidth (in Hz). + readout_type + Readout type ('symmetric' or 'flyback'). + oversampling + Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. + ramp_sampling + If True, ADC is active during gradient ramps for optimized timing. + partial_fourier_factor + Desired partial Fourier factor in "phase encoding" direction. + add_navigator_acq + If True, navigator acquisitions will be added for ghost corrections. + add_noise_acq + If True, noise acquisitions will be added at the beginning of the sequence. + show_plots + Toggles sequence plot. + test_report + Toggles advanced test report. + timing_check + Toggles timing check of the sequence. Returns ------- seq - Sequence object of radial FLASH sequence. + Sequence object of 2D EPI FID sequence. file_path Path to the sequence file. """ @@ -368,11 +410,11 @@ def main( enable_gradient_spoiling = True # define sequence filename - rs_string = 'rs' if ramp_sampling else 'nors' - pf_string = f'{partial_fourier_factor}pf'.replace('.', 'p') - readout_string = 'sym' if readout_type == 'symmetric' else 'flyb' - noise_string = 'withnoise' if add_noise_acq else 'nonoise' - nav_string = 'withnav' if add_navigator_acq else 'nonav' + rs_string = 'rs' if ramp_sampling else 'nors' # ramp sampling + pf_string = f'{partial_fourier_factor}pf'.replace('.', 'p') # partial fourier factor + readout_string = 'sym' if readout_type == 'symmetric' else 'flyb' # readout type + noise_string = 'withnoise' if add_noise_acq else 'nonoise' # noise acquisition + nav_string = 'withnav' if add_navigator_acq else 'nonav' # navigator acquisition filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}px' filename += f'_{readout_string}_fid_{oversampling}ro_{rs_string}_{pf_string}' @@ -420,7 +462,7 @@ def main( print('\nTiming check failed! Error listing follows\n') print(error_report) - # show advanced rest report + # show advanced test report if test_report: print('\nCreating advanced test report...') print(seq.test_report()) diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index ccd005b..47860a7 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -1,11 +1,10 @@ -"""2D Echo Planar Imaging (EPI) sequence.""" +"""2D Echo Planar Imaging (EPI) spin echo (SE) sequence.""" from math import floor from pathlib import Path from typing import Literal import ismrmrd -import matplotlib.pyplot as plt import numpy as np import pypulseq as pp @@ -73,20 +72,21 @@ def epi2d_se_kernel( Apodization factor of rf excitation pulse readout_type Readout type ('symmetric' or 'flyback'). - echo_type - Echo type ('FID' or 'SE'). oversampling - ADC oversampling factor. + Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. ramp_sampling - If True, ADC is active during gradient ramps (optimized timing). + If True, ADC is active during gradient ramps for optimized timing. partial_fourier_factor - Partial Fourier factor (0.5 to 1.0). + Desired partial Fourier factor in "phase encoding" direction. Must be larger than 0.5 and smaller or equal to 1. + The actual partial Fourier factor might slightly deviate from the desired value. add_spoiler - Enable spoiling gradients. + If True, a spoiler gradient will be added to the sequence after the EPI readout. add_noise_acq - Enable noise acquisition. + If True, noise acquisitions will be added at the beginning of the sequence. add_navigator_acq - Enable navigator acquisition. + If True, 3 navigator acquisitions will be added to the sequence to allow for ghost corrections. + The navigator acquisitions are added between the rf excitation pulse and the refocusing pulse. + Be aware that navigator acquisitions will increase the minimum echo and repetition times. mrd_header_file Filename of the ISMRMRD header file to be created. If None, no header file is created. @@ -266,10 +266,8 @@ def epi2d_se_kernel( acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) prot.append_acquisition(acq) - t_after_noise = sum(seq.block_durations.values()) - for slice_ in range(n_slices): - # define label(s) + # define slice label slice_label = pp.make_label(label='SLC', type='SET', value=slice_) # set frequency offset for current slice @@ -282,8 +280,8 @@ def epi2d_se_kernel( # add navigator scans for ghost correction if add_navigator_acq: # reverse the readout gradient and pre-winder in advance for navigator - gx_pre = pp.scale_grad(epi2d.gx_pre, 1) - gx = pp.scale_grad(epi2d.gx, 1) + gx_pre = pp.scale_grad(epi2d.gx_pre, -1) + gx = pp.scale_grad(epi2d.gx, -1) # add slice selection rewinder and readout pre-winder in x direction (gy_pre will be added after navigators) gzr, gx_pre = pp.align(left=[gzr], right=[gx_pre]) seq.add_block( @@ -292,10 +290,10 @@ def epi2d_se_kernel( pp.make_label(label='NAV', type='SET', value=1), pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), ) + # reverse gx_pre back after adding to sequence gx_pre = pp.scale_grad(gx_pre, -1) - # reverse gx_pre back after addBlock - gx_pre = pp.scale_grad(gx_pre, -1) + # add 3 navigator acquisitions for n in range(3): seq.add_block( gx, @@ -305,6 +303,7 @@ def epi2d_se_kernel( pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), ) gx = pp.scale_grad(gx, -1) + # add navigator acquisitions to ISMRMRD file if mrd_header_file: acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) @@ -390,11 +389,7 @@ def epi2d_se_kernel( seq.set_definition('TR', tr or float(min_tr)) seq.set_definition('ReadoutOversamplingFactor', oversampling) - t_exc = t_after_noise + rf.delay + rf.shape_dur / 2 - t_ref = t_exc + t_exc_to_ref - t_echo = t_exc + t_exc_to_ref + t_ref_to_kcenter - - return seq, min_te, min_tr, t_exc, t_ref, t_echo + return seq, min_te, min_tr def main( @@ -417,12 +412,51 @@ def main( test_report: bool = True, timing_check: bool = True, ) -> tuple[pp.Sequence, Path]: - """Generate an Echo-Planar Imaging (EPI) sequence. + """Generate a 2D Echo Planar Imaging (EPI) spin echo (SE) sequence. + + Parameters + ---------- + system + PyPulseq system limits object. + te + Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. + tr + Desired repetition time (TR) (in seconds). Minimum repetition time is used if set to None. + fov + Field of view in x and y direction (in meters). + n_readout + Number of frequency encoding steps. + n_phase_encoding + Number of phase encoding steps. + n_slices + Number of slices. + slice_thickness + Slice thickness of the 2D slice (in meters). + bandwidth + Total receiver bandwidth (in Hz). + readout_type + Readout type ('symmetric' or 'flyback'). + oversampling + Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. + ramp_sampling + If True, ADC is active during gradient ramps for optimized timing. + partial_fourier_factor + Desired partial Fourier factor in "phase encoding" direction. + add_navigator_acq + If True, navigator acquisitions will be added for ghost corrections. + add_noise_acq + If True, noise acquisitions will be added at the beginning of the sequence. + show_plots + Toggles sequence plot. + test_report + Toggles advanced test report. + timing_check + Toggles timing check of the sequence. Returns ------- seq - Sequence object of radial FLASH sequence. + Sequence object of 2D EPI SE sequence. file_path Path to the sequence file. """ @@ -439,11 +473,11 @@ def main( enable_gradient_spoiling = True # define sequence filename - rs_string = 'rs' if ramp_sampling else 'nors' - pf_string = f'{partial_fourier_factor}pf'.replace('.', 'p') - readout_string = 'sym' if readout_type == 'symmetric' else 'flyb' - noise_string = 'withnoise' if add_noise_acq else 'nonoise' - nav_string = 'withnav' if add_navigator_acq else 'nonav' + rs_string = 'rs' if ramp_sampling else 'nors' # ramp sampling + pf_string = f'{partial_fourier_factor}pf'.replace('.', 'p') # partial fourier factor + readout_string = 'sym' if readout_type == 'symmetric' else 'flyb' # readout type + noise_string = 'withnoise' if add_noise_acq else 'nonoise' # noise acquisition + nav_string = 'withnav' if add_navigator_acq else 'nonav' # navigator acquisition filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}px' filename += f'_{readout_string}_se_{oversampling}ro_{rs_string}_{pf_string}' @@ -458,7 +492,7 @@ def main( mrd_file = output_path / Path(filename + '_header.h5') - seq, _min_te, min_tr, t_exc, t_ref, t_echo = epi2d_se_kernel( + seq, _min_te, min_tr = epi2d_se_kernel( system=system, te=te, tr=tr, @@ -482,17 +516,6 @@ def main( mrd_header_file=mrd_file, ) - # print(f'min_te = {_min_te * 1000:.4f}') - # print(f'min_tr = {min_tr * 1000:.4f}') - # print(f't_exc = {t_exc * 1000:.4f}') - # print(f't_ref = {t_ref * 1000:.4f}') - # print(f't_echo = {t_echo * 1000:.4f}') - # print(f'Echo time (TE) = {(t_echo - t_exc) * 1000:.4f}') - - # print(f'Time from exc to ref: {(t_ref - t_exc) * 1000:.4f}') - # print(f'Time from ref to echo: {(t_echo - t_ref) * 1000:.4f}') - # print(f'added together = {((t_ref - t_exc) + (t_echo - t_ref)) * 1000:.4f}') - # check timing of the sequence if timing_check and not test_report: ok, error_report = seq.check_timing() @@ -502,7 +525,7 @@ def main( print('\nTiming check failed! Error listing follows\n') print(error_report) - # show advanced rest report + # show advanced test report if test_report: print('\nCreating advanced test report...') print(seq.test_report()) @@ -511,18 +534,8 @@ def main( print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") seq.write(str(output_path / filename), create_signature=True) - # # calculate k-space trajectory - # k_traj_adc, k_traj, _, _, _ = seq.calculate_kspace() - - # # plot trajectory - # fig = plt.figure() - # plt.plot(k_traj[0], k_traj[1], 'b') - # plt.plot(k_traj_adc[0], k_traj_adc[1], 'x', color='red', markersize=4) - # plt.grid() - # plt.show() - if show_plots: - seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=False) + seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=True) return seq, output_path / filename From 567cc09a70df2e02c1514a965b63d2d47117bd69 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 25 Feb 2026 12:45:37 +0100 Subject: [PATCH 09/32] fix mrd trajectory definitions for multi-slice --- src/mrseq/scripts/epi2d_fid.py | 27 +++++++++++++++------------ src/mrseq/scripts/epi2d_se.py | 31 +++++++++++++++++-------------- 2 files changed, 32 insertions(+), 26 deletions(-) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index c476ca1..e561ebf 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -103,6 +103,12 @@ def epi2d_fid_kernel( # create PyPulseq Sequence object and set system limits seq = pp.Sequence(system=system) + # define number of navigator acquisitions + n_navigator_acq = 3 if add_navigator_acq else 0 + + # define number of noise acquisitions + n_noise_acq = 1 if add_noise_acq else 0 + # create EpiReadout object epi2d = EpiReadout( system=system, @@ -231,7 +237,7 @@ def epi2d_fid_kernel( gx_pre = pp.scale_grad(gx_pre, -1) # add 3 navigator acquisitions - for n in range(3): + for n in range(n_navigator_acq): seq.add_block( gx, epi2d.adc, @@ -274,20 +280,17 @@ def epi2d_fid_kernel( # calculate k-space trajectory from Sequence to add this info to mrd file k_traj_adc, _, _, _, _ = seq.calculate_kspace() samples_per_acq = epi2d.adc.num_samples - number_of_total_acq = k_traj_adc.shape[-1] // samples_per_acq - number_of_epi_acq = epi2d.n_phase_enc_total - if 'data' in prot._dataset: - number_of_noise_acq = prot.number_of_acquisitions() - else: - number_of_noise_acq = 0 + n_total_acq = k_traj_adc.shape[-1] // samples_per_acq + n_epi_acq_per_slice = epi2d.n_phase_enc_total - if not number_of_epi_acq + number_of_noise_acq == number_of_total_acq: - raise (ValueError('Number of calculated acquisitions does not match expected number.')) + if not (n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq == n_total_acq: + raise ValueError('Number of calculated acquisitions does not match expected number.') - k_traj_adc_readout = k_traj_adc[:, number_of_noise_acq * samples_per_acq :] + # remove noise acquisitions from k-space trajectory + k_traj_adc_readout = k_traj_adc[:, n_noise_acq * samples_per_acq :] - # create mrd acquisition and add trajectory info for each EPI readout - for n in range(number_of_epi_acq): + # create mrd acquisition and add trajectory info for each EPI readout including navigator acquisitions + for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices): start = n * samples_per_acq end = start + samples_per_acq diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 47860a7..9a13135 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -103,6 +103,12 @@ def epi2d_se_kernel( # create PyPulseq Sequence object and set system limits seq = pp.Sequence(system=system) + # define number of navigator acquisitions + n_navigator_acq = 3 if add_navigator_acq else 0 + + # define number of noise acquisitions + n_noise_acq = 1 if add_noise_acq else 0 + # create EpiReadout object epi2d = EpiReadout( system=system, @@ -294,7 +300,7 @@ def epi2d_se_kernel( gx_pre = pp.scale_grad(gx_pre, -1) # add 3 navigator acquisitions - for n in range(3): + for n in range(n_navigator_acq): seq.add_block( gx, epi2d.adc, @@ -337,20 +343,17 @@ def epi2d_se_kernel( # calculate k-space trajectory from Sequence to add this info to mrd file k_traj_adc, _, _, _, _ = seq.calculate_kspace() samples_per_acq = epi2d.adc.num_samples - number_of_total_acq = k_traj_adc.shape[-1] // samples_per_acq - number_of_epi_acq = epi2d.n_phase_enc_total - if 'data' in prot._dataset: - number_of_noise_acq = prot.number_of_acquisitions() - else: - number_of_noise_acq = 0 + n_total_acq = k_traj_adc.shape[-1] // samples_per_acq + n_epi_acq_per_slice = epi2d.n_phase_enc_total - if not number_of_epi_acq + number_of_noise_acq == number_of_total_acq: - raise (ValueError('Number of calculated acquisitions does not match expected number.')) + if not (n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq == n_total_acq: + raise ValueError('Number of calculated acquisitions does not match expected number.') - k_traj_adc_readout = k_traj_adc[:, number_of_noise_acq * samples_per_acq :] + # remove noise acquisitions from k-space trajectory + k_traj_adc_readout = k_traj_adc[:, n_noise_acq * samples_per_acq :] - # create mrd acquisition and add trajectory info for each EPI readout - for n in range(number_of_epi_acq): + # create mrd acquisition and add trajectory info for each EPI readout including navigator acquisitions + for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices): start = n * samples_per_acq end = start + samples_per_acq @@ -399,13 +402,13 @@ def main( fov: float = 200e-3, n_readout: int = 64, n_phase_encoding: int = 64, - n_slices: int = 1, + n_slices: int = 3, slice_thickness: float = 8e-3, bandwidth: float = 100e3, readout_type: Literal['symmetric', 'flyback'] = 'symmetric', oversampling: Literal[1, 2, 4] = 2, ramp_sampling: bool = True, - partial_fourier_factor: float = 0.7, + partial_fourier_factor: float = 1.0, add_navigator_acq: bool = True, add_noise_acq: bool = True, show_plots: bool = True, From 0de08077ca548198723c427303580b119679aa1e Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Sat, 28 Feb 2026 20:53:32 +0000 Subject: [PATCH 10/32] tests added --- src/mrseq/utils/trajectory.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py index b183956..e4ab95a 100644 --- a/src/mrseq/utils/trajectory.py +++ b/src/mrseq/utils/trajectory.py @@ -714,11 +714,7 @@ def __init__( # recalculate the actual partial fourier factor actual_pf_factor = self.n_phase_enc_total / self.n_phase_encoding if actual_pf_factor != partial_fourier_factor: - warnings.warn( - f'Desired partial Fourier factor {partial_fourier_factor} adjusted to {actual_pf_factor}.', - UserWarning, - stacklevel=2, - ) + print(f'Adjusted partial Fourier factor from {partial_fourier_factor} to {actual_pf_factor:.2f}.') self.partial_fourier_factor = actual_pf_factor # Create pre-phaser gradients From 1ebea4dcd49ae7895ba7f36b4a549fe9d7292e3f Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Sat, 28 Feb 2026 22:27:52 +0000 Subject: [PATCH 11/32] tests and example added. ismrmrd header fixed --- examples/epi2d.ipynb | 410 ++++++++++++++++++++++++++++++++ src/mrseq/scripts/epi2d_fid.py | 69 +++--- src/mrseq/scripts/epi2d_se.py | 63 ++--- tests/scripts/test_epi2d_fid.py | 25 ++ tests/scripts/test_epi2d_se.py | 25 ++ 5 files changed, 513 insertions(+), 79 deletions(-) create mode 100644 examples/epi2d.ipynb create mode 100644 tests/scripts/test_epi2d_fid.py create mode 100644 tests/scripts/test_epi2d_se.py diff --git a/examples/epi2d.ipynb b/examples/epi2d.ipynb new file mode 100644 index 0000000..689ed76 --- /dev/null +++ b/examples/epi2d.ipynb @@ -0,0 +1,410 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Echo Planar Imaging\n", + "\n", + "Qualitative imaging with an EPI readout." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import tempfile\n", + "from pathlib import Path\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import MRzeroCore as mr0\n", + "import numpy as np\n", + "from mrpro.algorithms.reconstruction import DirectReconstruction\n", + "from mrpro.data import KData\n", + "from mrpro.data.traj_calculators import KTrajectoryIsmrmrd\n", + "from mrpro.operators import NonUniformFastFourierOp\n", + "\n", + "from mrseq.scripts.epi2d_fid import main as create_seq_fid\n", + "from mrseq.scripts.epi2d_se import main as create_seq_se\n", + "from mrseq.utils import combine_ismrmrd_files\n", + "from mrseq.utils import sys_defaults\n", + "from mrseq.utils.trajectory import EpiReadout" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Settings\n", + "We are going to use a numerical phantom with a matrix size of 128 x 128. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "image_matrix_size = [128, 128]\n", + "flip_angle_degree = 12\n", + "n_dummy_excitations = 200\n", + "\n", + "tmp = tempfile.TemporaryDirectory()" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "### Create the digital phantom\n", + "\n", + "We use the standard Brainweb phantom from [MRzero](https://github.com/MRsources/MRzero-Core), but we choose the B1-field to be constant everywhere." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "phantom = mr0.util.load_phantom(image_matrix_size)\n", + "phantom.B1[:] = 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "### EPI readouts\n", + "\n", + "Single-shot EPI readouts are a very efficient way to obtain a full image after a single excitation. \n", + "There are different options of how they can be carried out. " + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "We can have a symmetric readout where we alternate between going from left to right and right to left through k-space." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "epi_readout = EpiReadout(n_readout=12, n_phase_encoding=12, ramp_sampling=False, readout_type='symmetric')\n", + "fig = epi_readout.plot_trajectory()" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "Here we are sampling only during the flat time of the readout gradient. \n", + "We can make the sequence faster by also sampling during the ramp-up and ramp-down part of the readout gradient.\n", + "Nevertheless, this means that we don't have a strictly Cartesian trajectory anymore and gridding/nufft is required \n", + "during image reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "epi_readout = EpiReadout(n_readout=12, n_phase_encoding=12, ramp_sampling=True, readout_type='symmetric')\n", + "fig = epi_readout.plot_trajectory()" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "We can shorten the echo time by applying partial Fourier along the phase encoding direction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "epi_readout = EpiReadout(\n", + " n_readout=12, n_phase_encoding=12, ramp_sampling=True, partial_fourier_factor=0.75, readout_type='symmetric'\n", + ")\n", + "fig = epi_readout.plot_trajectory()" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "Another approach is to always go back to one side of k-space after each readout such that the readout direction is always from left to right.\n", + "This is less efficient but minimized any artifacts due to asymmetry between positive and negative readout gradients. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "epi_readout = EpiReadout(\n", + " n_readout=12, n_phase_encoding=12, partial_fourier_factor=0.75, ramp_sampling=True, readout_type='flyback'\n", + ")\n", + "fig = epi_readout.plot_trajectory()" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "### EPI sequences\n", + "\n", + "To create the EPI FID and SE sequences, we use the previously imported [epi2d_fid script](../src/mrseq/scripts/ei2d_fid.py) and [epi2d_se script](../src/mrseq/scripts/ei2d_se.py).\n", + "\n", + "We are going to simulate different types of EPI sequences and compare the obtained images." + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "#### EPI FID with symmetric readout without ramp sampling" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "sequence, fname_seq = create_seq_fid(\n", + " system=sys_defaults,\n", + " test_report=False,\n", + " timing_check=False,\n", + " show_plots=False,\n", + " fov_xy=float(phantom.size.numpy()[0]),\n", + " n_readout=image_matrix_size[0],\n", + " n_phase_encoding=image_matrix_size[0],\n", + " partial_fourier_factor=0.75,\n", + " ramp_sampling=False,\n", + " readout_type='symmetric',\n", + ")\n", + "\n", + "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)\n", + "fname_mrd = Path(tmp.name) / 'epi_fid_sym_no_ramp.mrd'\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)\n", + "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", + "\n", + "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", + "recon = DirectReconstruction(kdata, csm=None)\n", + "idata_fid_sym_no_ramp = recon(kdata).data" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "#### EPI FID with symmetric readout and ramp sampling" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "sequence, fname_seq = create_seq_fid(\n", + " system=sys_defaults,\n", + " test_report=False,\n", + " timing_check=False,\n", + " show_plots=False,\n", + " fov_xy=float(phantom.size.numpy()[0]),\n", + " n_readout=image_matrix_size[0],\n", + " n_phase_encoding=image_matrix_size[0],\n", + " partial_fourier_factor=0.75,\n", + " readout_type='symmetric',\n", + ")\n", + "\n", + "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)\n", + "fname_mrd = Path(tmp.name) / 'epi_fid_sym_ramp.mrd'\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)\n", + "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", + "\n", + "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", + "fourier_op = NonUniformFastFourierOp(\n", + " direction=(-2, -1),\n", + " recon_matrix=kdata.header.recon_matrix[-2:],\n", + " encoding_matrix=kdata.header.encoding_matrix[-2:],\n", + " traj=kdata.traj,\n", + ")\n", + "idata_fid_sym_ramp = fourier_op.H(kdata.data)[0]" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "#### EPI FID with flyback readout and ramp sampling" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "sequence, fname_seq = create_seq_fid(\n", + " system=sys_defaults,\n", + " test_report=False,\n", + " timing_check=False,\n", + " show_plots=False,\n", + " fov_xy=float(phantom.size.numpy()[0]),\n", + " n_readout=image_matrix_size[0],\n", + " n_phase_encoding=image_matrix_size[0],\n", + " partial_fourier_factor=0.75,\n", + " readout_type='flyback',\n", + ")\n", + "\n", + "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e0)\n", + "fname_mrd = Path(tmp.name) / 'epi_fid_flyback_ramp.mrd'\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)\n", + "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", + "\n", + "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", + "fourier_op = NonUniformFastFourierOp(\n", + " direction=(-2, -1),\n", + " recon_matrix=kdata.header.recon_matrix[-2:],\n", + " encoding_matrix=kdata.header.encoding_matrix[-2:],\n", + " traj=kdata.traj,\n", + ")\n", + "idata_fid_flybackramp = fourier_op.H(kdata.data)[0]" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "#### EPI SE with symmetric readout and ramp sampling" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "sequence, fname_seq = create_seq_se(\n", + " system=sys_defaults,\n", + " test_report=False,\n", + " timing_check=False,\n", + " show_plots=False,\n", + " fov_xy=float(phantom.size.numpy()[0]),\n", + " n_readout=image_matrix_size[0],\n", + " n_phase_encoding=image_matrix_size[0],\n", + " partial_fourier_factor=0.75,\n", + " readout_type='symmetric',\n", + ")\n", + "\n", + "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-1)\n", + "fname_mrd = Path(tmp.name) / 'epi_se_sym_ramp.mrd'\n", + "mr0.sig_to_mrd(fname_mrd, signal, sequence)\n", + "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", + "\n", + "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", + "fourier_op = NonUniformFastFourierOp(\n", + " direction=(-2, -1),\n", + " recon_matrix=kdata.header.recon_matrix[-2:],\n", + " encoding_matrix=kdata.header.encoding_matrix[-2:],\n", + " traj=kdata.traj,\n", + ")\n", + "idata_se_sym_ramp = fourier_op.H(kdata.data)[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 4, figsize=(16, 8))\n", + "for cax in ax.flatten():\n", + " cax.set_xticks([])\n", + " cax.set_yticks([])\n", + "\n", + "for idx, (idat, label) in enumerate(\n", + " zip(\n", + " [idata_fid_sym_no_ramp, idata_fid_sym_ramp, idata_fid_flybackramp],\n", + " ['symmetric no ramp', 'symmetric ramp', 'flyback ramp'],\n", + " strict=True,\n", + " )\n", + "):\n", + " idat = idat.abs().squeeze().numpy()\n", + " idat /= idat.max()\n", + "\n", + " if idx == 0:\n", + " idat_ref = idat\n", + " relative_error = 0.0\n", + " relative_error += np.sum(np.abs(idat - idat_ref)) / np.sum(np.abs(idat_ref))\n", + "\n", + " ax[0, idx].imshow(idat, vmin=0, vmax=1, cmap='gray')\n", + " ax[0, idx].set_title(label)\n", + " ax[1, idx].imshow(idat - idat_ref, vmin=-1, vmax=1, cmap='bwr')\n", + "\n", + "print(f'Relative error {relative_error}')\n", + "assert relative_error < 0.05" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mrseq_mrzero", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index e561ebf..0c45509 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -21,7 +21,7 @@ def epi2d_fid_kernel( system: pp.Opts, te: float | None, tr: float | None, - fov: float, + fov_xy: float, n_readout: int, n_phase_encoding: int, bandwidth: float, @@ -32,7 +32,7 @@ def epi2d_fid_kernel( rf_bwt: float, rf_apodization: float, readout_type: Literal['symmetric', 'flyback'], - oversampling: Literal[1, 2, 4], + readout_oversampling: Literal[1, 2, 4], ramp_sampling: bool, partial_fourier_factor: float, add_spoiler: bool, @@ -50,7 +50,7 @@ def epi2d_fid_kernel( Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. tr Desired repetition time (TR) (in seconds). - fov + fov_xy Field of view in x and y direction (in meters). n_readout Number of frequency encoding steps. @@ -72,7 +72,7 @@ def epi2d_fid_kernel( Apodization factor of rf excitation pulse. readout_type Readout type ('symmetric' or 'flyback'). - oversampling + readout_oversampling Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. ramp_sampling If True, ADC is active during gradient ramps for optimized timing. @@ -112,11 +112,11 @@ def epi2d_fid_kernel( # create EpiReadout object epi2d = EpiReadout( system=system, - fov=fov, + fov=fov_xy, n_readout=n_readout, n_phase_encoding=n_phase_encoding, bandwidth=bandwidth, - oversampling=oversampling, + oversampling=readout_oversampling, readout_type=readout_type, ramp_sampling=ramp_sampling, partial_fourier_factor=partial_fourier_factor, @@ -177,9 +177,9 @@ def epi2d_fid_kernel( if mrd_header_file: hdr = create_header( traj_type='other', - encoding_fov=Fov(x=fov * oversampling, y=fov, z=slice_thickness), - recon_fov=Fov(x=fov, y=fov, z=slice_thickness), - encoding_matrix=MatrixSize(n_x=n_readout, n_y=epi2d.n_phase_encoding, n_z=1), + encoding_fov=Fov(x=fov_xy * readout_oversampling, y=fov_xy, z=slice_thickness), + recon_fov=Fov(x=fov_xy, y=fov_xy, z=slice_thickness), + encoding_matrix=MatrixSize(n_x=n_readout * readout_oversampling, n_y=epi2d.n_phase_encoding, n_z=1), recon_matrix=MatrixSize(n_x=n_readout, n_y=n_phase_encoding, n_z=1), dwell_time=epi2d.adc.dwell, slice_limits=Limits(min=0, max=n_slices, center=n_slices // 2), @@ -204,11 +204,6 @@ def epi2d_fid_kernel( seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) seq.add_block(pp.make_delay(system.rf_dead_time)) - if mrd_header_file: - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) - prot.append_acquisition(acq) - for slice_ in range(n_slices): # define slice label slice_label = pp.make_label(label='SLC', type='SET', value=slice_) @@ -236,7 +231,7 @@ def epi2d_fid_kernel( # reverse gx_pre back after adding to sequence gx_pre = pp.scale_grad(gx_pre, -1) - # add 3 navigator acquisitions + # add navigator scans for ghost correction for n in range(n_navigator_acq): seq.add_block( gx, @@ -246,11 +241,6 @@ def epi2d_fid_kernel( pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), ) gx = pp.scale_grad(gx, -1) - # add navigator acquisitions to ISMRMRD file - if mrd_header_file: - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) - prot.append_acquisition(acq) # add echo time delay seq.add_block(pp.make_delay(te_delay)) @@ -286,17 +276,14 @@ def epi2d_fid_kernel( if not (n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq == n_total_acq: raise ValueError('Number of calculated acquisitions does not match expected number.') - # remove noise acquisitions from k-space trajectory - k_traj_adc_readout = k_traj_adc[:, n_noise_acq * samples_per_acq :] - - # create mrd acquisition and add trajectory info for each EPI readout including navigator acquisitions - for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices): + # create mrd acquisition and add trajectory info for each EPI readout including navigator and noise acquisitions + for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq): start = n * samples_per_acq end = start + samples_per_acq traj = np.zeros((samples_per_acq, 2), dtype=np.float32) - traj[:, 0] = np.round(k_traj_adc_readout[0, start:end] * fov, 3) - traj[:, 1] = np.round(k_traj_adc_readout[1, start:end] * fov, 3) + traj[:, 0] = np.round(k_traj_adc[0, start:end] * fov_xy * readout_oversampling, 3) + traj[:, 1] = np.round(k_traj_adc[1, start:end] * fov_xy, 3) acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=samples_per_acq) @@ -318,32 +305,32 @@ def epi2d_fid_kernel( epi2d.adc.num_samples * epi2d.adc.dwell, ] - seq.set_definition(key='TargetGriddedSamples', value=n_readout * oversampling) + seq.set_definition(key='TargetGriddedSamples', value=n_readout * readout_oversampling) seq.set_definition(key='TrapezoidGriddingParameters', value=gridding_params) # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov, fov, slice_thickness * n_slices]) + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness * n_slices]) seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, n_slices)) seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr or min_tr) - seq.set_definition('ReadoutOversamplingFactor', oversampling) + seq.set_definition('TE', te or float(min_te)) + seq.set_definition('TR', tr or float(min_tr)) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) - return seq, min_te, min_tr + return seq, float(min_te), float(min_tr) def main( system: pp.Opts | None = None, te: float | None = None, tr: float | None = None, - fov: float = 200e-3, + fov_xy: float = 200e-3, n_readout: int = 64, n_phase_encoding: int = 64, n_slices: int = 1, slice_thickness: float = 8e-3, bandwidth: float = 100e3, readout_type: Literal['symmetric', 'flyback'] = 'symmetric', - oversampling: Literal[1, 2, 4] = 2, + readout_oversampling: Literal[1, 2, 4] = 2, ramp_sampling: bool = True, partial_fourier_factor: float = 0.7, add_navigator_acq: bool = True, @@ -362,7 +349,7 @@ def main( Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. tr Desired repetition time (TR) (in seconds). Minimum repetition time is used if set to None. - fov + fov_xy Field of view in x and y direction (in meters). n_readout Number of frequency encoding steps. @@ -376,7 +363,7 @@ def main( Total receiver bandwidth (in Hz). readout_type Readout type ('symmetric' or 'flyback'). - oversampling + readout_oversampling Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. ramp_sampling If True, ADC is active during gradient ramps for optimized timing. @@ -419,8 +406,8 @@ def main( noise_string = 'withnoise' if add_noise_acq else 'nonoise' # noise acquisition nav_string = 'withnav' if add_navigator_acq else 'nonav' # navigator acquisition - filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}px' - filename += f'_{readout_string}_fid_{oversampling}ro_{rs_string}_{pf_string}' + filename = f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}px' + filename += f'_{readout_string}_fid_{readout_oversampling}ro_{rs_string}_{pf_string}' filename += f'_{noise_string}_{nav_string}' output_path = Path.cwd() / 'output' @@ -436,11 +423,11 @@ def main( system=system, te=te, tr=tr, - fov=fov, + fov_xy=fov_xy, n_readout=n_readout, n_phase_encoding=n_phase_encoding, bandwidth=bandwidth, - oversampling=oversampling, + readout_oversampling=readout_oversampling, slice_thickness=slice_thickness, n_slices=n_slices, rf_duration=rf_duration, diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 9a13135..2d89c63 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -21,7 +21,7 @@ def epi2d_se_kernel( system: pp.Opts, te: float | None, tr: float | None, - fov: float, + fov_xy: float, n_readout: int, n_phase_encoding: int, bandwidth: float, @@ -32,7 +32,7 @@ def epi2d_se_kernel( rf_bwt: float, rf_apodization: float, readout_type: Literal['symmetric', 'flyback'], - oversampling: Literal[1, 2, 4], + readout_oversampling: Literal[1, 2, 4], ramp_sampling: bool, partial_fourier_factor: float, add_spoiler: bool, @@ -50,7 +50,7 @@ def epi2d_se_kernel( Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. tr Desired repetition time (TR) (in seconds). - fov + fov_xy Field of view in x and y direction (in meters). n_readout Number of frequency encoding steps. @@ -72,7 +72,7 @@ def epi2d_se_kernel( Apodization factor of rf excitation pulse readout_type Readout type ('symmetric' or 'flyback'). - oversampling + readout_oversampling Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. ramp_sampling If True, ADC is active during gradient ramps for optimized timing. @@ -112,11 +112,11 @@ def epi2d_se_kernel( # create EpiReadout object epi2d = EpiReadout( system=system, - fov=fov, + fov=fov_xy, n_readout=n_readout, n_phase_encoding=n_phase_encoding, bandwidth=bandwidth, - oversampling=oversampling, + oversampling=readout_oversampling, readout_type=readout_type, ramp_sampling=ramp_sampling, partial_fourier_factor=partial_fourier_factor, @@ -240,9 +240,9 @@ def epi2d_se_kernel( if mrd_header_file: hdr = create_header( traj_type='other', - encoding_fov=Fov(x=fov * oversampling, y=fov, z=slice_thickness), - recon_fov=Fov(x=fov, y=fov, z=slice_thickness), - encoding_matrix=MatrixSize(n_x=n_readout, n_y=epi2d.n_phase_encoding, n_z=1), + encoding_fov=Fov(x=fov_xy * readout_oversampling, y=fov_xy, z=slice_thickness), + recon_fov=Fov(x=fov_xy, y=fov_xy, z=slice_thickness), + encoding_matrix=MatrixSize(n_x=n_readout * readout_oversampling, n_y=epi2d.n_phase_encoding, n_z=1), recon_matrix=MatrixSize(n_x=n_readout, n_y=n_phase_encoding, n_z=1), dwell_time=epi2d.adc.dwell, slice_limits=Limits(min=0, max=n_slices, center=n_slices // 2), @@ -267,11 +267,6 @@ def epi2d_se_kernel( seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) seq.add_block(pp.make_delay(system.rf_dead_time)) - if mrd_header_file: - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) - prot.append_acquisition(acq) - for slice_ in range(n_slices): # define slice label slice_label = pp.make_label(label='SLC', type='SET', value=slice_) @@ -309,11 +304,6 @@ def epi2d_se_kernel( pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), ) gx = pp.scale_grad(gx, -1) - # add navigator acquisitions to ISMRMRD file - if mrd_header_file: - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=epi2d.adc.num_samples) - prot.append_acquisition(acq) # add gy_pre and reset labels seq.add_block( @@ -349,17 +339,14 @@ def epi2d_se_kernel( if not (n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq == n_total_acq: raise ValueError('Number of calculated acquisitions does not match expected number.') - # remove noise acquisitions from k-space trajectory - k_traj_adc_readout = k_traj_adc[:, n_noise_acq * samples_per_acq :] - - # create mrd acquisition and add trajectory info for each EPI readout including navigator acquisitions - for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices): + # create mrd acquisition and add trajectory info for each EPI readout including navigator and noise acquisitions + for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq): start = n * samples_per_acq end = start + samples_per_acq traj = np.zeros((samples_per_acq, 2), dtype=np.float32) - traj[:, 0] = np.round(k_traj_adc_readout[0, start:end] * fov, 3) - traj[:, 1] = np.round(k_traj_adc_readout[1, start:end] * fov, 3) + traj[:, 0] = np.round(k_traj_adc[0, start:end] * fov_xy * readout_oversampling, 3) + traj[:, 1] = np.round(k_traj_adc[1, start:end] * fov_xy, 3) acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=samples_per_acq) @@ -381,32 +368,32 @@ def epi2d_se_kernel( epi2d.adc.num_samples * epi2d.adc.dwell, ] - seq.set_definition(key='TargetGriddedSamples', value=n_readout * oversampling) + seq.set_definition(key='TargetGriddedSamples', value=n_readout * readout_oversampling) seq.set_definition(key='TrapezoidGriddingParameters', value=gridding_params) # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov, fov, slice_thickness * n_slices]) + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness * n_slices]) seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, n_slices)) seq.set_definition('SliceThickness', slice_thickness) seq.set_definition('TE', float(t_exc_to_ref + t_ref_to_kcenter)) seq.set_definition('TR', tr or float(min_tr)) - seq.set_definition('ReadoutOversamplingFactor', oversampling) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) - return seq, min_te, min_tr + return seq, float(min_te), float(min_tr) def main( system: pp.Opts | None = None, te: float | None = None, tr: float | None = None, - fov: float = 200e-3, + fov_xy: float = 200e-3, n_readout: int = 64, n_phase_encoding: int = 64, n_slices: int = 3, slice_thickness: float = 8e-3, bandwidth: float = 100e3, readout_type: Literal['symmetric', 'flyback'] = 'symmetric', - oversampling: Literal[1, 2, 4] = 2, + readout_oversampling: Literal[1, 2, 4] = 2, ramp_sampling: bool = True, partial_fourier_factor: float = 1.0, add_navigator_acq: bool = True, @@ -425,7 +412,7 @@ def main( Desired echo time (TE) (in seconds). Minimum echo time is used if set to None. tr Desired repetition time (TR) (in seconds). Minimum repetition time is used if set to None. - fov + fov_xy Field of view in x and y direction (in meters). n_readout Number of frequency encoding steps. @@ -439,7 +426,7 @@ def main( Total receiver bandwidth (in Hz). readout_type Readout type ('symmetric' or 'flyback'). - oversampling + readout_oversampling Readout oversampling factor. Can be 1 (no oversampling), 2, or 4. ramp_sampling If True, ADC is active during gradient ramps for optimized timing. @@ -482,8 +469,8 @@ def main( noise_string = 'withnoise' if add_noise_acq else 'nonoise' # noise acquisition nav_string = 'withnav' if add_navigator_acq else 'nonav' # navigator acquisition - filename = f'{Path(__file__).stem}_{int(fov * 1000)}fov_{n_readout}px' - filename += f'_{readout_string}_se_{oversampling}ro_{rs_string}_{pf_string}' + filename = f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}px' + filename += f'_{readout_string}_se_{readout_oversampling}ro_{rs_string}_{pf_string}' filename += f'_{noise_string}_{nav_string}' output_path = Path.cwd() / 'output' @@ -499,11 +486,11 @@ def main( system=system, te=te, tr=tr, - fov=fov, + fov_xy=fov_xy, n_readout=n_readout, n_phase_encoding=n_phase_encoding, bandwidth=bandwidth, - oversampling=oversampling, + readout_oversampling=readout_oversampling, slice_thickness=slice_thickness, n_slices=n_slices, rf_duration=rf_duration, diff --git a/tests/scripts/test_epi2d_fid.py b/tests/scripts/test_epi2d_fid.py new file mode 100644 index 0000000..36953e8 --- /dev/null +++ b/tests/scripts/test_epi2d_fid.py @@ -0,0 +1,25 @@ +"""Tests for EPI FID sequence.""" + +import pytest +from mrseq.scripts.epi2d_fid import main as create_seq + +EXPECTED_DUR = 0.0389 # defined 2026-02-28 + + +def test_default_seq_duration(system_defaults): + """Test if default values result in expected sequence duration.""" + seq, _ = create_seq(system=system_defaults, show_plots=False) + duration = seq.duration()[0] + assert duration == pytest.approx(EXPECTED_DUR) + + +def test_seq_creation_error_on_short_te(system_defaults): + """Test if error is raised on too short echo time.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, te=1e-3, show_plots=False) + + +def test_seq_creation_error_on_short_tr(system_defaults): + """Test if error is raised on too short repetition time.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, tr=2e-3, show_plots=False) diff --git a/tests/scripts/test_epi2d_se.py b/tests/scripts/test_epi2d_se.py new file mode 100644 index 0000000..9a2000b --- /dev/null +++ b/tests/scripts/test_epi2d_se.py @@ -0,0 +1,25 @@ +"""Tests for EPI SE sequence.""" + +import pytest +from mrseq.scripts.epi2d_se import main as create_seq + +EXPECTED_DUR = 0.22586 # defined 2026-02-28 + + +def test_default_seq_duration(system_defaults): + """Test if default values result in expected sequence duration.""" + seq, _ = create_seq(system=system_defaults, show_plots=False) + duration = seq.duration()[0] + assert duration == pytest.approx(EXPECTED_DUR) + + +def test_seq_creation_error_on_short_te(system_defaults): + """Test if error is raised on too short echo time.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, te=1e-3, show_plots=False) + + +def test_seq_creation_error_on_short_tr(system_defaults): + """Test if error is raised on too short repetition time.""" + with pytest.raises(ValueError): + create_seq(system=system_defaults, tr=2e-3, show_plots=False) From 81338c8e219ceadf0855d7b7dade5ecd3cdb52a0 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Sun, 1 Mar 2026 22:16:25 +0000 Subject: [PATCH 12/32] separate epi readout + tests --- examples/epi2d.ipynb | 20 +- src/mrseq/scripts/epi2d_fid.py | 2 +- src/mrseq/scripts/epi2d_se.py | 2 +- src/mrseq/utils/EpiReadout.py | 368 ++++++++++++++++++++++++++++++ src/mrseq/utils/__init__.py | 1 + src/mrseq/utils/trajectory.py | 359 ----------------------------- tests/utils/gradient_waveforms.py | 42 ++++ tests/utils/test_epi_readout.py | 67 ++++++ tests/utils/test_trajectory.py | 40 +--- 9 files changed, 498 insertions(+), 403 deletions(-) create mode 100644 src/mrseq/utils/EpiReadout.py create mode 100644 tests/utils/gradient_waveforms.py create mode 100644 tests/utils/test_epi_readout.py diff --git a/examples/epi2d.ipynb b/examples/epi2d.ipynb index 689ed76..05646cf 100644 --- a/examples/epi2d.ipynb +++ b/examples/epi2d.ipynb @@ -214,6 +214,7 @@ " test_report=False,\n", " timing_check=False,\n", " show_plots=False,\n", + " te=0.073,\n", " fov_xy=float(phantom.size.numpy()[0]),\n", " n_readout=image_matrix_size[0],\n", " n_phase_encoding=image_matrix_size[0],\n", @@ -253,6 +254,7 @@ " test_report=False,\n", " timing_check=False,\n", " show_plots=False,\n", + " te=0.073,\n", " fov_xy=float(phantom.size.numpy()[0]),\n", " n_readout=image_matrix_size[0],\n", " n_phase_encoding=image_matrix_size[0],\n", @@ -296,6 +298,7 @@ " test_report=False,\n", " timing_check=False,\n", " show_plots=False,\n", + " te=0.073,\n", " fov_xy=float(phantom.size.numpy()[0]),\n", " n_readout=image_matrix_size[0],\n", " n_phase_encoding=image_matrix_size[0],\n", @@ -347,7 +350,7 @@ ")\n", "\n", "mr0_sequence = mr0.Sequence.import_file(str(fname_seq.with_suffix('.seq')))\n", - "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-1)\n", + "signal, ktraj_adc = mr0.util.simulate(mr0_sequence, phantom, accuracy=1e-6)\n", "fname_mrd = Path(tmp.name) / 'epi_se_sym_ramp.mrd'\n", "mr0.sig_to_mrd(fname_mrd, signal, sequence)\n", "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", @@ -382,8 +385,7 @@ " )\n", "):\n", " idat = idat.abs().squeeze().numpy()\n", - " idat /= idat.max()\n", - "\n", + " idat /= np.sort(idat.flatten())[int(idat.size * 0.99)]\n", " if idx == 0:\n", " idat_ref = idat\n", " relative_error = 0.0\n", @@ -393,9 +395,19 @@ " ax[0, idx].set_title(label)\n", " ax[1, idx].imshow(idat - idat_ref, vmin=-1, vmax=1, cmap='bwr')\n", "\n", + "relative_error /= 3\n", + "\n", "print(f'Relative error {relative_error}')\n", - "assert relative_error < 0.05" + "assert relative_error < 0.15" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index 0c45509..042d0a5 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -10,11 +10,11 @@ from mrseq.utils import round_to_raster from mrseq.utils import sys_defaults +from mrseq.utils.EpiReadout import EpiReadout from mrseq.utils.ismrmrd import Fov from mrseq.utils.ismrmrd import Limits from mrseq.utils.ismrmrd import MatrixSize from mrseq.utils.ismrmrd import create_header -from mrseq.utils.trajectory import EpiReadout def epi2d_fid_kernel( diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 2d89c63..1802436 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -10,11 +10,11 @@ from mrseq.utils import round_to_raster from mrseq.utils import sys_defaults +from mrseq.utils.EpiReadout import EpiReadout from mrseq.utils.ismrmrd import Fov from mrseq.utils.ismrmrd import Limits from mrseq.utils.ismrmrd import MatrixSize from mrseq.utils.ismrmrd import create_header -from mrseq.utils.trajectory import EpiReadout def epi2d_se_kernel( diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py new file mode 100644 index 0000000..e14395c --- /dev/null +++ b/src/mrseq/utils/EpiReadout.py @@ -0,0 +1,368 @@ +"""Echo Planar Imaging (EPI) readout sequence.""" + +from typing import Literal + +import ismrmrd +import matplotlib.pyplot as plt +import numpy as np +import pypulseq as pp + +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults + + +class EpiReadout: + """EPI readout module supporting flyback and symmetric readout, ramp sampling, oversampling, and partial Fourier. + + Attributes + ---------- + system + System limits. + ramp_sampling + If True, ADC is active during gradient ramps (optimized timing). + readout_type + Readout type ('flyback' or 'symmetric'). + spoiling_enable + Enable spoiling gradients (useful for calibration if False). + adc + ADC event. + gx_pre + Pre-phaser gradient in readout direction. + gy_pre + Pre-phaser gradient in phase encoding direction. + gx + Readout gradient. + gy + Phase encoding gradient. + gy_blip + Complete blip gradient used for flyback readout. + gy_blipup + Ramp-up part of blip gradient used for symmetric readout. + gy_blipdown + Ramp-down part of blip gradient used for symmetric readout. + gy_blipdownup + Combined gy_blipdown and gy_blipup gradient used for symmetric readout. + gz_spoil + Spoiler gradient in z-direction. + gx_flyback + Flyback gradient in x-direction. + """ + + def __init__( + self, + system: pp.Opts | None = None, + fov: float = 128e-3, + n_readout: int = 64, + n_phase_encoding: int = 64, + bandwidth: float = 50e3, + oversampling: int = 2, + ramp_sampling: bool = True, + readout_type: Literal['flyback', 'symmetric'] = 'symmetric', + partial_fourier_factor: float = 1.0, + adc_freq_offset: float = 0, + pe_enable: bool = True, + spoiling_enable: bool = True, + ): + """Initialize EPI Readout. + + Parameters + ---------- + system + System limits. + fov + Field of view in meters (square). + n_readout + Number of readout points. + n_phase_encoding + Number of phase encoding points. + bandwidth + Total receiver bandwidth in Hz. + oversampling + ADC oversampling factor. + ramp_sampling + If True, ADC is active during gradient ramps (optimized timing). + readout_type + Readout type ('flyback' or 'symmetric'). + partial_fourier_factor + Partial Fourier factor (0.5 to 1.0). + adc_freq_offset + Frequency offset for the ADC. + pe_enable + Enable phase encoding (useful for calibration scans if False). + spoiling_enable + Enable spoiling gradients. + """ + if not 0.5 < partial_fourier_factor <= 1.0: + raise ValueError('Desired partial Fourier factor must be larger than 0.5 and smaller or equal to 1.0.') + + # set system to default if not provided + if system is None: + system = sys_defaults + + self.system = system + self.fov = fov + self.n_readout = n_readout + self.n_phase_encoding = n_phase_encoding + self.bandwidth = bandwidth + self.oversampling = oversampling + self.ramp_sampling = ramp_sampling + self.readout_type = readout_type + self.partial_fourier_factor = partial_fourier_factor + self.adc_freq_offset = adc_freq_offset + self.pe_enable = pe_enable + self.spoiling_enable = spoiling_enable + + # Derived parameters + readout_time = n_readout / bandwidth + delta_kx = 1 / (fov * oversampling) + delta_ky = 1 / fov + + # Initiate all optional gradients as None + self.gx_flyback = None + self.gy_blipdown = None + self.gy_blipup = None + self.gy_blipdownup = None + self.gz_spoil = None + + # Create blip gradient with shortest possible timing + gy_blip_duration = np.ceil(2 * np.sqrt(delta_ky / system.max_slew) / 10e-6 / 2) * 10e-6 * 2 + gy_blip_half_dur = gy_blip_duration / 2 + self.gy_blip = pp.make_trapezoid(channel='y', system=self.system, area=-delta_ky, duration=gy_blip_duration) + + # Create readout gradient + gx_encoding_area = n_readout * delta_kx * oversampling + if self.ramp_sampling: + # Calculate additional gradient area from gy_blip assuming maximum slew rate + extra_area = np.power(gy_blip_half_dur, 2) * self.system.max_slew + + # Create gradient with additional area + gx = pp.make_trapezoid( + channel='x', + system=self.system, + area=gx_encoding_area + extra_area, + duration=gy_blip_half_dur + readout_time + gy_blip_half_dur, + ) + + if not gx.fall_time == gx.rise_time: + raise ValueError('Gradient fall time must be equal to rise time for ramp sampling.') + + # Second, correct area taking actual slew rate into account + gx_slew = gx.amplitude / gx.rise_time + gx_area_reduced_slew = gx.area - gx_slew * np.power(gy_blip_half_dur, 2) + + gx.amplitude = float(gx.amplitude * gx_encoding_area / gx_area_reduced_slew) + gx.area = float(gx.amplitude * (gx.rise_time / 2 + gx.flat_time + gx.fall_time / 2)) + gx.flat_area = float(gx.amplitude * gx.flat_time) + self.gx = gx + else: + self.gx = pp.make_trapezoid( + channel='x', + system=self.system, + flat_area=gx_encoding_area, + flat_time=readout_time, + ) + + # Create ADC event + adc_dwell = delta_kx / self.gx.amplitude + adc_dwell = round_to_raster(adc_dwell, self.system.adc_raster_time) + + adc_samples = int(round(readout_time / adc_dwell / 4) * 4) + + adc_time_to_center = adc_dwell * (adc_samples / 2 + 0.5) + adc_delay = self.gx.rise_time + self.gx.flat_time / 2 - adc_time_to_center + adc_dwell / 2 + # the adc delay has to be rounded to the rf raster time (not the adc raster time) + adc_delay = round_to_raster(adc_delay, self.system.rf_raster_time) + + self.adc = pp.make_adc( + num_samples=adc_samples, + dwell=adc_dwell, + freq_offset=adc_freq_offset, + delay=adc_delay, + system=self.system, + ) + + # Create and align pre-phaser gradients considering partial fourier factor + # determine the number of "PE" lines after (and including) k-space center (independent of partial fourier) + self.n_phase_enc_post_center = int(np.ceil(self.n_phase_encoding / 2 + 1)) + # find the closest number of lines to the desired partial fourier factor + valid_n_phase_total = int(np.ceil(partial_fourier_factor * self.n_phase_encoding)) + # ensure that at least one line before the center is acquired + self.n_phase_enc_pre_center = max(1, valid_n_phase_total - self.n_phase_enc_post_center) + # update the total number of "PE" lines + self.n_phase_enc_total = self.n_phase_enc_pre_center + self.n_phase_enc_post_center + # recalculate the actual partial fourier factor + actual_pf_factor = self.n_phase_enc_total / self.n_phase_encoding + if actual_pf_factor != partial_fourier_factor: + print(f'Adjusted partial Fourier factor from {partial_fourier_factor} to {actual_pf_factor:.2f}.') + self.partial_fourier_factor = actual_pf_factor + + # Create pre-phaser gradients + self.gx_pre = pp.make_trapezoid( + channel='x', + system=self.system, + area=-self.gx.area / 2 - delta_kx / 2, + ) + self.gy_pre = pp.make_trapezoid( + channel='y', + system=self.system, + area=self.n_phase_enc_pre_center * delta_ky, + ) + self.gx_pre, self.gy_pre = pp.align(right=[self.gx_pre, self.gy_pre]) + + # Create and align "phase encoding" gradients + if self.readout_type == 'flyback': + self.gx_flyback = pp.make_trapezoid(channel='x', system=self.system, area=-self.gx.area) + elif self.readout_type == 'symmetric': + # Split and align blip gradient in case of symmetric readout + gy_parts = pp.split_gradient_at(grad=self.gy_blip, time_point=gy_blip_duration / 2, system=self.system) + self.gy_blipup, self.gy_blipdown, _ = pp.align(right=gy_parts[0], left=[gy_parts[1], self.gx]) + self.gy_blipdownup = pp.add_gradients((self.gy_blipdown, self.gy_blipup), system=self.system) + else: + raise NotImplementedError('Currently, only "symmetric" and "flyback" readout types are supported.') + + # Disable phase encoding if self.pe_enable is False + if not self.pe_enable: + self.gy_pre.amplitude = 0 + if ( + self.readout_type == 'symmetric' + and self.gy_blipup is not None + and self.gy_blipdown is not None + and self.gy_blipdownup is not None + ): + self.gy_blipup.waveform *= 0 + self.gy_blipdown.waveform *= 0 + self.gy_blipdownup.waveform *= 0 + elif self.readout_type == 'flyback': + self.gy_blip.waveform *= 0 + + # Create spoiler gradient if spoiling is enabled + if self.spoiling_enable: + self.gz_spoil = pp.make_trapezoid(channel='z', system=self.system, area=4 * delta_ky * n_readout) + + @property + def time_to_center(self) -> float: + """Return time from beginning of readout to center of k-space (needed for TE calculations).""" + # i) add time for pre phaser + time_to_center = pp.calc_duration(self.gx_pre, self.gy_pre) + # ii) add time for completed k-space lines before central (ky = 0) line + if self.readout_type == 'flyback': + time_to_center += self.n_phase_enc_pre_center * ( + pp.calc_duration(self.gx) + pp.calc_duration(self.gx_flyback, self.gy_blip) + ) + elif self.readout_type == 'symmetric': + time_to_center += self.n_phase_enc_pre_center * pp.calc_duration(self.gx) + # iii) add time before start of ADC of central k-space line + if self.ramp_sampling: + time_to_center += self.adc.delay + else: + time_to_center += self.gx.rise_time + # iv) add time from start of ADC (for ky = 0) to timepoint when kx = 0 as well + time_to_center += self.adc.dwell * (self.adc.num_samples / 2 + 0.5) + + return float(time_to_center) + + @property + def time_to_center_without_prephaser(self) -> float: + """Return time from after pre-phasers to center of k-space (needed for TE calculations).""" + return self.time_to_center - pp.calc_duration(self.gx_pre, self.gy_pre) + + @property + def total_duration(self) -> float: + """Return total duration of readout including pre-phaser and optional spoiler.""" + total_duration = pp.calc_duration(self.gx_pre, self.gy_pre) + if self.readout_type == 'flyback': + total_duration += self.n_phase_enc_total * ( + pp.calc_duration(self.gx) + pp.calc_duration(self.gx_flyback, self.gy_blip) + ) + total_duration -= pp.calc_duration(self.gx_flyback, self.gy_blip) + elif self.readout_type == 'symmetric': + total_duration += self.n_phase_enc_total * pp.calc_duration(self.gx) + # add time for spoiler if enabled + if self.spoiling_enable: + total_duration += pp.calc_duration(self.gz_spoil) + + return float(total_duration) + + @property + def total_duration_without_prephaser(self) -> float: + """Return total duration of readout excluding pre-phasers but including optional spoiler.""" + return self.total_duration - pp.calc_duration(self.gx_pre, self.gy_pre) + + def add_to_seq( + self, + seq: pp.Sequence, + add_prephaser: bool = True, + mrd_dataset: ismrmrd.Dataset | None = None, + ) -> tuple[pp.Sequence, ismrmrd.Dataset | None]: + """Add EPI readout blocks to the sequence.""" + # (Re)set phase encoding (LIN) label + lin_label = pp.make_label(label='LIN', type='SET', value=0) + + # Add pre-phaser gradients if enabled + if add_prephaser: + seq.add_block(self.gx_pre, self.gy_pre) + + for pe_idx in range(self.n_phase_enc_total): + rev_label = pp.make_label(type='SET', label='REV', value=self.gx.amplitude < 0) + seg_label = pp.make_label(type='SET', label='SEG', value=self.gx.amplitude < 0) + if self.readout_type == 'symmetric': + # Select blip gradient based on phase encoding index + if pe_idx == 0: + gy_blip = self.gy_blipup + elif pe_idx == self.n_phase_enc_total - 1: + gy_blip = self.gy_blipdown + else: + gy_blip = self.gy_blipdownup + + # Add readout block and reverse polarity of readout gradient + seq.add_block(self.gx, gy_blip, self.adc, lin_label, rev_label, seg_label) + self.gx.amplitude = -self.gx.amplitude + + elif self.readout_type == 'flyback': + seq.add_block(self.gx, self.adc, lin_label, rev_label) + if pe_idx != self.n_phase_enc_total - 1: + seq.add_block(self.gx_flyback, self.gy_blip) + + lin_label = pp.make_label(label='LIN', type='INC', value=1) + + if self.spoiling_enable: + seq.add_block(self.gz_spoil) + + return seq, mrd_dataset + + def plot_sequence(self): + """Plot the sequence.""" + seq = pp.Sequence(self.system) + seq, _ = self.add_to_seq(seq) + _, axs1, _, axs2 = seq.plot(grad_disp='mT/m', plot_now=False) + # add vertical line at time to center to all plots in both figures + for ax in axs1 + axs2: + ax.axvline(x=self.time_to_center, color='r', linestyle='--') + plt.show() + + def plot_trajectory(self): + """Plot k-space trajectory.""" + seq = pp.Sequence(self.system) + # add dummy excitation pulse for trajectory calculation + seq.add_block( + pp.make_block_pulse( + flip_angle=np.pi / 2, + duration=2e-3, + delay=self.system.rf_dead_time, + use='excitation', + system=self.system, + ) + ) + seq, _ = self.add_to_seq(seq) + + # calculate k-space trajectory + k_traj_adc, k_traj, _, _, _ = seq.calculate_kspace() + + # plot trajectory + fig = plt.figure() + plt.plot(k_traj[0], k_traj[1], 'b') + plt.plot(k_traj_adc[0], k_traj_adc[1], 'x', color='red', markersize=4) + plt.grid() + plt.show() + + return fig diff --git a/src/mrseq/utils/__init__.py b/src/mrseq/utils/__init__.py index 00d413e..740fa48 100644 --- a/src/mrseq/utils/__init__.py +++ b/src/mrseq/utils/__init__.py @@ -3,3 +3,4 @@ from mrseq.utils.vds import variable_density_spiral_trajectory from mrseq.utils.trajectory import spiral_acquisition from mrseq.utils.ismrmrd import create_header, combine_ismrmrd_files +from mrseq.utils.EpiReadout import EpiReadout diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py index e4ab95a..91beacd 100644 --- a/src/mrseq/utils/trajectory.py +++ b/src/mrseq/utils/trajectory.py @@ -4,8 +4,6 @@ from typing import Any from typing import Literal -import ismrmrd -import matplotlib.pyplot as plt import numpy as np import pypulseq as pp @@ -530,360 +528,3 @@ def combine_gradients(*grad_objects, channel): time_to_echo = max_pre_duration + n_samples_to_echo * readout_oversampling * adc.dwell return gx_combined, gy_combined, adc, trajectory, time_to_echo - - -class EpiReadout: - """EPI readout module supporting flyback and symmetric readout, ramp sampling, oversampling, and partial Fourier. - - Attributes - ---------- - system - System limits. - ramp_sampling - If True, ADC is active during gradient ramps (optimized timing). - readout_type - Readout type ('flyback' or 'symmetric'). - spoiling_enable - Enable spoiling gradients (useful for calibration if False). - adc - ADC event. - gx_pre - Pre-phaser gradient in readout direction. - gy_pre - Pre-phaser gradient in phase encoding direction. - gx - Readout gradient. - gy - Phase encoding gradient. - gy_blip - Complete blip gradient used for flyback readout. - gy_blipup - Ramp-up part of blip gradient used for symmetric readout. - gy_blipdown - Ramp-down part of blip gradient used for symmetric readout. - gy_blipdownup - Combined gy_blipdown and gy_blipup gradient used for symmetric readout. - gz_spoil - Spoiler gradient in z-direction. - gx_flyback - Flyback gradient in x-direction. - """ - - def __init__( - self, - system: pp.Opts | None = None, - fov: float = 128e-3, - n_readout: int = 64, - n_phase_encoding: int = 64, - bandwidth: float = 50e3, - oversampling: int = 2, - ramp_sampling: bool = True, - readout_type: Literal['flyback', 'symmetric'] = 'symmetric', - partial_fourier_factor: float = 1.0, - adc_freq_offset: float = 0, - pe_enable: bool = True, - spoiling_enable: bool = True, - ): - """Initialize EPI Readout. - - Parameters - ---------- - system - System limits. - fov - Field of view in meters (square). - n_readout - Number of readout points. - n_phase_encoding - Number of phase encoding points. - bandwidth - Total receiver bandwidth in Hz. - oversampling - ADC oversampling factor. - ramp_sampling - If True, ADC is active during gradient ramps (optimized timing). - readout_type - Readout type ('flyback' or 'symmetric'). - partial_fourier_factor - Partial Fourier factor (0.5 to 1.0). - adc_freq_offset - Frequency offset for the ADC. - pe_enable - Enable phase encoding (useful for calibration scans if False). - spoiling_enable - Enable spoiling gradients. - """ - if not 0.5 < partial_fourier_factor <= 1.0: - raise ValueError('Desired partial Fourier factor must be larger than 0.5 and smaller or equal to 1.0.') - - # set system to default if not provided - if system is None: - system = sys_defaults - - self.system = system - self.fov = fov - self.n_readout = n_readout - self.n_phase_encoding = n_phase_encoding - self.bandwidth = bandwidth - self.oversampling = oversampling - self.ramp_sampling = ramp_sampling - self.readout_type = readout_type - self.partial_fourier_factor = partial_fourier_factor - self.adc_freq_offset = adc_freq_offset - self.pe_enable = pe_enable - self.spoiling_enable = spoiling_enable - - # Derived parameters - readout_time = n_readout / bandwidth - delta_kx = 1 / (fov * oversampling) - delta_ky = 1 / fov - - # Initiate all optional gradients as None - self.gx_flyback = None - self.gy_blipdown = None - self.gy_blipup = None - self.gy_blipdownup = None - self.gz_spoil = None - - # Create blip gradient with shortest possible timing - gy_blip_duration = np.ceil(2 * np.sqrt(delta_ky / system.max_slew) / 10e-6 / 2) * 10e-6 * 2 - gy_blip_half_dur = gy_blip_duration / 2 - self.gy_blip = pp.make_trapezoid(channel='y', system=self.system, area=-delta_ky, duration=gy_blip_duration) - - # Create readout gradient - gx_encoding_area = n_readout * delta_kx * oversampling - if self.ramp_sampling: - # Calculate additional gradient area from gy_blip assuming maximum slew rate - extra_area = np.power(gy_blip_half_dur, 2) * self.system.max_slew - - # Create gradient with additional area - gx = pp.make_trapezoid( - channel='x', - system=self.system, - area=gx_encoding_area + extra_area, - duration=gy_blip_half_dur + readout_time + gy_blip_half_dur, - ) - - if not gx.fall_time == gx.rise_time: - raise ValueError('Gradient fall time must be equal to rise time for ramp sampling.') - - # Second, correct area taking actual slew rate into account - gx_slew = gx.amplitude / gx.rise_time - gx_area_reduced_slew = gx.area - gx_slew * np.power(gy_blip_half_dur, 2) - - gx.amplitude = float(gx.amplitude * gx_encoding_area / gx_area_reduced_slew) - gx.area = float(gx.amplitude * (gx.rise_time / 2 + gx.flat_time + gx.fall_time / 2)) - gx.flat_area = float(gx.amplitude * gx.flat_time) - self.gx = gx - else: - self.gx = pp.make_trapezoid( - channel='x', - system=self.system, - flat_area=gx_encoding_area, - flat_time=readout_time, - ) - - # Create ADC event - adc_dwell = delta_kx / self.gx.amplitude - adc_dwell = round_to_raster(adc_dwell, self.system.adc_raster_time) - - adc_samples = int(round(readout_time / adc_dwell / 4) * 4) - - adc_time_to_center = adc_dwell * (adc_samples / 2 + 0.5) - adc_delay = self.gx.rise_time + self.gx.flat_time / 2 - adc_time_to_center + adc_dwell / 2 - # the adc delay has to be rounded to the rf raster time (not the adc raster time) - adc_delay = round_to_raster(adc_delay, self.system.rf_raster_time) - - self.adc = pp.make_adc( - num_samples=adc_samples, - dwell=adc_dwell, - freq_offset=adc_freq_offset, - delay=adc_delay, - system=self.system, - ) - - # Create and align pre-phaser gradients considering partial fourier factor - # determine the number of "PE" lines after (and including) k-space center (independent of partial fourier) - self.n_phase_enc_post_center = int(np.ceil(self.n_phase_encoding / 2 + 1)) - # find the closest number of lines to the desired partial fourier factor - valid_n_phase_total = int(np.ceil(partial_fourier_factor * self.n_phase_encoding)) - # ensure that at least one line before the center is acquired - self.n_phase_enc_pre_center = max(1, valid_n_phase_total - self.n_phase_enc_post_center) - # update the total number of "PE" lines - self.n_phase_enc_total = self.n_phase_enc_pre_center + self.n_phase_enc_post_center - # recalculate the actual partial fourier factor - actual_pf_factor = self.n_phase_enc_total / self.n_phase_encoding - if actual_pf_factor != partial_fourier_factor: - print(f'Adjusted partial Fourier factor from {partial_fourier_factor} to {actual_pf_factor:.2f}.') - self.partial_fourier_factor = actual_pf_factor - - # Create pre-phaser gradients - self.gx_pre = pp.make_trapezoid( - channel='x', - system=self.system, - area=-self.gx.area / 2 - delta_kx / 2, - ) - self.gy_pre = pp.make_trapezoid( - channel='y', - system=self.system, - area=self.n_phase_enc_pre_center * delta_ky, - ) - self.gx_pre, self.gy_pre = pp.align(right=[self.gx_pre, self.gy_pre]) - - # Create and align "phase encoding" gradients - if self.readout_type == 'flyback': - self.gx_flyback = pp.make_trapezoid(channel='x', system=self.system, area=-self.gx.area) - elif self.readout_type == 'symmetric': - # Split and align blip gradient in case of symmetric readout - gy_parts = pp.split_gradient_at(grad=self.gy_blip, time_point=gy_blip_duration / 2, system=self.system) - self.gy_blipup, self.gy_blipdown, _ = pp.align(right=gy_parts[0], left=[gy_parts[1], self.gx]) - self.gy_blipdownup = pp.add_gradients((self.gy_blipdown, self.gy_blipup), system=self.system) - else: - raise NotImplementedError('Currently, only "symmetric" and "flyback" readout types are supported.') - - # Disable phase encoding if self.pe_enable is False - if not self.pe_enable: - self.gy_pre.amplitude = 0 - if ( - self.readout_type == 'symmetric' - and self.gy_blipup is not None - and self.gy_blipdown is not None - and self.gy_blipdownup is not None - ): - self.gy_blipup.waveform *= 0 - self.gy_blipdown.waveform *= 0 - self.gy_blipdownup.waveform *= 0 - elif self.readout_type == 'flyback': - self.gy_blip.waveform *= 0 - - # Create spoiler gradient if spoiling is enabled - if self.spoiling_enable: - self.gz_spoil = pp.make_trapezoid(channel='z', system=self.system, area=4 * delta_ky * n_readout) - - @property - def time_to_center(self) -> float: - """Return time from beginning of readout to center of k-space (needed for TE calculations).""" - # i) add time for pre phaser - time_to_center = pp.calc_duration(self.gx_pre, self.gy_pre) - # ii) add time for completed k-space lines before central (ky = 0) line - if self.readout_type == 'flyback': - time_to_center += self.n_phase_enc_pre_center * ( - pp.calc_duration(self.gx) + pp.calc_duration(self.gx_flyback, self.gy_blip) - ) - elif self.readout_type == 'symmetric': - time_to_center += self.n_phase_enc_pre_center * pp.calc_duration(self.gx) - # iii) add time before start of ADC of central k-space line - if self.ramp_sampling: - time_to_center += self.adc.delay - else: - time_to_center += self.gx.rise_time - # iv) add time from start of ADC (for ky = 0) to timepoint when kx = 0 as well - time_to_center += self.adc.dwell * (self.adc.num_samples / 2 + 0.5) - - return float(time_to_center) - - @property - def time_to_center_without_prephaser(self) -> float: - """Return time from after pre-phasers to center of k-space (needed for TE calculations).""" - return self.time_to_center - pp.calc_duration(self.gx_pre, self.gy_pre) - - @property - def total_duration(self) -> float: - """Return total duration of readout including pre-phaser and optional spoiler.""" - total_duration = pp.calc_duration(self.gx_pre, self.gy_pre) - if self.readout_type == 'flyback': - total_duration += self.n_phase_enc_total * ( - pp.calc_duration(self.gx) + pp.calc_duration(self.gx_flyback, self.gy_blip) - ) - total_duration -= pp.calc_duration(self.gx_flyback, self.gy_blip) - elif self.readout_type == 'symmetric': - total_duration += self.n_phase_enc_total * pp.calc_duration(self.gx) - # add time for spoiler if enabled - if self.spoiling_enable: - total_duration += pp.calc_duration(self.gz_spoil) - - return float(total_duration) - - @property - def total_duration_without_prephaser(self) -> float: - """Return total duration of readout excluding pre-phasers but including optional spoiler.""" - return self.total_duration - pp.calc_duration(self.gx_pre, self.gy_pre) - - def add_to_seq( - self, - seq: pp.Sequence, - add_prephaser: bool = True, - mrd_dataset: ismrmrd.Dataset | None = None, - ) -> tuple[pp.Sequence, ismrmrd.Dataset | None]: - """Add EPI readout blocks to the sequence.""" - # (Re)set phase encoding (LIN) label - lin_label = pp.make_label(label='LIN', type='SET', value=0) - - # Add pre-phaser gradients if enabled - if add_prephaser: - seq.add_block(self.gx_pre, self.gy_pre) - - for pe_idx in range(self.n_phase_enc_total): - rev_label = pp.make_label(type='SET', label='REV', value=self.gx.amplitude < 0) - seg_label = pp.make_label(type='SET', label='SEG', value=self.gx.amplitude < 0) - if self.readout_type == 'symmetric': - # Select blip gradient based on phase encoding index - if pe_idx == 0: - gy_blip = self.gy_blipup - elif pe_idx == self.n_phase_enc_total - 1: - gy_blip = self.gy_blipdown - else: - gy_blip = self.gy_blipdownup - - # Add readout block and reverse polarity of readout gradient - seq.add_block(self.gx, gy_blip, self.adc, lin_label, rev_label, seg_label) - self.gx.amplitude = -self.gx.amplitude - - elif self.readout_type == 'flyback': - seq.add_block(self.gx, self.adc, lin_label, rev_label) - if pe_idx != self.n_phase_enc_total - 1: - seq.add_block(self.gx_flyback, self.gy_blip) - - lin_label = pp.make_label(label='LIN', type='INC', value=1) - - if self.spoiling_enable: - seq.add_block(self.gz_spoil) - - return seq, mrd_dataset - - def plot_sequence(self): - """Plot the sequence.""" - seq = pp.Sequence(self.system) - seq, _ = self.add_to_seq(seq) - _, axs1, _, axs2 = seq.plot(grad_disp='mT/m', plot_now=False) - # add vertical line at time to center to all plots in both figures - for ax in axs1 + axs2: - ax.axvline(x=self.time_to_center, color='r', linestyle='--') - plt.show() - - def plot_trajectory(self): - """Plot k-space trajectory.""" - seq = pp.Sequence(self.system) - # add dummy excitation pulse for trajectory calculation - seq.add_block( - pp.make_block_pulse( - flip_angle=np.pi / 2, - duration=2e-3, - delay=self.system.rf_dead_time, - use='excitation', - system=self.system, - ) - ) - seq, _ = self.add_to_seq(seq) - - # calculate k-space trajectory - k_traj_adc, k_traj, _, _, _ = seq.calculate_kspace() - - # plot trajectory - fig = plt.figure() - plt.plot(k_traj[0], k_traj[1], 'b') - plt.plot(k_traj_adc[0], k_traj_adc[1], 'x', color='red', markersize=4) - plt.grid() - plt.show() - - return fig diff --git a/tests/utils/gradient_waveforms.py b/tests/utils/gradient_waveforms.py new file mode 100644 index 0000000..10ed2a8 --- /dev/null +++ b/tests/utils/gradient_waveforms.py @@ -0,0 +1,42 @@ +"""Helper function to get interpolate gradient waveforms.""" + +import numpy as np +import pypulseq as pp + + +def get_interp_waveform_for_gx_gy(seq: pp.Sequence, dt: np.ndarray | None = None, scale: float = 1.0): + """Interpolate gradient waveforms for the x and y axes. + + Parameters + ---------- + seq + The PyPulseq sequence object containing gradient waveforms. + dt + Desired time points for interpolation. If None, a default time array is generated. + scale + Scaling factor for the gradient waveforms. Default is 1. + + Returns + ------- + gx_waveform_intp + Interpolated gradient waveform for the x-axis. + gy_waveform_intp + Interpolated gradient waveform for the y-axis. + dt + Time points corresponding to the interpolated waveforms. + """ + w = seq.waveforms_and_times() + gx_waveform = w[0][0][1] * scale + gx_waveform_time = w[0][0][0] + + gy_waveform = w[0][1][1] * scale + gy_waveform_time = w[0][1][0] + + if dt is None: + dt = np.arange( + min(gx_waveform_time[0], gy_waveform_time[0]), max(gx_waveform_time[-1], gy_waveform_time[-1]), step=1e-7 + ) + gx_waveform_intp = np.interp(dt, gx_waveform_time, gx_waveform) + gy_waveform_intp = np.interp(dt, gy_waveform_time, gy_waveform) + + return gx_waveform_intp, gy_waveform_intp, dt diff --git a/tests/utils/test_epi_readout.py b/tests/utils/test_epi_readout.py new file mode 100644 index 0000000..98bea7a --- /dev/null +++ b/tests/utils/test_epi_readout.py @@ -0,0 +1,67 @@ +"""Tests for EPI readout.""" + +from typing import Literal + +import numpy as np +import pypulseq as pp +import pytest +from mrseq.utils.EpiReadout import EpiReadout + + +@pytest.mark.parametrize('fov', (64e-3, 128e-3)) +@pytest.mark.parametrize('n_readout', (32, 64)) +@pytest.mark.parametrize('n_phase_encoding', (32, 64)) +@pytest.mark.parametrize('bandwidth', (50e3, 40e3)) +@pytest.mark.parametrize('oversampling', (1, 2)) +@pytest.mark.parametrize('ramp_sampling', (True, False)) +@pytest.mark.parametrize('readout_type', ('flyback', 'symmetric')) +@pytest.mark.parametrize('partial_fourier_factor', (0.75, 1.0)) +@pytest.mark.parametrize('spoiling_enable', (True, False)) +def test_epi_time_to_center( + system_defaults: pp.Opts | None, + fov: float, + n_readout, + n_phase_encoding, + bandwidth: float, + oversampling: int, + ramp_sampling: bool, + readout_type: Literal['flyback', 'symmetric'], + partial_fourier_factor: float, + spoiling_enable: bool, +): + """Test epi readout for different parameter combinations.""" + + def get_time_to_k0(seq): + """Find time when gradients are 0 (k-space center).""" + k_traj_adc, _, _, _, t_adc = seq.calculate_kspace() + m0 = np.sqrt(k_traj_adc[0] ** 2 + k_traj_adc[1] ** 2) + k0_idx = np.argmin(m0[100:-100]) + 100 + return t_adc[k0_idx], np.median(np.diff(t_adc)) + + # create EpiReadout object + epi2d = EpiReadout( + system=system_defaults, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=bandwidth, + oversampling=oversampling, + readout_type=readout_type, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, + adc_freq_offset=0.0, + pe_enable=True, + spoiling_enable=spoiling_enable, + ) + + # Complete EPI readout + seq = pp.Sequence(system=system_defaults) + seq, _ = epi2d.add_to_seq(seq) + assert sum(seq.block_durations.values()) == pytest.approx(epi2d.total_duration, abs=1e-6) + time_to_k0, dwell_time = get_time_to_k0(seq) + assert np.isclose(time_to_k0, epi2d.time_to_center, atol=dwell_time) + + # EPI readout without prephaseser gradients + seq = pp.Sequence(system=system_defaults) + seq, _ = epi2d.add_to_seq(seq, add_prephaser=False) + assert sum(seq.block_durations.values()) == pytest.approx(epi2d.total_duration_without_prephaser, abs=1e-6) diff --git a/tests/utils/test_trajectory.py b/tests/utils/test_trajectory.py index 6a3cfdf..1fe97d8 100644 --- a/tests/utils/test_trajectory.py +++ b/tests/utils/test_trajectory.py @@ -11,6 +11,8 @@ from mrseq.utils.trajectory import undersampled_variable_density_spiral from scipy.signal import argrelextrema +from .gradient_waveforms import get_interp_waveform_for_gx_gy + @pytest.mark.parametrize('n_phase_encoding', [50, 51, 100]) @pytest.mark.parametrize('acceleration', [1, 2, 3, 4, 6]) @@ -87,44 +89,6 @@ def test_multi_gradient_echo(system_defaults): seq, _ = mecho.add_to_seq(seq, n_echoes=3) -def get_interp_waveform_for_gx_gy(seq: pp.Sequence, dt: np.ndarray | None = None, scale: float = 1.0): - """Interpolate gradient waveforms for the x and y axes. - - Parameters - ---------- - seq - The PyPulseq sequence object containing gradient waveforms. - dt - Desired time points for interpolation. If None, a default time array is generated. - scale - Scaling factor for the gradient waveforms. Default is 1. - - Returns - ------- - gx_waveform_intp - Interpolated gradient waveform for the x-axis. - gy_waveform_intp - Interpolated gradient waveform for the y-axis. - dt - Time points corresponding to the interpolated waveforms. - """ - w = seq.waveforms_and_times() - gx_waveform = w[0][0][1] * scale - gx_waveform_time = w[0][0][0] - - gy_waveform = w[0][1][1] * scale - gy_waveform_time = w[0][1][0] - - if dt is None: - dt = np.arange( - min(gx_waveform_time[0], gy_waveform_time[0]), max(gx_waveform_time[-1], gy_waveform_time[-1]), step=1e-7 - ) - gx_waveform_intp = np.interp(dt, gx_waveform_time, gx_waveform) - gy_waveform_intp = np.interp(dt, gy_waveform_time, gy_waveform) - - return gx_waveform_intp, gy_waveform_intp, dt - - @pytest.mark.parametrize('n_readout', (128, 256)) @pytest.mark.parametrize('fov', (128e-3, 320e-3)) @pytest.mark.parametrize('undersampling_factor', (1, 2, 3, 4)) From 2543378cb79306b42ad6182efe2cecf0676467c1 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 3 Mar 2026 11:30:51 +0100 Subject: [PATCH 13/32] replace traj calculation for mrd file with analytical expression --- src/mrseq/scripts/epi2d_fid.py | 73 ++++++----- src/mrseq/scripts/epi2d_se.py | 68 ++++++----- src/mrseq/utils/EpiReadout.py | 214 ++++++++++++++++++++++++++++++--- 3 files changed, 278 insertions(+), 77 deletions(-) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index 042d0a5..72e9d23 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -106,9 +106,6 @@ def epi2d_fid_kernel( # define number of navigator acquisitions n_navigator_acq = 3 if add_navigator_acq else 0 - # define number of noise acquisitions - n_noise_acq = 1 if add_noise_acq else 0 - # create EpiReadout object epi2d = EpiReadout( system=system, @@ -174,6 +171,7 @@ def epi2d_fid_kernel( raise ValueError(f'TR must be larger than {min_tr * 1000:.2f} ms. Current value is {tr * 1000:.3f} ms.') # create header + prot = None if mrd_header_file: hdr = create_header( traj_type='other', @@ -191,6 +189,15 @@ def epi2d_fid_kernel( prot = ismrmrd.Dataset(mrd_header_file, 'w') prot.write_xml_header(hdr.toXML('utf-8')) + # Precompute analytical navigator trajectory (single kx line, no ky blip) + if add_navigator_acq: + from mrseq.utils.EpiReadout import _trapezoid_area_at_times + + nav_sample_times = epi2d.adc.delay + (np.arange(epi2d.adc.num_samples) + 0.5) * epi2d.adc.dwell + nav_kx_forward = _trapezoid_area_at_times( + epi2d.gx.rise_time, epi2d.gx.flat_time, epi2d.gx.fall_time, abs(epi2d.gx.amplitude), nav_sample_times + ) + # obtain noise samples if selected if add_noise_acq: seq.add_block( @@ -204,6 +211,14 @@ def epi2d_fid_kernel( seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) seq.add_block(pp.make_delay(system.rf_dead_time)) + # Write noise trajectory to MRD (zero trajectory — no gradients active) + if mrd_header_file: + n_samples = epi2d.adc.num_samples + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) + acq.traj[:] = np.zeros((n_samples, 2), dtype=np.float32) + prot.append_acquisition(acq) + for slice_ in range(n_slices): # define slice label slice_label = pp.make_label(label='SLC', type='SET', value=slice_) @@ -231,6 +246,10 @@ def epi2d_fid_kernel( # reverse gx_pre back after adding to sequence gx_pre = pp.scale_grad(gx_pre, -1) + # Navigator kx offset: starts from -gx_pre.area (reversed pre-winder) + nav_kx_offset = -epi2d.gx_pre.area + nav_gx_sign = -1.0 # first navigator uses reversed gx + # add navigator scans for ghost correction for n in range(n_navigator_acq): seq.add_block( @@ -240,7 +259,25 @@ def epi2d_fid_kernel( pp.make_label(label='SEG', type='SET', value=gx.amplitude < 0), pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), ) + + # Write navigator trajectory to MRD + if mrd_header_file: + n_samples = epi2d.adc.num_samples + traj = np.zeros((n_samples, 2), dtype=np.float32) + if nav_gx_sign > 0: + nav_kx = nav_kx_offset + nav_kx_forward + else: + nav_kx = nav_kx_offset - nav_kx_forward + traj[:, 0] = np.round(nav_kx * fov_xy * readout_oversampling, 3) + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) + acq.traj[:] = traj + prot.append_acquisition(acq) + + # Update kx offset for next navigator (accumulate the signed gx area) + nav_kx_offset += nav_gx_sign * epi2d.gx.area gx = pp.scale_grad(gx, -1) + nav_gx_sign = -nav_gx_sign # add echo time delay seq.add_block(pp.make_delay(te_delay)) @@ -260,40 +297,16 @@ def epi2d_fid_kernel( seq.add_block(gzr, gx_pre, gy_pre) # add EPI readout block without pre-phaser gradients + # (trajectory is written per-readout inside add_to_seq when mrd_dataset is provided) seq, prot = epi2d.add_to_seq(seq, add_prephaser=False, mrd_dataset=prot) # add repetition time delay if tr_delay > 0: seq.add_block(pp.make_delay(tr_delay)) + # close ISMRMRD file if mrd_header_file: - # calculate k-space trajectory from Sequence to add this info to mrd file - k_traj_adc, _, _, _, _ = seq.calculate_kspace() - samples_per_acq = epi2d.adc.num_samples - n_total_acq = k_traj_adc.shape[-1] // samples_per_acq - n_epi_acq_per_slice = epi2d.n_phase_enc_total - - if not (n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq == n_total_acq: - raise ValueError('Number of calculated acquisitions does not match expected number.') - - # create mrd acquisition and add trajectory info for each EPI readout including navigator and noise acquisitions - for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq): - start = n * samples_per_acq - end = start + samples_per_acq - - traj = np.zeros((samples_per_acq, 2), dtype=np.float32) - traj[:, 0] = np.round(k_traj_adc[0, start:end] * fov_xy * readout_oversampling, 3) - traj[:, 1] = np.round(k_traj_adc[1, start:end] * fov_xy, 3) - - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=samples_per_acq) - acq.traj[:] = traj - - prot.append_acquisition(acq) - - # close ISMRMRD file - if mrd_header_file: - prot.close() + prot.close() # set gridding definitions extracted from EpiReadout if ramp_sampling: diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 1802436..26a943d 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -106,9 +106,6 @@ def epi2d_se_kernel( # define number of navigator acquisitions n_navigator_acq = 3 if add_navigator_acq else 0 - # define number of noise acquisitions - n_noise_acq = 1 if add_noise_acq else 0 - # create EpiReadout object epi2d = EpiReadout( system=system, @@ -237,6 +234,7 @@ def epi2d_se_kernel( print(f'Current repetition time = {(min_tr + tr_delay) * 1000:.4f} ms') # create header + prot = None if mrd_header_file: hdr = create_header( traj_type='other', @@ -254,6 +252,15 @@ def epi2d_se_kernel( prot = ismrmrd.Dataset(mrd_header_file, 'w') prot.write_xml_header(hdr.toXML('utf-8')) + # Precompute analytical navigator trajectory (single kx line, no ky blip) + if add_navigator_acq: + from mrseq.utils.EpiReadout import _trapezoid_area_at_times + + nav_sample_times = epi2d.adc.delay + (np.arange(epi2d.adc.num_samples) + 0.5) * epi2d.adc.dwell + nav_kx_forward = _trapezoid_area_at_times( + epi2d.gx.rise_time, epi2d.gx.flat_time, epi2d.gx.fall_time, abs(epi2d.gx.amplitude), nav_sample_times + ) + # obtain noise samples if selected if add_noise_acq: seq.add_block( @@ -267,6 +274,14 @@ def epi2d_se_kernel( seq.add_block(pp.make_label(label='NOISE', type='SET', value=False)) seq.add_block(pp.make_delay(system.rf_dead_time)) + # Write noise trajectory to MRD (zero trajectory — no gradients active) + if mrd_header_file: + n_samples = epi2d.adc.num_samples + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) + acq.traj[:] = np.zeros((n_samples, 2), dtype=np.float32) + prot.append_acquisition(acq) + for slice_ in range(n_slices): # define slice label slice_label = pp.make_label(label='SLC', type='SET', value=slice_) @@ -294,6 +309,10 @@ def epi2d_se_kernel( # reverse gx_pre back after adding to sequence gx_pre = pp.scale_grad(gx_pre, -1) + # Navigator kx offset: starts from -gx_pre.area (reversed pre-winder) + nav_kx_offset = -epi2d.gx_pre.area + nav_gx_sign = -1.0 # first navigator uses reversed gx + # add 3 navigator acquisitions for n in range(n_navigator_acq): seq.add_block( @@ -303,7 +322,25 @@ def epi2d_se_kernel( pp.make_label(label='SEG', type='SET', value=gx.amplitude < 0), pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), ) + + # Write navigator trajectory to MRD + if mrd_header_file: + n_samples = epi2d.adc.num_samples + traj = np.zeros((n_samples, 2), dtype=np.float32) + if nav_gx_sign > 0: + nav_kx = nav_kx_offset + nav_kx_forward + else: + nav_kx = nav_kx_offset - nav_kx_forward + traj[:, 0] = np.round(nav_kx * fov_xy * readout_oversampling, 3) + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) + acq.traj[:] = traj + prot.append_acquisition(acq) + + # Update kx offset for next navigator (accumulate the signed gx area) + nav_kx_offset += nav_gx_sign * epi2d.gx.area gx = pp.scale_grad(gx, -1) + nav_gx_sign = -nav_gx_sign # add gy_pre and reset labels seq.add_block( @@ -324,36 +361,13 @@ def epi2d_se_kernel( seq.add_block(pp.make_delay(te_delay2)) # add EPI readout block without pre-phaser gradients + # (trajectory is written per-readout inside add_to_seq when mrd_dataset is provided) seq, prot = epi2d.add_to_seq(seq, add_prephaser=False, mrd_dataset=prot) # add repetition time delay if tr_delay > 0: seq.add_block(pp.make_delay(tr_delay)) - # calculate k-space trajectory from Sequence to add this info to mrd file - k_traj_adc, _, _, _, _ = seq.calculate_kspace() - samples_per_acq = epi2d.adc.num_samples - n_total_acq = k_traj_adc.shape[-1] // samples_per_acq - n_epi_acq_per_slice = epi2d.n_phase_enc_total - - if not (n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq == n_total_acq: - raise ValueError('Number of calculated acquisitions does not match expected number.') - - # create mrd acquisition and add trajectory info for each EPI readout including navigator and noise acquisitions - for n in range((n_epi_acq_per_slice + n_navigator_acq) * n_slices + n_noise_acq): - start = n * samples_per_acq - end = start + samples_per_acq - - traj = np.zeros((samples_per_acq, 2), dtype=np.float32) - traj[:, 0] = np.round(k_traj_adc[0, start:end] * fov_xy * readout_oversampling, 3) - traj[:, 1] = np.round(k_traj_adc[1, start:end] * fov_xy, 3) - - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=samples_per_acq) - acq.traj[:] = traj - - prot.append_acquisition(acq) - # close ISMRMRD file if mrd_header_file: prot.close() diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index e14395c..f820933 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -11,6 +11,99 @@ from mrseq.utils import sys_defaults +def _trapezoid_area_at_times( + rise_time: float, + flat_time: float, + fall_time: float, + amplitude: float, + sample_times: np.ndarray, +) -> np.ndarray: + """Calculate accumulated gradient area analytically for a trapezoidal gradient. + + Parameters + ---------- + rise_time + Gradient rise time (in seconds). + flat_time + Gradient flat time (in seconds). + fall_time + Gradient fall time (in seconds). + amplitude + Gradient amplitude (in Hz/m). + sample_times + Array of time points (in seconds) relative to the start of the gradient. + + Returns + ------- + area + Accumulated gradient area at each sample time (in 1/m). + """ + t1 = rise_time + t2 = rise_time + flat_time + + area_at_t1 = amplitude * rise_time / 2 + area_at_t2 = area_at_t1 + amplitude * flat_time + + area = np.zeros_like(sample_times) + for i, t in enumerate(sample_times): + if t <= 0: + area[i] = 0.0 + elif t <= t1: + area[i] = amplitude * t**2 / (2 * rise_time) + elif t <= t2: + area[i] = area_at_t1 + amplitude * (t - t1) + else: + dt = t - t2 + area[i] = area_at_t2 + amplitude * dt - amplitude * dt**2 / (2 * fall_time) + + return area + + +def _piecewise_linear_area_at_times( + wf_times: np.ndarray, + wf_values: np.ndarray, + sample_times: np.ndarray, +) -> np.ndarray: + """Compute cumulative integral of a piecewise-linear waveform at given sample times. + + Parameters + ---------- + wf_times + Time points of the piecewise-linear waveform (in seconds), sorted ascending. + wf_values + Waveform amplitude values at each time point (in Hz/m). + sample_times + Times at which to evaluate the cumulative area (in seconds). + + Returns + ------- + area + Cumulative area at each sample time (in 1/m). + """ + n_knots = len(wf_times) + cum_area = np.zeros(n_knots) + for k in range(1, n_knots): + dt = wf_times[k] - wf_times[k - 1] + cum_area[k] = cum_area[k - 1] + (wf_values[k] + wf_values[k - 1]) / 2 * dt + + area = np.zeros_like(sample_times) + for i, t in enumerate(sample_times): + if t <= wf_times[0]: + area[i] = 0.0 + elif t >= wf_times[-1]: + area[i] = cum_area[-1] + else: + k = np.searchsorted(wf_times, t, side='right') + base_area = cum_area[k - 1] + dt_seg = wf_times[k] - wf_times[k - 1] + dt = t - wf_times[k - 1] + frac = dt / dt_seg + wf_at_t = wf_values[k - 1] + frac * (wf_values[k] - wf_values[k - 1]) + area[i] = base_area + (wf_values[k - 1] + wf_at_t) / 2 * dt + + return area + + class EpiReadout: """EPI readout module supporting flyback and symmetric readout, ramp sampling, oversampling, and partial Fourier. @@ -288,16 +381,98 @@ def total_duration_without_prephaser(self) -> float: """Return total duration of readout excluding pre-phasers but including optional spoiler.""" return self.total_duration - pp.calc_duration(self.gx_pre, self.gy_pre) + def calculate_trajectory(self) -> tuple[np.ndarray, np.ndarray]: + """Compute analytical k-space trajectory for the EPI readout. + + Returns the kx and ky positions at each ADC sample for each phase encoding line. + This replaces the need to call ``seq.calculate_kspace()`` for trajectory information. + + Returns + ------- + kx + kx trajectory in 1/m, shape ``(n_phase_enc_total, adc.num_samples)``. + ky + ky trajectory in 1/m, shape ``(n_phase_enc_total, adc.num_samples)``. + """ + n_samples = self.adc.num_samples + n_pe = self.n_phase_enc_total + + # ADC sample times relative to the start of the gradient block (center of each dwell) + sample_times = self.adc.delay + (np.arange(n_samples) + 0.5) * self.adc.dwell + + # kx for a single readout line with positive gradient amplitude + kx_forward = _trapezoid_area_at_times( + self.gx.rise_time, self.gx.flat_time, self.gx.fall_time, abs(self.gx.amplitude), sample_times + ) + + # After gx_pre, kx starts at gx_pre.area + kx_offset = self.gx_pre.area + + # ky: pre-phaser moves to top of k-space + ky_start = self.gy_pre.area + + # For symmetric EPI, compute per-sample ky offset from the blip waveforms. + # tt is relative to the waveform start; add delay to get block-relative times. + if self.readout_type == 'symmetric': + bu_tt = self.gy_blipup.tt + self.gy_blipup.delay + bd_tt = self.gy_blipdown.tt + self.gy_blipdown.delay + bdu_tt = self.gy_blipdownup.tt + self.gy_blipdownup.delay + + ky_offset_blipup = _piecewise_linear_area_at_times(bu_tt, self.gy_blipup.waveform, sample_times) + ky_offset_blipdown = _piecewise_linear_area_at_times(bd_tt, self.gy_blipdown.waveform, sample_times) + ky_offset_blipdownup = _piecewise_linear_area_at_times(bdu_tt, self.gy_blipdownup.waveform, sample_times) + + blipup_total = _piecewise_linear_area_at_times(bu_tt, self.gy_blipup.waveform, np.array([bu_tt[-1]]))[0] + blipdownup_total = _piecewise_linear_area_at_times( + bdu_tt, self.gy_blipdownup.waveform, np.array([bdu_tt[-1]]) + )[0] + + kx = np.zeros((n_pe, n_samples)) + ky = np.zeros((n_pe, n_samples)) + + # Track cumulative ky base (ky value at the START of each block) + ky_base = ky_start + + for pe_idx in range(n_pe): + if self.readout_type == 'flyback': + kx[pe_idx, :] = kx_offset + kx_forward + ky[pe_idx, :] = ky_start + pe_idx * self.gy_blip.area + + elif self.readout_type == 'symmetric': + if pe_idx % 2 == 0: + kx[pe_idx, :] = kx_offset + kx_forward + else: + kx[pe_idx, :] = kx_offset + self.gx.area - kx_forward + + if pe_idx == 0: + ky[pe_idx, :] = ky_base + ky_offset_blipup + ky_base += blipup_total + elif pe_idx == n_pe - 1: + ky[pe_idx, :] = ky_base + ky_offset_blipdown + else: + ky[pe_idx, :] = ky_base + ky_offset_blipdownup + ky_base += blipdownup_total + + return kx, ky + def add_to_seq( self, seq: pp.Sequence, add_prephaser: bool = True, mrd_dataset: ismrmrd.Dataset | None = None, ) -> tuple[pp.Sequence, ismrmrd.Dataset | None]: - """Add EPI readout blocks to the sequence.""" + """Add EPI readout blocks to the sequence. + + If `mrd_dataset` is provided, the analytical k-space trajectory is + computed and written to the MRD file for each readout line. + """ # (Re)set phase encoding (LIN) label lin_label = pp.make_label(label='LIN', type='SET', value=0) + # Compute analytical trajectory for MRD writing + if mrd_dataset is not None: + kx_traj, ky_traj = self.calculate_trajectory() + # Add pre-phaser gradients if enabled if add_prephaser: seq.add_block(self.gx_pre, self.gy_pre) @@ -323,6 +498,18 @@ def add_to_seq( if pe_idx != self.n_phase_enc_total - 1: seq.add_block(self.gx_flyback, self.gy_blip) + # Write per-readout trajectory to MRD + if mrd_dataset is not None: + n_samples = self.adc.num_samples + traj = np.zeros((n_samples, 2), dtype=np.float32) + traj[:, 0] = np.round(kx_traj[pe_idx] * self.fov * self.oversampling, 3) + traj[:, 1] = np.round(ky_traj[pe_idx] * self.fov, 3) + + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) + acq.traj[:] = traj + mrd_dataset.append_acquisition(acq) + lin_label = pp.make_label(label='LIN', type='INC', value=1) if self.spoiling_enable: @@ -341,27 +528,14 @@ def plot_sequence(self): plt.show() def plot_trajectory(self): - """Plot k-space trajectory.""" - seq = pp.Sequence(self.system) - # add dummy excitation pulse for trajectory calculation - seq.add_block( - pp.make_block_pulse( - flip_angle=np.pi / 2, - duration=2e-3, - delay=self.system.rf_dead_time, - use='excitation', - system=self.system, - ) - ) - seq, _ = self.add_to_seq(seq) - - # calculate k-space trajectory - k_traj_adc, k_traj, _, _, _ = seq.calculate_kspace() + """Plot k-space trajectory using the analytical expression.""" + kx, ky = self.calculate_trajectory() - # plot trajectory fig = plt.figure() - plt.plot(k_traj[0], k_traj[1], 'b') - plt.plot(k_traj_adc[0], k_traj_adc[1], 'x', color='red', markersize=4) + plt.plot(kx.T, ky.T, 'b') + plt.plot(kx.ravel(), ky.ravel(), 'x', color='red', markersize=4) + plt.xlabel('kx (1/m)') + plt.ylabel('ky (1/m)') plt.grid() plt.show() From e744aeea417f35b25e754a2fe7cc1aaca2da18e1 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 3 Mar 2026 13:12:13 +0100 Subject: [PATCH 14/32] remove rounding of traj values in mrd file --- src/mrseq/scripts/epi2d_fid.py | 2 +- src/mrseq/scripts/epi2d_se.py | 2 +- src/mrseq/utils/EpiReadout.py | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index 72e9d23..12047d0 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -268,7 +268,7 @@ def epi2d_fid_kernel( nav_kx = nav_kx_offset + nav_kx_forward else: nav_kx = nav_kx_offset - nav_kx_forward - traj[:, 0] = np.round(nav_kx * fov_xy * readout_oversampling, 3) + traj[:, 0] = nav_kx * fov_xy * readout_oversampling acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) acq.traj[:] = traj diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 26a943d..27619e0 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -331,7 +331,7 @@ def epi2d_se_kernel( nav_kx = nav_kx_offset + nav_kx_forward else: nav_kx = nav_kx_offset - nav_kx_forward - traj[:, 0] = np.round(nav_kx * fov_xy * readout_oversampling, 3) + traj[:, 0] = nav_kx * fov_xy * readout_oversampling acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) acq.traj[:] = traj diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index f820933..df13af3 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -502,8 +502,8 @@ def add_to_seq( if mrd_dataset is not None: n_samples = self.adc.num_samples traj = np.zeros((n_samples, 2), dtype=np.float32) - traj[:, 0] = np.round(kx_traj[pe_idx] * self.fov * self.oversampling, 3) - traj[:, 1] = np.round(ky_traj[pe_idx] * self.fov, 3) + traj[:, 0] = kx_traj[pe_idx] * self.fov * self.oversampling + traj[:, 1] = ky_traj[pe_idx] * self.fov acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) From 2c1d2d433fc7b7a7bf173a1a6d7a8a7d0c952feb Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 3 Mar 2026 13:12:42 +0100 Subject: [PATCH 15/32] add tests for analytical epi traj --- tests/scripts/test_epi2d_fid.py | 104 ++++++++++++++++++++++++++++ tests/scripts/test_epi2d_se.py | 116 ++++++++++++++++++++++++++++++++ tests/utils/test_epi_readout.py | 57 ++++++++++++++++ 3 files changed, 277 insertions(+) diff --git a/tests/scripts/test_epi2d_fid.py b/tests/scripts/test_epi2d_fid.py index 36953e8..0aba17f 100644 --- a/tests/scripts/test_epi2d_fid.py +++ b/tests/scripts/test_epi2d_fid.py @@ -1,7 +1,15 @@ """Tests for EPI FID sequence.""" +import tempfile +from pathlib import Path +from typing import Literal + +import ismrmrd +import numpy as np import pytest +from mrseq.scripts.epi2d_fid import epi2d_fid_kernel from mrseq.scripts.epi2d_fid import main as create_seq +from mrseq.utils.EpiReadout import EpiReadout EXPECTED_DUR = 0.0389 # defined 2026-02-28 @@ -23,3 +31,99 @@ def test_seq_creation_error_on_short_tr(system_defaults): """Test if error is raised on too short repetition time.""" with pytest.raises(ValueError): create_seq(system=system_defaults, tr=2e-3, show_plots=False) + + +@pytest.mark.parametrize('readout_type', ('flyback', 'symmetric')) +@pytest.mark.parametrize( + 'add_navigator_acq, add_noise_acq', + [(False, False), (True, False), (False, True), (True, True)], +) +def test_mrd_trajectory( + system_defaults, + readout_type: Literal['flyback', 'symmetric'], + add_navigator_acq: bool, + add_noise_acq: bool, +): + """Test that the MRD trajectory matches the analytical and PyPulseq trajectories.""" + fov = 200e-3 + n_readout = 64 + n_phase_encoding = 64 + oversampling = 2 + n_slices = 1 + + mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_fid_traj.h5' + if mrd_file.exists(): + mrd_file.unlink() + + seq, _, _ = epi2d_fid_kernel( + system=system_defaults, + te=None, + tr=None, + fov_xy=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=100e3, + slice_thickness=5e-3, + n_slices=n_slices, + rf_duration=2e-3, + rf_flip_angle=90, + rf_bwt=4, + rf_apodization=0.5, + readout_type=readout_type, + readout_oversampling=oversampling, + ramp_sampling=True, + partial_fourier_factor=1.0, + add_spoiler=True, + add_noise_acq=add_noise_acq, + add_navigator_acq=add_navigator_acq, + mrd_header_file=mrd_file, + ) + + # Read MRD trajectories + ds = ismrmrd.Dataset(mrd_file, 'w', False) + n_acq = ds.number_of_acquisitions() + traj_mrd = [ds.read_acquisition(i).traj.copy() for i in range(n_acq)] + ds.close() + mrd_file.unlink() + + # Recreate EpiReadout for reference + epi2d = EpiReadout( + system=system_defaults, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=100e3, + oversampling=oversampling, + ramp_sampling=True, + readout_type=readout_type, + partial_fourier_factor=1.0, + pe_enable=True, + spoiling_enable=True, + ) + n_samples = epi2d.adc.num_samples + n_nav = 3 if add_navigator_acq else 0 + n_noise = 1 if add_noise_acq else 0 + n_pre_per_slice = n_nav + n_epi = epi2d.n_phase_enc_total + assert n_acq == n_noise + (n_pre_per_slice + n_epi) * n_slices + + # Compare noise + navigator against calculate_kspace + k_traj_adc, _, _, _, _ = seq.calculate_kspace() + n_pre_total = n_noise + n_pre_per_slice * n_slices + for i in range(n_pre_total): + start = i * n_samples + end = start + n_samples + kx_calc = (k_traj_adc[0, start:end] * fov * oversampling).astype(np.float32) + ky_calc = (k_traj_adc[1, start:end] * fov).astype(np.float32) + np.testing.assert_array_equal(traj_mrd[i][:, 0], kx_calc) + np.testing.assert_array_equal(traj_mrd[i][:, 1], ky_calc) + + # Compare EPI readout lines against analytical trajectory + kx_analytical, ky_analytical = epi2d.calculate_trajectory() + for s in range(n_slices): + for pe_idx in range(n_epi): + mrd_idx = n_noise + s * (n_nav + n_epi) + n_nav + pe_idx + kx_expected = (kx_analytical[pe_idx] * fov * oversampling).astype(np.float32) + ky_expected = (ky_analytical[pe_idx] * fov).astype(np.float32) + np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 0], kx_expected) + np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 1], ky_expected) diff --git a/tests/scripts/test_epi2d_se.py b/tests/scripts/test_epi2d_se.py index 9a2000b..a0e3921 100644 --- a/tests/scripts/test_epi2d_se.py +++ b/tests/scripts/test_epi2d_se.py @@ -1,7 +1,15 @@ """Tests for EPI SE sequence.""" +import tempfile +from pathlib import Path +from typing import Literal + +import ismrmrd +import numpy as np import pytest +from mrseq.scripts.epi2d_se import epi2d_se_kernel from mrseq.scripts.epi2d_se import main as create_seq +from mrseq.utils.EpiReadout import EpiReadout EXPECTED_DUR = 0.22586 # defined 2026-02-28 @@ -23,3 +31,111 @@ def test_seq_creation_error_on_short_tr(system_defaults): """Test if error is raised on too short repetition time.""" with pytest.raises(ValueError): create_seq(system=system_defaults, tr=2e-3, show_plots=False) + + +@pytest.mark.parametrize('readout_type', ('flyback', 'symmetric')) +@pytest.mark.parametrize( + 'add_navigator_acq, add_noise_acq', + [(False, False), (True, False), (False, True), (True, True)], +) +def test_mrd_trajectory( + system_defaults, + readout_type: Literal['flyback', 'symmetric'], + add_navigator_acq: bool, + add_noise_acq: bool, +): + """Test that the MRD trajectory matches the analytical and PyPulseq trajectories.""" + fov = 200e-3 + n_readout = 64 + n_phase_encoding = 64 + oversampling = 2 + n_slices = 1 + + mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_se_traj.h5' + if mrd_file.exists(): + mrd_file.unlink() + + seq, _, _ = epi2d_se_kernel( + system=system_defaults, + te=None, + tr=None, + fov_xy=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=100e3, + slice_thickness=5e-3, + n_slices=n_slices, + rf_duration=2e-3, + rf_flip_angle=90, + rf_bwt=4, + rf_apodization=0.5, + readout_type=readout_type, + readout_oversampling=oversampling, + ramp_sampling=True, + partial_fourier_factor=1.0, + add_spoiler=True, + add_noise_acq=add_noise_acq, + add_navigator_acq=add_navigator_acq, + mrd_header_file=mrd_file, + ) + + # Read MRD trajectories + ds = ismrmrd.Dataset(mrd_file, 'w', False) + n_acq = ds.number_of_acquisitions() + traj_mrd = [ds.read_acquisition(i).traj.copy() for i in range(n_acq)] + ds.close() + mrd_file.unlink() + + # Recreate EpiReadout for reference + epi2d = EpiReadout( + system=system_defaults, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=100e3, + oversampling=oversampling, + ramp_sampling=True, + readout_type=readout_type, + partial_fourier_factor=1.0, + pe_enable=True, + spoiling_enable=True, + ) + n_samples = epi2d.adc.num_samples + n_nav = 3 if add_navigator_acq else 0 + n_noise = 1 if add_noise_acq else 0 + n_pre_per_slice = n_nav + n_epi = epi2d.n_phase_enc_total + assert n_acq == n_noise + (n_pre_per_slice + n_epi) * n_slices + + # Compare noise + navigator against calculate_kspace + k_traj_adc, _, _, _, _ = seq.calculate_kspace() + n_pre_total = n_noise + n_pre_per_slice * n_slices + for i in range(n_pre_total): + start = i * n_samples + end = start + n_samples + kx_calc = (k_traj_adc[0, start:end] * fov * oversampling).astype(np.float32) + ky_calc = (k_traj_adc[1, start:end] * fov).astype(np.float32) + np.testing.assert_array_equal(traj_mrd[i][:, 0], kx_calc) + np.testing.assert_array_equal(traj_mrd[i][:, 1], ky_calc) + + # Compare navigator trajectories per slice against calculate_kspace + for s in range(n_slices): + for n in range(n_nav): + mrd_idx = n_noise + s * (n_nav + n_epi) + n + calc_idx = n_noise + s * (n_nav + n_epi) + n + start = calc_idx * n_samples + end = start + n_samples + kx_calc = (k_traj_adc[0, start:end] * fov * oversampling).astype(np.float32) + ky_calc = (k_traj_adc[1, start:end] * fov).astype(np.float32) + np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 0], kx_calc) + np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 1], ky_calc) + + # Compare EPI readout lines against analytical trajectory + kx_analytical, ky_analytical = epi2d.calculate_trajectory() + for s in range(n_slices): + for pe_idx in range(n_epi): + mrd_idx = n_noise + s * (n_nav + n_epi) + n_nav + pe_idx + kx_expected = (kx_analytical[pe_idx] * fov * oversampling).astype(np.float32) + ky_expected = (ky_analytical[pe_idx] * fov).astype(np.float32) + np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 0], kx_expected) + np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 1], ky_expected) diff --git a/tests/utils/test_epi_readout.py b/tests/utils/test_epi_readout.py index 98bea7a..82c29f3 100644 --- a/tests/utils/test_epi_readout.py +++ b/tests/utils/test_epi_readout.py @@ -65,3 +65,60 @@ def get_time_to_k0(seq): seq = pp.Sequence(system=system_defaults) seq, _ = epi2d.add_to_seq(seq, add_prephaser=False) assert sum(seq.block_durations.values()) == pytest.approx(epi2d.total_duration_without_prephaser, abs=1e-6) + + +@pytest.mark.parametrize('fov', (64e-3, 128e-3)) +@pytest.mark.parametrize('n_readout', (32, 64)) +@pytest.mark.parametrize('n_phase_encoding', (32, 64)) +@pytest.mark.parametrize('bandwidth', (50e3, 40e3)) +@pytest.mark.parametrize('oversampling', (1, 2)) +@pytest.mark.parametrize('ramp_sampling', (True, False)) +@pytest.mark.parametrize('readout_type', ('flyback', 'symmetric')) +@pytest.mark.parametrize('partial_fourier_factor', (0.75, 1.0)) +def test_epi_analytical_trajectory( + system_defaults: pp.Opts | None, + fov: float, + n_readout: int, + n_phase_encoding: int, + bandwidth: float, + oversampling: int, + ramp_sampling: bool, + readout_type: Literal['flyback', 'symmetric'], + partial_fourier_factor: float, +): + """Test that the analytical trajectory matches calculate_kspace().""" + epi2d = EpiReadout( + system=system_defaults, + fov=fov, + n_readout=n_readout, + n_phase_encoding=n_phase_encoding, + bandwidth=bandwidth, + oversampling=oversampling, + readout_type=readout_type, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, + adc_freq_offset=0.0, + pe_enable=True, + spoiling_enable=True, + ) + + # Build sequence and get PyPulseq trajectory + seq = pp.Sequence(system=system_defaults) + seq, _ = epi2d.add_to_seq(seq) + k_traj_adc, _, _, _, _ = seq.calculate_kspace() + + # Get analytical trajectory + kx_analytical, ky_analytical = epi2d.calculate_trajectory() + + n_samples = epi2d.adc.num_samples + n_pe = epi2d.n_phase_enc_total + + # The prephaser block has no ADC, so calculate_kspace ADC samples start at the first readout + for pe_idx in range(n_pe): + start = pe_idx * n_samples + end = start + n_samples + kx_pulseq = k_traj_adc[0, start:end] + ky_pulseq = k_traj_adc[1, start:end] + + np.testing.assert_allclose(kx_analytical[pe_idx], kx_pulseq, atol=1e-9) + np.testing.assert_allclose(ky_analytical[pe_idx], ky_pulseq, atol=1e-9) From 8403cb38cf14b8c887305cd3ff061f491e994d10 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 3 Mar 2026 13:20:53 +0100 Subject: [PATCH 16/32] fix mypy errors --- src/mrseq/scripts/epi2d_fid.py | 3 +++ src/mrseq/scripts/epi2d_se.py | 3 +++ src/mrseq/utils/EpiReadout.py | 3 +++ tests/scripts/test_epi2d_fid.py | 2 +- tests/scripts/test_epi2d_se.py | 2 +- 5 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index 12047d0..8461714 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -213,6 +213,7 @@ def epi2d_fid_kernel( # Write noise trajectory to MRD (zero trajectory — no gradients active) if mrd_header_file: + assert prot is not None n_samples = epi2d.adc.num_samples acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) @@ -262,6 +263,7 @@ def epi2d_fid_kernel( # Write navigator trajectory to MRD if mrd_header_file: + assert prot is not None n_samples = epi2d.adc.num_samples traj = np.zeros((n_samples, 2), dtype=np.float32) if nav_gx_sign > 0: @@ -306,6 +308,7 @@ def epi2d_fid_kernel( # close ISMRMRD file if mrd_header_file: + assert prot is not None prot.close() # set gridding definitions extracted from EpiReadout diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 27619e0..a032f6f 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -276,6 +276,7 @@ def epi2d_se_kernel( # Write noise trajectory to MRD (zero trajectory — no gradients active) if mrd_header_file: + assert prot is not None n_samples = epi2d.adc.num_samples acq = ismrmrd.Acquisition() acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) @@ -325,6 +326,7 @@ def epi2d_se_kernel( # Write navigator trajectory to MRD if mrd_header_file: + assert prot is not None n_samples = epi2d.adc.num_samples traj = np.zeros((n_samples, 2), dtype=np.float32) if nav_gx_sign > 0: @@ -370,6 +372,7 @@ def epi2d_se_kernel( # close ISMRMRD file if mrd_header_file: + assert prot is not None prot.close() # set gridding definitions extracted from EpiReadout diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index df13af3..fcfd632 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -414,6 +414,9 @@ def calculate_trajectory(self) -> tuple[np.ndarray, np.ndarray]: # For symmetric EPI, compute per-sample ky offset from the blip waveforms. # tt is relative to the waveform start; add delay to get block-relative times. if self.readout_type == 'symmetric': + assert self.gy_blipup is not None + assert self.gy_blipdown is not None + assert self.gy_blipdownup is not None bu_tt = self.gy_blipup.tt + self.gy_blipup.delay bd_tt = self.gy_blipdown.tt + self.gy_blipdown.delay bdu_tt = self.gy_blipdownup.tt + self.gy_blipdownup.delay diff --git a/tests/scripts/test_epi2d_fid.py b/tests/scripts/test_epi2d_fid.py index 0aba17f..758ccf5 100644 --- a/tests/scripts/test_epi2d_fid.py +++ b/tests/scripts/test_epi2d_fid.py @@ -48,7 +48,7 @@ def test_mrd_trajectory( fov = 200e-3 n_readout = 64 n_phase_encoding = 64 - oversampling = 2 + oversampling: Literal[1, 2, 4] = 2 n_slices = 1 mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_fid_traj.h5' diff --git a/tests/scripts/test_epi2d_se.py b/tests/scripts/test_epi2d_se.py index a0e3921..9eb694d 100644 --- a/tests/scripts/test_epi2d_se.py +++ b/tests/scripts/test_epi2d_se.py @@ -48,7 +48,7 @@ def test_mrd_trajectory( fov = 200e-3 n_readout = 64 n_phase_encoding = 64 - oversampling = 2 + oversampling: Literal[1, 2, 4] = 2 n_slices = 1 mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_se_traj.h5' From a53c9b0c9a70e2c995129306bd95e83da60a1b9b Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Tue, 3 Mar 2026 13:40:18 +0000 Subject: [PATCH 17/32] add direct recon for all --- examples/epi2d.ipynb | 38 +++++++++++--------------------------- 1 file changed, 11 insertions(+), 27 deletions(-) diff --git a/examples/epi2d.ipynb b/examples/epi2d.ipynb index 05646cf..d53f39c 100644 --- a/examples/epi2d.ipynb +++ b/examples/epi2d.ipynb @@ -34,13 +34,12 @@ "from mrpro.algorithms.reconstruction import DirectReconstruction\n", "from mrpro.data import KData\n", "from mrpro.data.traj_calculators import KTrajectoryIsmrmrd\n", - "from mrpro.operators import NonUniformFastFourierOp\n", "\n", "from mrseq.scripts.epi2d_fid import main as create_seq_fid\n", "from mrseq.scripts.epi2d_se import main as create_seq_se\n", "from mrseq.utils import combine_ismrmrd_files\n", "from mrseq.utils import sys_defaults\n", - "from mrseq.utils.trajectory import EpiReadout" + "from mrseq.utils.EpiReadout import EpiReadout" ] }, { @@ -231,7 +230,7 @@ "\n", "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", "recon = DirectReconstruction(kdata, csm=None)\n", - "idata_fid_sym_no_ramp = recon(kdata).data" + "idata_fid_sym_no_ramp = recon(kdata)" ] }, { @@ -269,13 +268,8 @@ "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", "\n", "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", - "fourier_op = NonUniformFastFourierOp(\n", - " direction=(-2, -1),\n", - " recon_matrix=kdata.header.recon_matrix[-2:],\n", - " encoding_matrix=kdata.header.encoding_matrix[-2:],\n", - " traj=kdata.traj,\n", - ")\n", - "idata_fid_sym_ramp = fourier_op.H(kdata.data)[0]" + "recon = DirectReconstruction(kdata, csm=None)\n", + "idata_fid_sym_ramp = recon(kdata)" ] }, { @@ -313,13 +307,8 @@ "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", "\n", "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", - "fourier_op = NonUniformFastFourierOp(\n", - " direction=(-2, -1),\n", - " recon_matrix=kdata.header.recon_matrix[-2:],\n", - " encoding_matrix=kdata.header.encoding_matrix[-2:],\n", - " traj=kdata.traj,\n", - ")\n", - "idata_fid_flybackramp = fourier_op.H(kdata.data)[0]" + "recon = DirectReconstruction(kdata, csm=None)\n", + "idata_fid_flyback_ramp = recon(kdata)" ] }, { @@ -356,13 +345,8 @@ "combine_ismrmrd_files(fname_mrd, Path(f'{fname_seq}_header.h5'))\n", "\n", "kdata = KData.from_file(str(fname_mrd).replace('.mrd', '_with_traj.mrd'), trajectory=KTrajectoryIsmrmrd())\n", - "fourier_op = NonUniformFastFourierOp(\n", - " direction=(-2, -1),\n", - " recon_matrix=kdata.header.recon_matrix[-2:],\n", - " encoding_matrix=kdata.header.encoding_matrix[-2:],\n", - " traj=kdata.traj,\n", - ")\n", - "idata_se_sym_ramp = fourier_op.H(kdata.data)[0]" + "recon = DirectReconstruction(kdata, csm=None)\n", + "idata_se_sym_ramp = recon(kdata)" ] }, { @@ -377,14 +361,14 @@ " cax.set_xticks([])\n", " cax.set_yticks([])\n", "\n", - "for idx, (idat, label) in enumerate(\n", + "for idx, (idata, label) in enumerate(\n", " zip(\n", - " [idata_fid_sym_no_ramp, idata_fid_sym_ramp, idata_fid_flybackramp],\n", + " [idata_fid_sym_no_ramp, idata_fid_sym_ramp, idata_fid_flyback_ramp],\n", " ['symmetric no ramp', 'symmetric ramp', 'flyback ramp'],\n", " strict=True,\n", " )\n", "):\n", - " idat = idat.abs().squeeze().numpy()\n", + " idat = idata.data.abs().squeeze().numpy()\n", " idat /= np.sort(idat.flatten())[int(idat.size * 0.99)]\n", " if idx == 0:\n", " idat_ref = idat\n", From 01f9dc59c309dee2fc6328b370adb5733cffe9c3 Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Tue, 3 Mar 2026 21:56:04 +0000 Subject: [PATCH 18/32] fix se epi2d --- examples/epi2d.ipynb | 15 ++++----------- src/mrseq/scripts/epi2d_se.py | 15 +++++---------- 2 files changed, 9 insertions(+), 21 deletions(-) diff --git a/examples/epi2d.ipynb b/examples/epi2d.ipynb index d53f39c..31b76ac 100644 --- a/examples/epi2d.ipynb +++ b/examples/epi2d.ipynb @@ -331,6 +331,7 @@ " test_report=False,\n", " timing_check=False,\n", " show_plots=False,\n", + " n_slices=1,\n", " fov_xy=float(phantom.size.numpy()[0]),\n", " n_readout=image_matrix_size[0],\n", " n_phase_encoding=image_matrix_size[0],\n", @@ -363,8 +364,8 @@ "\n", "for idx, (idata, label) in enumerate(\n", " zip(\n", - " [idata_fid_sym_no_ramp, idata_fid_sym_ramp, idata_fid_flyback_ramp],\n", - " ['symmetric no ramp', 'symmetric ramp', 'flyback ramp'],\n", + " [idata_fid_sym_no_ramp, idata_fid_sym_ramp, idata_fid_flyback_ramp, idata_se_sym_ramp],\n", + " ['symmetric no ramp', 'symmetric ramp', 'flyback ramp', 'se symmetric ramp'],\n", " strict=True,\n", " )\n", "):\n", @@ -382,16 +383,8 @@ "relative_error /= 3\n", "\n", "print(f'Relative error {relative_error}')\n", - "assert relative_error < 0.15" + "assert relative_error < 0.19" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "26", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 1802436..e59e8a8 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -280,30 +280,25 @@ def epi2d_se_kernel( # add navigator scans for ghost correction if add_navigator_acq: - # reverse the readout gradient and pre-winder in advance for navigator - gx_pre = pp.scale_grad(epi2d.gx_pre, -1) - gx = pp.scale_grad(epi2d.gx, -1) # add slice selection rewinder and readout pre-winder in x direction (gy_pre will be added after navigators) - gzr, gx_pre = pp.align(left=[gzr], right=[gx_pre]) + gzr, gx_pre = pp.align(left=[gzr], right=[epi2d.gx_pre]) seq.add_block( gzr, gx_pre, pp.make_label(label='NAV', type='SET', value=1), pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), ) - # reverse gx_pre back after adding to sequence - gx_pre = pp.scale_grad(gx_pre, -1) # add 3 navigator acquisitions for n in range(n_navigator_acq): + gx_sign = (-1) ** n seq.add_block( - gx, + pp.scale_grad(epi2d.gx, gx_sign), epi2d.adc, - pp.make_label(label='REV', type='SET', value=gx.amplitude < 0), - pp.make_label(label='SEG', type='SET', value=gx.amplitude < 0), + pp.make_label(label='REV', type='SET', value=gx_sign < 0), + pp.make_label(label='SEG', type='SET', value=gx_sign < 0), pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), ) - gx = pp.scale_grad(gx, -1) # add gy_pre and reset labels seq.add_block( From 468cac4857d99364ba830a9501e66b91201f8478 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 4 Mar 2026 08:37:00 +0100 Subject: [PATCH 19/32] adjust for v141 compatibility --- src/mrseq/scripts/epi2d_fid.py | 6 +++++- src/mrseq/scripts/epi2d_se.py | 6 +++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index 8461714..f67b7d4 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -10,6 +10,7 @@ from mrseq.utils import round_to_raster from mrseq.utils import sys_defaults +from mrseq.utils import write_sequence from mrseq.utils.EpiReadout import EpiReadout from mrseq.utils.ismrmrd import Fov from mrseq.utils.ismrmrd import Limits @@ -354,6 +355,7 @@ def main( show_plots: bool = True, test_report: bool = True, timing_check: bool = True, + v141_compatibility: bool = True, ) -> tuple[pp.Sequence, Path]: """Generate a 2D Echo Planar Imaging (EPI) FID sequence. @@ -395,6 +397,8 @@ def main( Toggles advanced test report. timing_check Toggles timing check of the sequence. + v141_compatibility + Save the sequence in pulseq v1.4.1 for backwards compatibility. Returns ------- @@ -475,7 +479,7 @@ def main( # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") - seq.write(str(output_path / filename), create_signature=True) + write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) if show_plots: seq.plot(time_range=(0, tr or _min_tr), plot_now=True) diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index dc6497a..b3ec22d 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -10,6 +10,7 @@ from mrseq.utils import round_to_raster from mrseq.utils import sys_defaults +from mrseq.utils import write_sequence from mrseq.utils.EpiReadout import EpiReadout from mrseq.utils.ismrmrd import Fov from mrseq.utils.ismrmrd import Limits @@ -412,6 +413,7 @@ def main( show_plots: bool = True, test_report: bool = True, timing_check: bool = True, + v141_compatibility: bool = True, ) -> tuple[pp.Sequence, Path]: """Generate a 2D Echo Planar Imaging (EPI) spin echo (SE) sequence. @@ -453,6 +455,8 @@ def main( Toggles advanced test report. timing_check Toggles timing check of the sequence. + v141_compatibility + Save the sequence in pulseq v1.4.1 for backwards compatibility. Returns ------- @@ -533,7 +537,7 @@ def main( # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") - seq.write(str(output_path / filename), create_signature=True) + write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) if show_plots: seq.plot(time_range=(0, 10 * (tr or min_tr)), plot_now=True) From 7e4cd030c0df3c16344abdb7abfcd10947426a01 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 4 Mar 2026 10:23:58 +0100 Subject: [PATCH 20/32] add calculate_full_trajectory method for better visualization in plot_trajectory --- src/mrseq/utils/EpiReadout.py | 328 ++++++++++++++++++++++++++-------- 1 file changed, 256 insertions(+), 72 deletions(-) diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index fcfd632..49d1f43 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -39,10 +39,12 @@ def _trapezoid_area_at_times( Accumulated gradient area at each sample time (in 1/m). """ t1 = rise_time - t2 = rise_time + flat_time + t2 = t1 + flat_time + t3 = t2 + fall_time area_at_t1 = amplitude * rise_time / 2 area_at_t2 = area_at_t1 + amplitude * flat_time + area_at_t3 = area_at_t2 + amplitude * fall_time / 2 area = np.zeros_like(sample_times) for i, t in enumerate(sample_times): @@ -52,25 +54,27 @@ def _trapezoid_area_at_times( area[i] = amplitude * t**2 / (2 * rise_time) elif t <= t2: area[i] = area_at_t1 + amplitude * (t - t1) - else: + elif t <= t3: dt = t - t2 area[i] = area_at_t2 + amplitude * dt - amplitude * dt**2 / (2 * fall_time) + else: + area[i] = area_at_t3 return area def _piecewise_linear_area_at_times( - wf_times: np.ndarray, - wf_values: np.ndarray, + waveform_times: np.ndarray, + waveform_values: np.ndarray, sample_times: np.ndarray, ) -> np.ndarray: """Compute cumulative integral of a piecewise-linear waveform at given sample times. Parameters ---------- - wf_times + waveform_times Time points of the piecewise-linear waveform (in seconds), sorted ascending. - wf_values + waveform_values Waveform amplitude values at each time point (in Hz/m). sample_times Times at which to evaluate the cumulative area (in seconds). @@ -80,26 +84,26 @@ def _piecewise_linear_area_at_times( area Cumulative area at each sample time (in 1/m). """ - n_knots = len(wf_times) + n_knots = len(waveform_times) cum_area = np.zeros(n_knots) for k in range(1, n_knots): - dt = wf_times[k] - wf_times[k - 1] - cum_area[k] = cum_area[k - 1] + (wf_values[k] + wf_values[k - 1]) / 2 * dt + dt = waveform_times[k] - waveform_times[k - 1] + cum_area[k] = cum_area[k - 1] + (waveform_values[k] + waveform_values[k - 1]) / 2 * dt area = np.zeros_like(sample_times) for i, t in enumerate(sample_times): - if t <= wf_times[0]: + if t <= waveform_times[0]: area[i] = 0.0 - elif t >= wf_times[-1]: + elif t >= waveform_times[-1]: area[i] = cum_area[-1] else: - k = np.searchsorted(wf_times, t, side='right') + k = np.searchsorted(waveform_times, t, side='right') base_area = cum_area[k - 1] - dt_seg = wf_times[k] - wf_times[k - 1] - dt = t - wf_times[k - 1] + dt_seg = waveform_times[k] - waveform_times[k - 1] + dt = t - waveform_times[k - 1] frac = dt / dt_seg - wf_at_t = wf_values[k - 1] + frac * (wf_values[k] - wf_values[k - 1]) - area[i] = base_area + (wf_values[k - 1] + wf_at_t) / 2 * dt + wf_at_t = waveform_values[k - 1] + frac * (waveform_values[k] - waveform_values[k - 1]) + area[i] = base_area + (waveform_values[k - 1] + wf_at_t) / 2 * dt return area @@ -381,6 +385,77 @@ def total_duration_without_prephaser(self) -> float: """Return total duration of readout excluding pre-phasers but including optional spoiler.""" return self.total_duration - pp.calc_duration(self.gx_pre, self.gy_pre) + def _readout_block_kxy( + self, + pe_idx: int, + sample_times: np.ndarray, + ky_base: float, + ) -> tuple[np.ndarray, np.ndarray, float]: + """Compute kx/ky at given time points for a single readout block. + + Parameters + ---------- + pe_idx + Phase encoding index (0-based). + sample_times + Time points relative to the start of the readout block (in seconds). + ky_base + Cumulative ky offset at the start of this block (in 1/m). + + Returns + ------- + kx + kx positions at each sample time (in 1/m). + ky + ky positions at each sample time (in 1/m). + ky_base_next + Updated ky_base for the next block (in 1/m). + """ + n_pe = self.n_phase_enc_total + kx_offset = self.gx_pre.area + + # kx: accumulated area under the readout gradient + kx_forward = _trapezoid_area_at_times( + self.gx.rise_time, self.gx.flat_time, self.gx.fall_time, abs(self.gx.amplitude), sample_times + ) + + if self.readout_type == 'flyback': + kx = kx_offset + kx_forward + ky = np.full_like(sample_times, ky_base) + ky_base_next = ky_base + self.gy_blip.area + + elif self.readout_type == 'symmetric': + if pe_idx % 2 == 0: + kx = kx_offset + kx_forward + else: + kx = kx_offset + self.gx.area - kx_forward + + # ky changes during readout due to blip overlap + if pe_idx == 0: + assert self.gy_blipup is not None + bu_tt = self.gy_blipup.tt + self.gy_blipup.delay + ky_blip = _piecewise_linear_area_at_times(bu_tt, self.gy_blipup.waveform, sample_times) + blip_total = _piecewise_linear_area_at_times(bu_tt, self.gy_blipup.waveform, np.array([bu_tt[-1]]))[0] + elif pe_idx == n_pe - 1: + assert self.gy_blipdown is not None + bd_tt = self.gy_blipdown.tt + self.gy_blipdown.delay + ky_blip = _piecewise_linear_area_at_times(bd_tt, self.gy_blipdown.waveform, sample_times) + blip_total = 0.0 + else: + assert self.gy_blipdownup is not None + bdu_tt = self.gy_blipdownup.tt + self.gy_blipdownup.delay + ky_blip = _piecewise_linear_area_at_times(bdu_tt, self.gy_blipdownup.waveform, sample_times) + blip_total = _piecewise_linear_area_at_times( + bdu_tt, self.gy_blipdownup.waveform, np.array([bdu_tt[-1]]) + )[0] + + ky = ky_base + ky_blip + ky_base_next = ky_base + blip_total + else: + raise NotImplementedError(f'Readout type "{self.readout_type}" is not supported.') + + return kx, ky, ky_base_next + def calculate_trajectory(self) -> tuple[np.ndarray, np.ndarray]: """Compute analytical k-space trajectory for the EPI readout. @@ -400,63 +475,121 @@ def calculate_trajectory(self) -> tuple[np.ndarray, np.ndarray]: # ADC sample times relative to the start of the gradient block (center of each dwell) sample_times = self.adc.delay + (np.arange(n_samples) + 0.5) * self.adc.dwell - # kx for a single readout line with positive gradient amplitude - kx_forward = _trapezoid_area_at_times( - self.gx.rise_time, self.gx.flat_time, self.gx.fall_time, abs(self.gx.amplitude), sample_times - ) - - # After gx_pre, kx starts at gx_pre.area - kx_offset = self.gx_pre.area + kx = np.zeros((n_pe, n_samples)) + ky = np.zeros((n_pe, n_samples)) + ky_base = self.gy_pre.area - # ky: pre-phaser moves to top of k-space - ky_start = self.gy_pre.area + for pe_idx in range(n_pe): + kx[pe_idx, :], ky[pe_idx, :], ky_base = self._readout_block_kxy(pe_idx, sample_times, ky_base) - # For symmetric EPI, compute per-sample ky offset from the blip waveforms. - # tt is relative to the waveform start; add delay to get block-relative times. - if self.readout_type == 'symmetric': - assert self.gy_blipup is not None - assert self.gy_blipdown is not None - assert self.gy_blipdownup is not None - bu_tt = self.gy_blipup.tt + self.gy_blipup.delay - bd_tt = self.gy_blipdown.tt + self.gy_blipdown.delay - bdu_tt = self.gy_blipdownup.tt + self.gy_blipdownup.delay + return kx, ky - ky_offset_blipup = _piecewise_linear_area_at_times(bu_tt, self.gy_blipup.waveform, sample_times) - ky_offset_blipdown = _piecewise_linear_area_at_times(bd_tt, self.gy_blipdown.waveform, sample_times) - ky_offset_blipdownup = _piecewise_linear_area_at_times(bdu_tt, self.gy_blipdownup.waveform, sample_times) + def calculate_full_trajectory(self, n_points_per_segment: int = 200) -> dict: + """Compute full analytical k-space trajectory including prephaser and transitions. - blipup_total = _piecewise_linear_area_at_times(bu_tt, self.gy_blipup.waveform, np.array([bu_tt[-1]]))[0] - blipdownup_total = _piecewise_linear_area_at_times( - bdu_tt, self.gy_blipdownup.waveform, np.array([bdu_tt[-1]]) - )[0] + This extends ``calculate_trajectory()`` by also computing the k-space path during + the prephaser gradients and the transitions between readout lines (flyback gradients + or symmetric ramp transitions). Useful for visualization. - kx = np.zeros((n_pe, n_samples)) - ky = np.zeros((n_pe, n_samples)) + Parameters + ---------- + n_points_per_segment + Number of densely sampled time points per gradient block for smooth trajectory lines. - # Track cumulative ky base (ky value at the START of each block) - ky_base = ky_start + Returns + ------- + result + Dictionary with keys: - for pe_idx in range(n_pe): - if self.readout_type == 'flyback': - kx[pe_idx, :] = kx_offset + kx_forward - ky[pe_idx, :] = ky_start + pe_idx * self.gy_blip.area + - ``'kx'``: full kx trajectory as 1-D array (1/m). + - ``'ky'``: full ky trajectory as 1-D array (1/m). + - ``'adc_indices'``: indices into kx/ky arrays that correspond to ADC samples. + """ + n_samples = self.adc.num_samples + n_pe = self.n_phase_enc_total - elif self.readout_type == 'symmetric': - if pe_idx % 2 == 0: - kx[pe_idx, :] = kx_offset + kx_forward - else: - kx[pe_idx, :] = kx_offset + self.gx.area - kx_forward + all_kx: list[np.ndarray] = [] + all_ky: list[np.ndarray] = [] + adc_indices: list[int] = [] + current_idx = 0 + + # --- Prephaser block --- + pre_dur = pp.calc_duration(self.gx_pre, self.gy_pre) + t_pre = np.linspace(0, pre_dur, n_points_per_segment, endpoint=False) + kx_pre = _trapezoid_area_at_times( + self.gx_pre.rise_time, + self.gx_pre.flat_time, + self.gx_pre.fall_time, + self.gx_pre.amplitude, + t_pre - self.gx_pre.delay, + ) + ky_pre = _trapezoid_area_at_times( + self.gy_pre.rise_time, + self.gy_pre.flat_time, + self.gy_pre.fall_time, + self.gy_pre.amplitude, + t_pre - self.gy_pre.delay, + ) + all_kx.append(kx_pre) + all_ky.append(ky_pre) + current_idx += len(t_pre) - if pe_idx == 0: - ky[pe_idx, :] = ky_base + ky_offset_blipup - ky_base += blipup_total - elif pe_idx == n_pe - 1: - ky[pe_idx, :] = ky_base + ky_offset_blipdown - else: - ky[pe_idx, :] = ky_base + ky_offset_blipdownup - ky_base += blipdownup_total + # Current kx/ky offsets at the end of the prephaser + kx_offset = self.gx_pre.area + ky_base = self.gy_pre.area - return kx, ky + for pe_idx in range(n_pe): + # --- Readout block (densely sampled) --- + block_dur = pp.calc_duration(self.gx) + t_block = np.linspace(0, block_dur, n_points_per_segment, endpoint=False) + + kx_block, ky_block, ky_base = self._readout_block_kxy(pe_idx, t_block, ky_base) + + all_kx.append(kx_block) + all_ky.append(ky_block) + + # Record ADC sample indices + for s in range(n_samples): + adc_indices.append(current_idx + np.searchsorted(t_block, self.adc.delay + (s + 0.5) * self.adc.dwell)) + current_idx += len(t_block) + + # Flyback + blip transition block (flyback only) + if self.readout_type == 'flyback' and pe_idx < n_pe - 1: + assert self.gx_flyback is not None + fb_dur = pp.calc_duration(self.gx_flyback, self.gy_blip) + t_fb = np.linspace(0, fb_dur, n_points_per_segment // 2, endpoint=False) + # kx during flyback + kx_fb = ( + kx_offset + + self.gx.area + + _trapezoid_area_at_times( + self.gx_flyback.rise_time, + self.gx_flyback.flat_time, + self.gx_flyback.fall_time, + self.gx_flyback.amplitude, + t_fb - self.gx_flyback.delay, + ) + ) + # ky during flyback + ky_pre_blip = ky_base - self.gy_blip.area + ky_fb = ky_pre_blip + _trapezoid_area_at_times( + self.gy_blip.rise_time, + self.gy_blip.flat_time, + self.gy_blip.fall_time, + self.gy_blip.amplitude, + t_fb - self.gy_blip.delay, + ) + all_kx.append(kx_fb) + all_ky.append(ky_fb) + current_idx += len(t_fb) + + kx_full = np.concatenate(all_kx) + ky_full = np.concatenate(all_ky) + adc_idx = np.array(adc_indices, dtype=int) + # Clip to valid range + adc_idx = np.clip(adc_idx, 0, len(kx_full) - 1) + + return {'kx': kx_full, 'ky': ky_full, 'adc_indices': adc_idx} def add_to_seq( self, @@ -530,16 +663,67 @@ def plot_sequence(self): ax.axvline(x=self.time_to_center, color='r', linestyle='--') plt.show() - def plot_trajectory(self): - """Plot k-space trajectory using the analytical expression.""" - kx, ky = self.calculate_trajectory() + def plot_trajectory(self, show_adc: bool = True): + """Plot full k-space trajectory with color-coded direction. + + The trajectory is colored from start (dark) to end (bright) using the + copper colormap. Small arrows at the midpoint of each readout line + indicate the traversal direction. ADC sample positions are shown as + red markers. + + Parameters + ---------- + show_adc + If True, ADC sample positions are shown as red markers. + + Returns + ------- + fig + Matplotlib figure object. + """ + from matplotlib.collections import LineCollection + + traj = self.calculate_full_trajectory() + kx = traj['kx'] + ky = traj['ky'] + adc_idx = traj['adc_indices'] + + fig, ax = plt.subplots() + + # Create line segments colored by traversal order + points = np.column_stack([kx, ky]).reshape(-1, 1, 2) + segments = np.concatenate([points[:-1], points[1:]], axis=1) + t = np.linspace(0, 1, len(segments)) + lc = LineCollection(segments, cmap='copper', linewidth=1.0) + lc.set_array(t) + ax.add_collection(lc) + + # ADC samples + if show_adc: + ax.plot(kx[adc_idx], ky[adc_idx], 'x', color='red', markersize=3, zorder=3) + + # Direction arrows at the midpoint of each readout line + n_samples = self.adc.num_samples + for pe_idx in range(self.n_phase_enc_total): + start = pe_idx * n_samples + mid = start + n_samples // 2 + kx_adc = kx[adc_idx] + ky_adc = ky[adc_idx] + dx = kx_adc[mid] - kx_adc[mid - 1] + dy = ky_adc[mid] - ky_adc[mid - 1] + ax.annotate( + '', + xy=(kx_adc[mid] + dx, ky_adc[mid] + dy), + xytext=(kx_adc[mid], ky_adc[mid]), + arrowprops={'arrowstyle': '->', 'color': 'blue', 'lw': 1.5}, + ) - fig = plt.figure() - plt.plot(kx.T, ky.T, 'b') - plt.plot(kx.ravel(), ky.ravel(), 'x', color='red', markersize=4) - plt.xlabel('kx (1/m)') - plt.ylabel('ky (1/m)') - plt.grid() + ax.set_xlabel('kx (1/m)') + ax.set_ylabel('ky (1/m)') + ax.set_aspect('equal', adjustable='datalim') + ax.autoscale_view() + ax.grid(True, alpha=0.3) + plt.tight_layout() plt.show() return fig From 74f42beef13e0c0da74a29adf26806881ff13723 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 4 Mar 2026 10:27:30 +0100 Subject: [PATCH 21/32] fix mypy error --- src/mrseq/utils/EpiReadout.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index 49d1f43..d740d00 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -694,7 +694,7 @@ def plot_trajectory(self, show_adc: bool = True): points = np.column_stack([kx, ky]).reshape(-1, 1, 2) segments = np.concatenate([points[:-1], points[1:]], axis=1) t = np.linspace(0, 1, len(segments)) - lc = LineCollection(segments, cmap='copper', linewidth=1.0) + lc = LineCollection(segments.tolist(), cmap='copper', linewidth=1.0) lc.set_array(t) ax.add_collection(lc) From dc500da9efb8df00c3c586f990821ad613c0654e Mon Sep 17 00:00:00 2001 From: Christoph Kolbitsch Date: Fri, 20 Mar 2026 19:00:45 +0000 Subject: [PATCH 22/32] ensures epi object stays unchanged --- src/mrseq/utils/EpiReadout.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index d740d00..d7f112c 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -614,8 +614,6 @@ def add_to_seq( seq.add_block(self.gx_pre, self.gy_pre) for pe_idx in range(self.n_phase_enc_total): - rev_label = pp.make_label(type='SET', label='REV', value=self.gx.amplitude < 0) - seg_label = pp.make_label(type='SET', label='SEG', value=self.gx.amplitude < 0) if self.readout_type == 'symmetric': # Select blip gradient based on phase encoding index if pe_idx == 0: @@ -626,11 +624,13 @@ def add_to_seq( gy_blip = self.gy_blipdownup # Add readout block and reverse polarity of readout gradient - seq.add_block(self.gx, gy_blip, self.adc, lin_label, rev_label, seg_label) - self.gx.amplitude = -self.gx.amplitude + gx_scale = -1 if np.mod(pe_idx, 2) else 1 + rev_label = pp.make_label(type='SET', label='REV', value=gx_scale == -1) + seg_label = pp.make_label(type='SET', label='SEG', value=gx_scale == -1) + seq.add_block(pp.scale_grad(self.gx, gx_scale), gy_blip, self.adc, lin_label, rev_label, seg_label) elif self.readout_type == 'flyback': - seq.add_block(self.gx, self.adc, lin_label, rev_label) + seq.add_block(self.gx, self.adc, lin_label) if pe_idx != self.n_phase_enc_total - 1: seq.add_block(self.gx_flyback, self.gy_blip) From 97dd12dc7a141c82b1eb8fdb6a2ddea065ae4e13 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 7 Apr 2026 14:15:48 +0200 Subject: [PATCH 23/32] adjust epi2d_se default values to match epi2d_fid --- src/mrseq/scripts/epi2d_se.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index b3ec22d..8ef43a0 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -401,13 +401,13 @@ def main( fov_xy: float = 200e-3, n_readout: int = 64, n_phase_encoding: int = 64, - n_slices: int = 3, + n_slices: int = 1, slice_thickness: float = 8e-3, bandwidth: float = 100e3, readout_type: Literal['symmetric', 'flyback'] = 'symmetric', readout_oversampling: Literal[1, 2, 4] = 2, ramp_sampling: bool = True, - partial_fourier_factor: float = 1.0, + partial_fourier_factor: float = 0.7, add_navigator_acq: bool = True, add_noise_acq: bool = True, show_plots: bool = True, From cbaf711711b43b62773affdc8e64681b59c29faa Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 19 May 2026 10:51:27 +0200 Subject: [PATCH 24/32] adjust expected dur to new default settings --- tests/scripts/test_epi2d_se.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/scripts/test_epi2d_se.py b/tests/scripts/test_epi2d_se.py index 9eb694d..78162d8 100644 --- a/tests/scripts/test_epi2d_se.py +++ b/tests/scripts/test_epi2d_se.py @@ -11,7 +11,7 @@ from mrseq.scripts.epi2d_se import main as create_seq from mrseq.utils.EpiReadout import EpiReadout -EXPECTED_DUR = 0.22586 # defined 2026-02-28 +EXPECTED_DUR = 0.04846 # defined 2026-05-19 def test_default_seq_duration(system_defaults): From 469449622e16c4fea431dff18562cf55576d3537 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 19 May 2026 10:57:44 +0200 Subject: [PATCH 25/32] change phase encoding to start with negative ky values --- src/mrseq/utils/EpiReadout.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index d7f112c..cae7c14 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -224,7 +224,7 @@ def __init__( # Create blip gradient with shortest possible timing gy_blip_duration = np.ceil(2 * np.sqrt(delta_ky / system.max_slew) / 10e-6 / 2) * 10e-6 * 2 gy_blip_half_dur = gy_blip_duration / 2 - self.gy_blip = pp.make_trapezoid(channel='y', system=self.system, area=-delta_ky, duration=gy_blip_duration) + self.gy_blip = pp.make_trapezoid(channel='y', system=self.system, area=delta_ky, duration=gy_blip_duration) # Create readout gradient gx_encoding_area = n_readout * delta_kx * oversampling @@ -302,7 +302,7 @@ def __init__( self.gy_pre = pp.make_trapezoid( channel='y', system=self.system, - area=self.n_phase_enc_pre_center * delta_ky, + area=-self.n_phase_enc_pre_center * delta_ky, ) self.gx_pre, self.gy_pre = pp.align(right=[self.gx_pre, self.gy_pre]) From f1cb387bd4ecc63f2f1bfd86cfa3859add5fb124 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 19 May 2026 11:01:12 +0200 Subject: [PATCH 26/32] refactor navigator code: - unify code for FID and SE sequence - move navigator code from kernels to EpiReadout class - allow arbitrary number of navigator scans --- src/mrseq/scripts/epi2d_fid.py | 67 +++++----------------------- src/mrseq/scripts/epi2d_se.py | 57 ++++-------------------- src/mrseq/utils/EpiReadout.py | 79 ++++++++++++++++++++++++++++++++++ 3 files changed, 97 insertions(+), 106 deletions(-) diff --git a/src/mrseq/scripts/epi2d_fid.py b/src/mrseq/scripts/epi2d_fid.py index f67b7d4..0addece 100644 --- a/src/mrseq/scripts/epi2d_fid.py +++ b/src/mrseq/scripts/epi2d_fid.py @@ -1,6 +1,5 @@ """2D Echo Planar Imaging (EPI) FID sequence.""" -from math import floor from pathlib import Path from typing import Literal @@ -140,7 +139,7 @@ def epi2d_fid_kernel( min_te += max(rf.ringdown_time, gz.fall_time) # RF ringdown time or gradient fall time if add_navigator_acq: min_te += pp.calc_duration(gzr, epi2d.gx_pre) - min_te += 3 * pp.calc_duration(epi2d.gx) + min_te += n_navigator_acq * pp.calc_duration(epi2d.gx) min_te += pp.calc_duration(epi2d.gy_pre) else: min_te += pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gx_pre) @@ -157,7 +156,7 @@ def epi2d_fid_kernel( min_tr = pp.calc_duration(rf, gz) if add_navigator_acq: min_tr += pp.calc_duration(gzr, epi2d.gx_pre) - min_tr += 3 * pp.calc_duration(epi2d.gx) + min_tr += n_navigator_acq * pp.calc_duration(epi2d.gx) min_tr += pp.calc_duration(epi2d.gy_pre) else: min_tr += pp.calc_duration(gzr, epi2d.gx_pre, epi2d.gx_pre) @@ -190,15 +189,6 @@ def epi2d_fid_kernel( prot = ismrmrd.Dataset(mrd_header_file, 'w') prot.write_xml_header(hdr.toXML('utf-8')) - # Precompute analytical navigator trajectory (single kx line, no ky blip) - if add_navigator_acq: - from mrseq.utils.EpiReadout import _trapezoid_area_at_times - - nav_sample_times = epi2d.adc.delay + (np.arange(epi2d.adc.num_samples) + 0.5) * epi2d.adc.dwell - nav_kx_forward = _trapezoid_area_at_times( - epi2d.gx.rise_time, epi2d.gx.flat_time, epi2d.gx.fall_time, abs(epi2d.gx.amplitude), nav_sample_times - ) - # obtain noise samples if selected if add_noise_acq: seq.add_block( @@ -234,53 +224,16 @@ def epi2d_fid_kernel( # add navigator scans for ghost correction if add_navigator_acq: - # reverse the readout gradient and pre-winder in advance for navigator - gx_pre = pp.scale_grad(epi2d.gx_pre, -1) - gx = pp.scale_grad(epi2d.gx, -1) - # add slice selection rewinder and readout pre-winder in x direction (gy_pre will be added after navigators) - gzr, gx_pre = pp.align(left=[gzr], right=[gx_pre]) - seq.add_block( + # FID: invert pre-winder for odd number of navigators + gx_pre_factor = -1 if n_navigator_acq % 2 == 1 else 1 + seq, prot = epi2d.add_navigator_to_seq( + seq, gzr, - gx_pre, - pp.make_label(label='NAV', type='SET', value=1), - pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), + n_navigator_acq, + n_phase_encoding, + gx_pre_factor, + mrd_dataset=prot, ) - # reverse gx_pre back after adding to sequence - gx_pre = pp.scale_grad(gx_pre, -1) - - # Navigator kx offset: starts from -gx_pre.area (reversed pre-winder) - nav_kx_offset = -epi2d.gx_pre.area - nav_gx_sign = -1.0 # first navigator uses reversed gx - - # add navigator scans for ghost correction - for n in range(n_navigator_acq): - seq.add_block( - gx, - epi2d.adc, - pp.make_label(label='REV', type='SET', value=gx.amplitude < 0), - pp.make_label(label='SEG', type='SET', value=gx.amplitude < 0), - pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), - ) - - # Write navigator trajectory to MRD - if mrd_header_file: - assert prot is not None - n_samples = epi2d.adc.num_samples - traj = np.zeros((n_samples, 2), dtype=np.float32) - if nav_gx_sign > 0: - nav_kx = nav_kx_offset + nav_kx_forward - else: - nav_kx = nav_kx_offset - nav_kx_forward - traj[:, 0] = nav_kx * fov_xy * readout_oversampling - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) - acq.traj[:] = traj - prot.append_acquisition(acq) - - # Update kx offset for next navigator (accumulate the signed gx area) - nav_kx_offset += nav_gx_sign * epi2d.gx.area - gx = pp.scale_grad(gx, -1) - nav_gx_sign = -nav_gx_sign # add echo time delay seq.add_block(pp.make_delay(te_delay)) diff --git a/src/mrseq/scripts/epi2d_se.py b/src/mrseq/scripts/epi2d_se.py index 8ef43a0..d02aaaf 100644 --- a/src/mrseq/scripts/epi2d_se.py +++ b/src/mrseq/scripts/epi2d_se.py @@ -1,6 +1,5 @@ """2D Echo Planar Imaging (EPI) spin echo (SE) sequence.""" -from math import floor from pathlib import Path from typing import Literal @@ -253,15 +252,6 @@ def epi2d_se_kernel( prot = ismrmrd.Dataset(mrd_header_file, 'w') prot.write_xml_header(hdr.toXML('utf-8')) - # Precompute analytical navigator trajectory (single kx line, no ky blip) - if add_navigator_acq: - from mrseq.utils.EpiReadout import _trapezoid_area_at_times - - nav_sample_times = epi2d.adc.delay + (np.arange(epi2d.adc.num_samples) + 0.5) * epi2d.adc.dwell - nav_kx_forward = _trapezoid_area_at_times( - epi2d.gx.rise_time, epi2d.gx.flat_time, epi2d.gx.fall_time, abs(epi2d.gx.amplitude), nav_sample_times - ) - # obtain noise samples if selected if add_noise_acq: seq.add_block( @@ -297,48 +287,17 @@ def epi2d_se_kernel( # add navigator scans for ghost correction if add_navigator_acq: - # add slice selection rewinder and readout pre-winder in x direction (gy_pre will be added after navigators) - gzr, gx_pre = pp.align(left=[gzr], right=[epi2d.gx_pre]) - seq.add_block( + # SE: invert pre-winder for even number of navigators + gx_pre_factor = -1 if n_navigator_acq % 2 == 0 else 1 + seq, prot = epi2d.add_navigator_to_seq( + seq, gzr, - gx_pre, - pp.make_label(label='NAV', type='SET', value=1), - pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), + n_navigator_acq, + n_phase_encoding, + gx_pre_factor, + mrd_dataset=prot, ) - # Navigator kx offset: starts from -gx_pre.area (reversed pre-winder) - nav_kx_offset = -epi2d.gx_pre.area - nav_gx_sign = -1.0 # first navigator uses reversed gx - - # add 3 navigator acquisitions - for n in range(n_navigator_acq): - gx_sign = (-1) ** n - seq.add_block( - pp.scale_grad(epi2d.gx, gx_sign), - epi2d.adc, - pp.make_label(label='REV', type='SET', value=gx_sign < 0), - pp.make_label(label='SEG', type='SET', value=gx_sign < 0), - pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), - ) - - # Write navigator trajectory to MRD - if mrd_header_file: - assert prot is not None - n_samples = epi2d.adc.num_samples - traj = np.zeros((n_samples, 2), dtype=np.float32) - if nav_gx_sign > 0: - nav_kx = nav_kx_offset + nav_kx_forward - else: - nav_kx = nav_kx_offset - nav_kx_forward - traj[:, 0] = nav_kx * fov_xy * readout_oversampling - acq = ismrmrd.Acquisition() - acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) - acq.traj[:] = traj - prot.append_acquisition(acq) - - # Update kx offset for next navigator (accumulate the signed gx area) - nav_kx_offset += gx_sign * epi2d.gx.area - # add gy_pre and reset labels seq.add_block( pp.scale_grad(epi2d.gy_pre, -1), diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index cae7c14..8be18c8 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -591,6 +591,85 @@ def calculate_full_trajectory(self, n_points_per_segment: int = 200) -> dict: return {'kx': kx_full, 'ky': ky_full, 'adc_indices': adc_idx} + def add_navigator_to_seq( + self, + seq: pp.Sequence, + gz_rewinder, + n_navigator_acq: int, + n_phase_encoding: int, + gx_pre_factor: int, + mrd_dataset: ismrmrd.Dataset | None = None, + ) -> tuple[pp.Sequence, ismrmrd.Dataset | None]: + """Add navigator acquisitions to the sequence for EPI ghost correction. + + Parameters + ---------- + seq + PyPulseq Sequence object. + gz_rewinder + Slice-selection rewinder gradient. + n_navigator_acq + Number of navigator acquisitions. + n_phase_encoding + Number of phase encoding steps (used for LIN label). + gx_pre_factor + Sign factor for the readout pre-winder. Use ``-1`` for odd navigator + count (FID) or even navigator count (SE) to ensure the kx position + is correct at the start of the imaging readout. + mrd_dataset + ISMRMRD dataset for writing navigator trajectories. If None, no + trajectory is written. + """ + from math import floor + + # Align slice-selection rewinder and readout pre-winder + gz_rewinder, gx_pre = pp.align(left=[gz_rewinder], right=[self.gx_pre]) + seq.add_block( + gz_rewinder, + pp.scale_grad(gx_pre, gx_pre_factor), + pp.make_label(label='NAV', type='SET', value=1), + pp.make_label(label='LIN', type='SET', value=floor(n_phase_encoding / 2)), + ) + + # Precompute analytical navigator kx trajectory (single line, no ky blip) + nav_sample_times = self.adc.delay + (np.arange(self.adc.num_samples) + 0.5) * self.adc.dwell + nav_kx_forward = _trapezoid_area_at_times( + self.gx.rise_time, self.gx.flat_time, self.gx.fall_time, abs(self.gx.amplitude), nav_sample_times + ) + + # Navigator kx offset: starts from gx_pre.area scaled by gx_pre_factor + nav_kx_offset = self.gx_pre.area * gx_pre_factor + + # Add navigator acquisitions + for n in range(n_navigator_acq): + gx_sign = gx_pre_factor * (-1) ** n + seq.add_block( + pp.scale_grad(self.gx, gx_sign), + self.adc, + pp.make_label(label='REV', type='SET', value=gx_sign < 0), + pp.make_label(label='SEG', type='SET', value=gx_sign < 0), + pp.make_label(label='AVG', type='SET', value=(n + 1) == 3), + ) + + # Write navigator trajectory to MRD + if mrd_dataset is not None: + n_samples = self.adc.num_samples + traj = np.zeros((n_samples, 2), dtype=np.float32) + if gx_sign > 0: + nav_kx = nav_kx_offset + nav_kx_forward + else: + nav_kx = nav_kx_offset - nav_kx_forward + traj[:, 0] = nav_kx * self.fov * self.oversampling + acq = ismrmrd.Acquisition() + acq.resize(trajectory_dimensions=2, number_of_samples=n_samples) + acq.traj[:] = traj + mrd_dataset.append_acquisition(acq) + + # Update kx offset for next navigator + nav_kx_offset += gx_sign * self.gx.area + + return seq, mrd_dataset + def add_to_seq( self, seq: pp.Sequence, From 79401ebc13a51effccfcff605f33c909f5fded87 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 19 May 2026 11:12:48 +0200 Subject: [PATCH 27/32] temporarily limit ismrmrd version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index dc236f9..84ef7e1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,7 +42,7 @@ dependencies = [ "numpy>=1.23,<3.0", "matplotlib<4.0", "pypulseq>=1.4.1", - "ismrmrd>=1.14.1", + "ismrmrd>=1.14.1,<1.15", ] [project.optional-dependencies] From 10627fce688c2f67690325746d3578f1745e9ef5 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 19 May 2026 11:37:46 +0200 Subject: [PATCH 28/32] update example notebook --- examples/epi2d.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/epi2d.ipynb b/examples/epi2d.ipynb index 31b76ac..5fc4a78 100644 --- a/examples/epi2d.ipynb +++ b/examples/epi2d.ipynb @@ -383,7 +383,7 @@ "relative_error /= 3\n", "\n", "print(f'Relative error {relative_error}')\n", - "assert relative_error < 0.19" + "assert relative_error < 0.20" ] } ], From 678740d69459d8331d2425e22f64b3dbe2305374 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 20 May 2026 11:18:35 +0200 Subject: [PATCH 29/32] add contact and acknowledgment section to README --- README.md | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 1ea29ca..594b5cb 100644 --- a/README.md +++ b/README.md @@ -15,9 +15,9 @@ We are looking forward to your contributions via Pull-Requests. -### Installation for developers +## Installation for developers -#### Prerequisites for Windows +### Prerequisites for Windows Before installing `MRseq` with development dependencies on Windows, you need: @@ -29,7 +29,7 @@ Before installing `MRseq` with development dependencies on Windows, you need: 2. **Rust toolchain** (automatically installed by MRzeroCore if not present) -#### Installation steps +### Installation steps 1. Clone the `MRseq` repository 2. Create/select a python environment @@ -37,7 +37,7 @@ Before installing `MRseq` with development dependencies on Windows, you need: 4. Setup pre-commit hook: ``` pre-commit install ``` -### Licensing and Dependencies +## Licensing and Dependencies The core source code of `MRseq` is licensed under the **Apache-2.0 license**. @@ -45,3 +45,21 @@ For quality control and demonstration purposes, the notebooks in the `examples/` **Important for Commercial Users:** The AGPL dependency is listed only as an optional dev dependency. Installing this package via `pip install mrseq` **does not** install the AGPL library. Your usage of the core library remains subject to the permissive Apache-2.0 terms, unaffected by the licensing of the test suite or examples. + +## Contact +For questions or support, please contact us: + +> Christoph Kolbitsch +> Physikalisch-Technische Bundesanstalt (PTB) +> Abbestraße 2-12 +> 10587 Berlin, Germany +> christoph.kolbitsch@ptb.de + +> Patrick Schünke +> NVision Quantum Technologies GmbH +> Wolfgang-Paul-Straße 2 +> 89081 Ulm, Germany +> patrick.schuenke@nvision-quantum.com + +## Acknowledgments +We would like to thank [NVision Quantum Technologies GmbH](https://www.nvision-quantum.com/) for their contributions to this project. From 2be0697223d3a63197904acc566803b4fc5346b2 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Fri, 22 May 2026 08:46:26 +0200 Subject: [PATCH 30/32] extend trajectory tests --- tests/scripts/test_epi2d_fid.py | 39 +++++++++++++------------ tests/scripts/test_epi2d_se.py | 51 ++++++++++++++------------------- 2 files changed, 42 insertions(+), 48 deletions(-) diff --git a/tests/scripts/test_epi2d_fid.py b/tests/scripts/test_epi2d_fid.py index 758ccf5..30de95f 100644 --- a/tests/scripts/test_epi2d_fid.py +++ b/tests/scripts/test_epi2d_fid.py @@ -33,23 +33,27 @@ def test_seq_creation_error_on_short_tr(system_defaults): create_seq(system=system_defaults, tr=2e-3, show_plots=False) +@pytest.mark.parametrize('add_noise_acq', [True, False]) @pytest.mark.parametrize('readout_type', ('flyback', 'symmetric')) -@pytest.mark.parametrize( - 'add_navigator_acq, add_noise_acq', - [(False, False), (True, False), (False, True), (True, True)], -) +@pytest.mark.parametrize('add_navigator_acq', [True, False]) +@pytest.mark.parametrize('partial_fourier_factor', [0.75, 1.0]) +@pytest.mark.parametrize('ramp_sampling', [True, False]) +@pytest.mark.parametrize('oversampling', [1, 2]) +@pytest.mark.parametrize('n_slices', [1, 4]) def test_mrd_trajectory( system_defaults, readout_type: Literal['flyback', 'symmetric'], add_navigator_acq: bool, add_noise_acq: bool, + partial_fourier_factor: float, + ramp_sampling: bool, + oversampling: Literal[1, 2], + n_slices: int, ): """Test that the MRD trajectory matches the analytical and PyPulseq trajectories.""" fov = 200e-3 - n_readout = 64 - n_phase_encoding = 64 - oversampling: Literal[1, 2, 4] = 2 - n_slices = 1 + n_readout = 32 + n_phase_encoding = 32 mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_fid_traj.h5' if mrd_file.exists(): @@ -71,8 +75,8 @@ def test_mrd_trajectory( rf_apodization=0.5, readout_type=readout_type, readout_oversampling=oversampling, - ramp_sampling=True, - partial_fourier_factor=1.0, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, add_spoiler=True, add_noise_acq=add_noise_acq, add_navigator_acq=add_navigator_acq, @@ -94,9 +98,9 @@ def test_mrd_trajectory( n_phase_encoding=n_phase_encoding, bandwidth=100e3, oversampling=oversampling, - ramp_sampling=True, + ramp_sampling=ramp_sampling, readout_type=readout_type, - partial_fourier_factor=1.0, + partial_fourier_factor=partial_fourier_factor, pe_enable=True, spoiling_enable=True, ) @@ -107,18 +111,17 @@ def test_mrd_trajectory( n_epi = epi2d.n_phase_enc_total assert n_acq == n_noise + (n_pre_per_slice + n_epi) * n_slices - # Compare noise + navigator against calculate_kspace + # Compare Pulseq trajectories against MRD trajectories including noise and navigator acquisitions k_traj_adc, _, _, _, _ = seq.calculate_kspace() - n_pre_total = n_noise + n_pre_per_slice * n_slices - for i in range(n_pre_total): + for i in range(n_acq): start = i * n_samples end = start + n_samples kx_calc = (k_traj_adc[0, start:end] * fov * oversampling).astype(np.float32) ky_calc = (k_traj_adc[1, start:end] * fov).astype(np.float32) - np.testing.assert_array_equal(traj_mrd[i][:, 0], kx_calc) - np.testing.assert_array_equal(traj_mrd[i][:, 1], ky_calc) + np.testing.assert_array_almost_equal(traj_mrd[i][:, 0], kx_calc, decimal=10) + np.testing.assert_array_almost_equal(traj_mrd[i][:, 1], ky_calc, decimal=10) - # Compare EPI readout lines against analytical trajectory + # Compare analytical trajectories against MRD trajectories without noise and navigator acquisitions kx_analytical, ky_analytical = epi2d.calculate_trajectory() for s in range(n_slices): for pe_idx in range(n_epi): diff --git a/tests/scripts/test_epi2d_se.py b/tests/scripts/test_epi2d_se.py index 78162d8..b716a5b 100644 --- a/tests/scripts/test_epi2d_se.py +++ b/tests/scripts/test_epi2d_se.py @@ -33,23 +33,27 @@ def test_seq_creation_error_on_short_tr(system_defaults): create_seq(system=system_defaults, tr=2e-3, show_plots=False) +@pytest.mark.parametrize('add_noise_acq', [True, False]) @pytest.mark.parametrize('readout_type', ('flyback', 'symmetric')) -@pytest.mark.parametrize( - 'add_navigator_acq, add_noise_acq', - [(False, False), (True, False), (False, True), (True, True)], -) +@pytest.mark.parametrize('add_navigator_acq', [True, False]) +@pytest.mark.parametrize('partial_fourier_factor', [0.75, 1.0]) +@pytest.mark.parametrize('ramp_sampling', [True, False]) +@pytest.mark.parametrize('oversampling', [1, 2]) +@pytest.mark.parametrize('n_slices', [1, 4]) def test_mrd_trajectory( system_defaults, readout_type: Literal['flyback', 'symmetric'], add_navigator_acq: bool, add_noise_acq: bool, + partial_fourier_factor: float, + ramp_sampling: bool, + oversampling: Literal[1, 2], + n_slices: int, ): """Test that the MRD trajectory matches the analytical and PyPulseq trajectories.""" fov = 200e-3 - n_readout = 64 - n_phase_encoding = 64 - oversampling: Literal[1, 2, 4] = 2 - n_slices = 1 + n_readout = 32 + n_phase_encoding = 32 mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_se_traj.h5' if mrd_file.exists(): @@ -71,8 +75,8 @@ def test_mrd_trajectory( rf_apodization=0.5, readout_type=readout_type, readout_oversampling=oversampling, - ramp_sampling=True, - partial_fourier_factor=1.0, + ramp_sampling=ramp_sampling, + partial_fourier_factor=partial_fourier_factor, add_spoiler=True, add_noise_acq=add_noise_acq, add_navigator_acq=add_navigator_acq, @@ -94,9 +98,9 @@ def test_mrd_trajectory( n_phase_encoding=n_phase_encoding, bandwidth=100e3, oversampling=oversampling, - ramp_sampling=True, + ramp_sampling=ramp_sampling, readout_type=readout_type, - partial_fourier_factor=1.0, + partial_fourier_factor=partial_fourier_factor, pe_enable=True, spoiling_enable=True, ) @@ -107,30 +111,17 @@ def test_mrd_trajectory( n_epi = epi2d.n_phase_enc_total assert n_acq == n_noise + (n_pre_per_slice + n_epi) * n_slices - # Compare noise + navigator against calculate_kspace + # Compare Pulseq trajectories against MRD trajectories including noise and navigator acquisitions k_traj_adc, _, _, _, _ = seq.calculate_kspace() - n_pre_total = n_noise + n_pre_per_slice * n_slices - for i in range(n_pre_total): + for i in range(n_acq): start = i * n_samples end = start + n_samples kx_calc = (k_traj_adc[0, start:end] * fov * oversampling).astype(np.float32) ky_calc = (k_traj_adc[1, start:end] * fov).astype(np.float32) - np.testing.assert_array_equal(traj_mrd[i][:, 0], kx_calc) - np.testing.assert_array_equal(traj_mrd[i][:, 1], ky_calc) + np.testing.assert_array_almost_equal(traj_mrd[i][:, 0], kx_calc, decimal=10) + np.testing.assert_array_almost_equal(traj_mrd[i][:, 1], ky_calc, decimal=10) - # Compare navigator trajectories per slice against calculate_kspace - for s in range(n_slices): - for n in range(n_nav): - mrd_idx = n_noise + s * (n_nav + n_epi) + n - calc_idx = n_noise + s * (n_nav + n_epi) + n - start = calc_idx * n_samples - end = start + n_samples - kx_calc = (k_traj_adc[0, start:end] * fov * oversampling).astype(np.float32) - ky_calc = (k_traj_adc[1, start:end] * fov).astype(np.float32) - np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 0], kx_calc) - np.testing.assert_array_equal(traj_mrd[mrd_idx][:, 1], ky_calc) - - # Compare EPI readout lines against analytical trajectory + # Compare analytical trajectories against MRD trajectories without noise and navigator acquisitions kx_analytical, ky_analytical = epi2d.calculate_trajectory() for s in range(n_slices): for pe_idx in range(n_epi): From 2c008449418e4856772398ecde5d6b25f34ab50b Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Fri, 22 May 2026 10:41:19 +0200 Subject: [PATCH 31/32] change to symmetric encoding in frequency direction --- src/mrseq/utils/EpiReadout.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mrseq/utils/EpiReadout.py b/src/mrseq/utils/EpiReadout.py index 8be18c8..252b430 100644 --- a/src/mrseq/utils/EpiReadout.py +++ b/src/mrseq/utils/EpiReadout.py @@ -297,7 +297,7 @@ def __init__( self.gx_pre = pp.make_trapezoid( channel='x', system=self.system, - area=-self.gx.area / 2 - delta_kx / 2, + area=-self.gx.area / 2, ) self.gy_pre = pp.make_trapezoid( channel='y', From 83d7f19362cf11a29bd22c4c49a67df92c7a3697 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Fri, 22 May 2026 11:09:59 +0200 Subject: [PATCH 32/32] use pytest tmp_path fixture to avoid file related errors --- tests/scripts/test_epi2d_fid.py | 8 ++------ tests/scripts/test_epi2d_se.py | 8 ++------ 2 files changed, 4 insertions(+), 12 deletions(-) diff --git a/tests/scripts/test_epi2d_fid.py b/tests/scripts/test_epi2d_fid.py index 30de95f..6e13ecb 100644 --- a/tests/scripts/test_epi2d_fid.py +++ b/tests/scripts/test_epi2d_fid.py @@ -1,7 +1,5 @@ """Tests for EPI FID sequence.""" -import tempfile -from pathlib import Path from typing import Literal import ismrmrd @@ -49,15 +47,14 @@ def test_mrd_trajectory( ramp_sampling: bool, oversampling: Literal[1, 2], n_slices: int, + tmp_path, ): """Test that the MRD trajectory matches the analytical and PyPulseq trajectories.""" fov = 200e-3 n_readout = 32 n_phase_encoding = 32 - mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_fid_traj.h5' - if mrd_file.exists(): - mrd_file.unlink() + mrd_file = tmp_path / 'test_epi2d_fid_traj.h5' seq, _, _ = epi2d_fid_kernel( system=system_defaults, @@ -88,7 +85,6 @@ def test_mrd_trajectory( n_acq = ds.number_of_acquisitions() traj_mrd = [ds.read_acquisition(i).traj.copy() for i in range(n_acq)] ds.close() - mrd_file.unlink() # Recreate EpiReadout for reference epi2d = EpiReadout( diff --git a/tests/scripts/test_epi2d_se.py b/tests/scripts/test_epi2d_se.py index b716a5b..6d8daf7 100644 --- a/tests/scripts/test_epi2d_se.py +++ b/tests/scripts/test_epi2d_se.py @@ -1,7 +1,5 @@ """Tests for EPI SE sequence.""" -import tempfile -from pathlib import Path from typing import Literal import ismrmrd @@ -49,15 +47,14 @@ def test_mrd_trajectory( ramp_sampling: bool, oversampling: Literal[1, 2], n_slices: int, + tmp_path, ): """Test that the MRD trajectory matches the analytical and PyPulseq trajectories.""" fov = 200e-3 n_readout = 32 n_phase_encoding = 32 - mrd_file = Path(tempfile.gettempdir()) / 'test_epi2d_se_traj.h5' - if mrd_file.exists(): - mrd_file.unlink() + mrd_file = tmp_path / 'test_epi2d_se_traj.h5' seq, _, _ = epi2d_se_kernel( system=system_defaults, @@ -88,7 +85,6 @@ def test_mrd_trajectory( n_acq = ds.number_of_acquisitions() traj_mrd = [ds.read_acquisition(i).traj.copy() for i in range(n_acq)] ds.close() - mrd_file.unlink() # Recreate EpiReadout for reference epi2d = EpiReadout(