diff --git a/examples/scripts/write_2Dt1_mprage.py b/examples/scripts/write_2Dt1_mprage.py index 62aab9d..5cb05e0 100644 --- a/examples/scripts/write_2Dt1_mprage.py +++ b/examples/scripts/write_2Dt1_mprage.py @@ -39,10 +39,11 @@ apodization=0.5, time_bw_product=4, return_gz=True, + use='excitation', ) flip90 = 90 * pi / 180 -rf90 = pp.make_block_pulse(flip_angle=flip90, system=system, duration=500e-6, time_bw_product=4) +rf90 = pp.make_block_pulse(flip_angle=flip90, system=system, duration=500e-6, time_bw_product=4, use='preparation') # ========= # Readout diff --git a/examples/scripts/write_3Dt1_mprage.py b/examples/scripts/write_3Dt1_mprage.py index 392a692..60ddf7c 100644 --- a/examples/scripts/write_3Dt1_mprage.py +++ b/examples/scripts/write_3Dt1_mprage.py @@ -28,10 +28,12 @@ # RF preparatory, excitation # ========= flip_exc = 12 * pi / 180 -rf = pp.make_block_pulse(flip_angle=flip_exc, system=system, duration=250e-6, time_bw_product=4) +rf = pp.make_block_pulse(flip_angle=flip_exc, system=system, duration=250e-6, time_bw_product=4, use='excitation') flip_prep = 90 * pi / 180 -rf_prep = pp.make_block_pulse(flip_angle=flip_prep, system=system, duration=500e-6, time_bw_product=4) +rf_prep = pp.make_block_pulse( + flip_angle=flip_prep, system=system, duration=500e-6, time_bw_product=4, use='preparation' +) # ========= # Readout diff --git a/examples/scripts/write_epi.py b/examples/scripts/write_epi.py index 3f8425c..aec5ade 100644 --- a/examples/scripts/write_epi.py +++ b/examples/scripts/write_epi.py @@ -43,6 +43,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'epi_p time_bw_product=4, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Define other gradients and ADC events diff --git a/examples/scripts/write_epi_label.py b/examples/scripts/write_epi_label.py index dcc0ded..af99e4f 100644 --- a/examples/scripts/write_epi_label.py +++ b/examples/scripts/write_epi_label.py @@ -47,6 +47,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'epi_l time_bw_product=4, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Define trigger diff --git a/examples/scripts/write_epi_se.py b/examples/scripts/write_epi_se.py index 4b6bac1..7f033de 100644 --- a/examples/scripts/write_epi_se.py +++ b/examples/scripts/write_epi_se.py @@ -39,6 +39,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'epi_s time_bw_product=4, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Define other gradients and ADC events diff --git a/examples/scripts/write_epi_se_rs.py b/examples/scripts/write_epi_se_rs.py index 99289e5..c2f7a11 100644 --- a/examples/scripts/write_epi_se_rs.py +++ b/examples/scripts/write_epi_se_rs.py @@ -57,6 +57,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'epi_s bandwidth=np.abs(sat_freq), freq_offset=sat_freq, delay=system.rf_dead_time, + use='saturation', ) gz_fs = pp.make_trapezoid(channel='z', system=system, delay=pp.calc_duration(rf_fs), area=1 / 1e-4) @@ -70,6 +71,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'epi_s time_bw_product=4, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Create 90 degree slice refocusing pulse and gradients diff --git a/examples/scripts/write_gre.py b/examples/scripts/write_gre.py index 78dfa02..47500df 100644 --- a/examples/scripts/write_gre.py +++ b/examples/scripts/write_gre.py @@ -5,7 +5,7 @@ import pypulseq as pp -def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_pypulseq.seq'): +def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_pypulseq.seq', paper_plot: bool = False): # ====== # SETUP # ====== @@ -44,6 +44,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_p system=system, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Define other gradients and ADC events delta_k = 1 / fov @@ -120,7 +121,10 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_p # VISUALIZATION # ====== if plot: - seq.plot() + if paper_plot: + seq.paper_plot() + else: + seq.plot() seq.calculate_kspace() @@ -143,4 +147,4 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_p if __name__ == '__main__': - main(plot=False, write_seq=True) + seq = main(plot=True, paper_plot=True, write_seq=False) diff --git a/examples/scripts/write_gre_label.py b/examples/scripts/write_gre_label.py index bf06bd8..05da178 100644 --- a/examples/scripts/write_gre_label.py +++ b/examples/scripts/write_gre_label.py @@ -47,6 +47,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_l system=system, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Define other gradients and ADC events diff --git a/examples/scripts/write_haste.py b/examples/scripts/write_haste.py index 6359e9c..7bcd52d 100644 --- a/examples/scripts/write_haste.py +++ b/examples/scripts/write_haste.py @@ -72,6 +72,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'haste phase_offset=rfex_phase, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) GS_ex = make_trapezoid( channel='z', diff --git a/examples/scripts/write_mprage.py b/examples/scripts/write_mprage.py index 6476af1..fbedd32 100644 --- a/examples/scripts/write_mprage.py +++ b/examples/scripts/write_mprage.py @@ -45,9 +45,16 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'mprag ax.n3 = xyz.index(ax.d3) # Create alpha-degree hard pulse and gradient - rf = pp.make_block_pulse(flip_angle=alpha * np.pi / 180, system=system, duration=rf_len, delay=system.rf_dead_time) + rf = pp.make_block_pulse( + flip_angle=alpha * np.pi / 180, system=system, duration=rf_len, delay=system.rf_dead_time, use='excitation' + ) rf180 = pp.make_adiabatic_pulse( - pulse_type='hypsec', system=system, duration=10.24e-3, dwell=1e-5, delay=system.rf_dead_time + pulse_type='hypsec', + system=system, + duration=10.24e-3, + dwell=1e-5, + delay=system.rf_dead_time, + use='inversion', ) # Define other gradients and ADC events diff --git a/examples/scripts/write_radial_gre.py b/examples/scripts/write_radial_gre.py index 8109e33..0385bc7 100644 --- a/examples/scripts/write_radial_gre.py +++ b/examples/scripts/write_radial_gre.py @@ -45,6 +45,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'gre_r time_bw_product=4, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Define other gradients and ADC events diff --git a/examples/scripts/write_tse.py b/examples/scripts/write_tse.py index 8255dea..75abfac 100644 --- a/examples/scripts/write_tse.py +++ b/examples/scripts/write_tse.py @@ -63,6 +63,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'tse_p phase_offset=rf_ex_phase, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) gs_ex = pp.make_trapezoid( channel='z', diff --git a/examples/scripts/write_ute.py b/examples/scripts/write_ute.py index 4817df2..10fc293 100644 --- a/examples/scripts/write_ute.py +++ b/examples/scripts/write_ute.py @@ -54,6 +54,7 @@ def main(plot: bool = False, write_seq: bool = False, seq_filename: str = 'ute_p system=system, return_gz=True, delay=system.rf_dead_time, + use='excitation', ) # Align RO asymmetry to ADC samples diff --git a/pyproject.toml b/pyproject.toml index eda2dfd..ee5a606 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "pypulseq" -version = "1.4.2.post1" +version = "1.5.0" authors = [{ name = "Keerthi Sravan Ravi", email = "ks3621@columbia.edu" }] maintainers = [ { name = "Bilal Tasdelen" }, @@ -28,6 +28,7 @@ dependencies = [ [project.optional-dependencies] sigpy = ["sigpy>=0.1.26"] +mplcursors = ["mplcursors"] test = [ "coverage", "codecov", diff --git a/src/pypulseq/Sequence/block.py b/src/pypulseq/Sequence/block.py index 7330ff6..435894a 100644 --- a/src/pypulseq/Sequence/block.py +++ b/src/pypulseq/Sequence/block.py @@ -125,7 +125,7 @@ def set_block(self, block_index: int, *args: Union[SimpleNamespace, float]) -> N if hasattr(event, 'id'): adc_id = event.id else: - adc_id = register_adc_event(self, event) + adc_id, _ = register_adc_event(self, event) new_block[5] = adc_id duration = max(duration, event.delay + event.num_samples * event.dwell + event.dead_time) @@ -216,6 +216,7 @@ def set_block(self, block_index: int, *args: Union[SimpleNamespace, float]) -> N # Now we add the ID new_block[6] = extension_id + # TODO: replace with independent function (#PR 289) # ========= # PERFORM GRADIENT CHECKS # ========= @@ -254,7 +255,7 @@ def set_block(self, block_index: int, *args: Union[SimpleNamespace, float]) -> N if prev_type == 't': last = 0 elif prev_type == 'g': - last = prev_lib['data'][5] + last = prev_lib['data'][2] # v150 # Check whether the difference between the last gradient value and # the first value of the new gradient is achievable with the @@ -273,7 +274,7 @@ def set_block(self, block_index: int, *args: Union[SimpleNamespace, float]) -> N if next_type == 't': first = 0 elif next_type == 'g': - first = next_lib['data'][4] + first = next_lib['data'][1] # v150 else: first = 0 @@ -288,15 +289,63 @@ def set_block(self, block_index: int, *args: Union[SimpleNamespace, float]) -> N raise RuntimeError('First gradient in the the first block has to start at 0.') if ( - grad_to_check.stop[1] > self.system.max_slew * self.system.grad_raster_time + abs(grad_to_check.stop[1]) > self.system.max_slew * self.system.grad_raster_time and abs(grad_to_check.stop[0] - duration) > 1e-7 ): raise RuntimeError("A gradient that doesn't end at zero needs to be aligned to the block boundary.") + # ========= + # END GRADIENT CHECKS + # ========= self.block_events[block_index] = new_block self.block_durations[block_index] = float(duration) +def get_raw_block_content_IDs(self, block_index: int) -> SimpleNamespace: + """ + Returns PyPulseq block content IDs at `block_index` position in `self.block_events`. + + No block events are created, only the IDs of the objects are returned. + + Parameters + ---------- + block_index : int + Index of PyPulseq block to be retrieved from `self.block_events`. + + Returns + ------- + block : SimpleNamespace + PyPulseq block content IDs at 'block_index' position in `self.block_events`. + """ + raw_block = SimpleNamespace(block_duration=0, rf=0, gx=0, gy=0, gz=0, adc=0, ext=[]) + event_ind = self.block_events[block_index] + + # Extensions + if event_ind[6] > 0: + next_ext_id = event_ind[6] + while next_ext_id != 0: + ext_data = self.extensions_library.data[next_ext_id] + raw_block.ext.append(ext_data[:2]) + next_ext_id = ext_data[2] + raw_block.ext = np.stack(raw_block.ext, axis=-1) + + # RF + if event_ind[1] > 0: + raw_block.rf = event_ind[1] + + # Gradients + grad_channels = ['gx', 'gy', 'gz'] + for i in range(len(grad_channels)): + if event_ind[2 + i] > 0: + setattr(raw_block, grad_channels[i], event_ind[2 + i]) + + # ADC + if event_ind[5] > 0: + raw_block.adc = event_ind[5] + + return raw_block + + def get_block(self, block_index: int) -> SimpleNamespace: """ Returns PyPulseq block at `block_index` position in `self.block_events`. @@ -342,6 +391,8 @@ def get_block(self, block_index: int) -> SimpleNamespace: else: block.rf = self.rf_from_lib_data(self.rf_library.data[event_ind[1]], 'u') # Undefined type/use + # TODO: add optional rf ID from raw_block + # Gradients grad_channels = ['gx', 'gy', 'gz'] for i in range(len(grad_channels)): @@ -353,9 +404,9 @@ def get_block(self, block_index: int) -> SimpleNamespace: grad.channel = grad_channels[i][1] if grad.type == 'grad': amplitude = lib_data[0] - shape_id = lib_data[1] - time_id = lib_data[2] - delay = lib_data[3] + shape_id = lib_data[3] # change in v150 + time_id = lib_data[4] # change in v150 + delay = lib_data[5] # change in v150 shape_data = self.shape_library.data[shape_id] compressed.num_samples = shape_data[0] compressed.data = shape_data[1:] @@ -365,22 +416,36 @@ def get_block(self, block_index: int) -> SimpleNamespace: if time_id == 0: grad.tt = (np.arange(1, len(g) + 1) - 0.5) * self.grad_raster_time t_end = len(g) * self.grad_raster_time + grad.area = sum(grad.waveform) * self.grad_raster_time + elif time_id == -1: + # Gradient with oversampling by a factor of 2 + grad.tt = 0.5 * (np.arange(1, len(g) + 1)) * self.grad_raster_time + if len(grad.tt) != len(grad.waveform): + raise ValueError( + f'Mismatch between time shape length ({len(grad.tt)}) and gradient shape length ({len(grad.waveform)}).' + ) + if len(grad.waveform) % 2 != 1: + raise ValueError('Oversampled gradient waveforms must have odd number of samples') + t_end = (len(g) + 1) * self.grad_raster_time + grad.area = sum(grad.waveform[::2]) * self.grad_raster_time # remove oversampling else: t_shape_data = self.shape_library.data[time_id] compressed.num_samples = t_shape_data[0] compressed.data = t_shape_data[1:] grad.tt = decompress_shape(compressed) * self.grad_raster_time - - assert len(grad.waveform) == len(grad.tt) + if len(grad.tt) != len(grad.waveform): + raise ValueError( + f'Mismatch between time shape length ({len(grad.tt)}) and gradient shape length ({len(grad.waveform)}).' + ) t_end = grad.tt[-1] + grad.area = 0.5 * sum((grad.tt[1:] - grad.tt[:-1]) * (grad.waveform[1:] + grad.waveform[:-1])) grad.shape_id = shape_id grad.time_id = time_id grad.delay = delay grad.shape_dur = t_end - if len(lib_data) > 5: - grad.first = lib_data[4] - grad.last = lib_data[5] + grad.first = lib_data[1] # change in v150 - we always have first/last now + grad.last = lib_data[2] # change in v150 - we always have first/last now else: grad.amplitude = lib_data[0] grad.rise_time = lib_data[1] @@ -390,23 +455,38 @@ def get_block(self, block_index: int) -> SimpleNamespace: grad.area = grad.amplitude * (grad.flat_time + grad.rise_time / 2 + grad.fall_time / 2) grad.flat_area = grad.amplitude * grad.flat_time + # TODO: add optional grad ID from raw_block + setattr(block, grad_channels[i], grad) # ADC if event_ind[5] > 0: lib_data = self.adc_library.data[event_ind[5]] + shape_id_phase_modulation = lib_data[-2] + if shape_id_phase_modulation: + shape_data = self.shape_library.data[shape_id_phase_modulation] + compressed = SimpleNamespace() + compressed.num_samples = shape_data[0] + compressed.data = shape_data[1:] + phase_shape = decompress_shape(compressed) + else: + phase_shape = np.array([], dtype=float) adc = SimpleNamespace() - ( - adc.num_samples, - adc.dwell, - adc.delay, - adc.freq_offset, - adc.phase_offset, - adc.dead_time, - ) = [lib_data[x] for x in range(6)] + adc.num_samples = lib_data[0] + adc.dwell = lib_data[1] + adc.delay = lib_data[2] + adc.freq_ppm = lib_data[3] + adc.phase_ppm = lib_data[4] + adc.freq_offset = lib_data[5] + adc.phase_offset = lib_data[6] + adc.phase_modulation = phase_shape + adc.dead_time = self.system.adc_dead_time adc.num_samples = int(adc.num_samples) adc.type = 'adc' + + # TODO: add optional adc ID from raw_block + block.adc = adc if event_ind[6] > 0: @@ -486,7 +566,7 @@ def get_block(self, block_index: int) -> SimpleNamespace: return block -def register_adc_event(self, event: EventLibrary) -> int: +def register_adc_event(self, event: EventLibrary) -> Tuple[int, List[int]]: """ Parameters @@ -496,25 +576,52 @@ def register_adc_event(self, event: EventLibrary) -> int: Returns ------- - int - ID of registered ADC event. + int, [int, ...] + ID of registered RF event, list of shape IDs """ + surely_new = False + + # Handle phase modulation + if not hasattr(event, 'phase_modulation') or event.phase_modulation is None or len(event.phase_modulation) == 0: + shape_id = 0 + else: + if hasattr(event, 'shape_id'): + shape_id = event.shape_id + else: + phase_shape = compress_shape(np.asarray(event.phase_modulation).flatten()) + shape_data = np.concatenate(([phase_shape.num_samples], phase_shape.data)) + shape_id, shape_found = self.shape_library.find_or_insert(shape_data) + if not shape_found: + surely_new = True + + # Construct the ADC event data data = ( event.num_samples, event.dwell, - event.delay, + max(event.delay, event.dead_time), + event.freq_ppm, + event.phase_ppm, event.freq_offset, event.phase_offset, + shape_id, event.dead_time, ) - adc_id, found = self.adc_library.find_or_insert(new_data=data) - # Clear block cache because ADC was overwritten - # TODO: Could find only the blocks that are affected by the changes - if self.use_block_cache and found: - self.block_cache.clear() + # Insert or find/insert into libraryAdd commentMore actions + if surely_new: + adc_id = self.adc_library.insert(0, data) + else: + adc_id, found = self.adc_library.find_or_insert(data) + + # Clear block cache if overwritten + if self.use_block_cache and found: + self.block_cache.clear() + + # Optional mapping + if hasattr(event, 'name'): + self.adc_id_to_name_map[adc_id] = event.name - return adc_id + return adc_id, shape_id def register_control_event(self, event: SimpleNamespace) -> int: @@ -567,37 +674,45 @@ def register_grad_event(self, event: SimpleNamespace) -> Union[int, Tuple[int, L """ may_exist = True any_changed = False + if event.type == 'grad': - amplitude = np.abs(event.waveform).max() + amplitude = np.max(np.abs(event.waveform)) if amplitude > 0: fnz = event.waveform[np.nonzero(event.waveform)[0][0]] - amplitude *= np.sign(fnz) if fnz != 0 else 1 # Workaround for np.sign(0) = 0 + amplitude *= np.sign(fnz) if fnz != 0 else 1 + # Shape ID initialization if hasattr(event, 'shape_IDs'): shape_IDs = event.shape_IDs else: shape_IDs = [0, 0] - if amplitude != 0: - g = event.waveform / amplitude - else: - g = event.waveform + + # Shape for waveform + g = event.waveform / amplitude if amplitude != 0 else event.waveform c_shape = compress_shape(g) s_data = np.concatenate(([c_shape.num_samples], c_shape.data)) shape_IDs[0], found = self.shape_library.find_or_insert(s_data) - may_exist = may_exist & found + may_exist = may_exist and found any_changed = any_changed or found - # Check whether tt == np.arange(len(event.tt)) * self.grad_raster_time + 0.5 - tt_regular = (np.floor(event.tt / self.grad_raster_time) == np.arange(len(event.tt))).all() + # Shape for timing + c_time = compress_shape(event.tt / self.grad_raster_time) + t_data = np.concatenate(([c_time.num_samples], c_time.data)) - if not tt_regular: - c_time = compress_shape(event.tt / self.grad_raster_time) - t_data = np.concatenate(([c_time.num_samples], c_time.data)) + if len(c_time.data) == 4 and np.allclose(c_time.data, [0.5, 1, 1, c_time.num_samples - 3]): + # Standard raster → leave shape_IDs[1] as 0 + pass + elif len(c_time.data) == 3 and np.allclose(c_time.data, [0.5, 0.5, c_time.num_samples - 2]): + # Half-raster → set to -1 as special flag + shape_IDs[1] = -1 + else: shape_IDs[1], found = self.shape_library.find_or_insert(t_data) - may_exist = may_exist & found + may_exist = may_exist and found any_changed = any_changed or found - data = (amplitude, *shape_IDs, event.delay, event.first, event.last) + # Updated data layout to match MATLAB v1.5.0 ordering + data = (amplitude, event.first, event.last, *shape_IDs, event.delay) + elif event.type == 'trap': data = ( event.amplitude, @@ -620,6 +735,9 @@ def register_grad_event(self, event: SimpleNamespace) -> Union[int, Tuple[int, L if self.use_block_cache and any_changed: self.block_cache.clear() + if hasattr(event, 'name'): + self.grad_id_to_name_map[grad_id] = event.name + if event.type == 'grad': return grad_id, shape_IDs elif event.type == 'trap': @@ -760,8 +878,19 @@ def register_rf_event(self, event: SimpleNamespace) -> Tuple[int, List[int]]: use = event.use[0] else: use = 'u' + else: + raise ValueError('Parameter "use" is not optional since v1.5.0') - data = (amplitude, *shape_IDs, event.delay, event.freq_offset, event.phase_offset) + data = ( + amplitude, + *shape_IDs, + event.center, + event.delay, + event.freq_ppm, + event.phase_ppm, + event.freq_offset, + event.phase_offset, + ) if may_exist: rf_id, found = self.rf_library.find_or_insert(new_data=data, data_type=use) @@ -773,4 +902,7 @@ def register_rf_event(self, event: SimpleNamespace) -> Tuple[int, List[int]]: else: rf_id = self.rf_library.insert(key_id=0, new_data=data, data_type=use) + if hasattr(event, 'name'): + self.rf_id_to_name_map[rf_id] = event.name + return rf_id, shape_IDs diff --git a/src/pypulseq/Sequence/read_seq.py b/src/pypulseq/Sequence/read_seq.py index e508fd4..efaba88 100644 --- a/src/pypulseq/Sequence/read_seq.py +++ b/src/pypulseq/Sequence/read_seq.py @@ -1,11 +1,12 @@ import re import warnings from types import SimpleNamespace -from typing import Dict, List, Tuple +from typing import Dict, List, Optional, Tuple import numpy as np from pypulseq.calc_duration import calc_duration +from pypulseq.calc_rf_center import calc_rf_center from pypulseq.compress_shape import compress_shape from pypulseq.decompress_shape import decompress_shape from pypulseq.event_lib import EventLibrary @@ -118,6 +119,11 @@ def read(self, path: str, detect_rf_use: bool = False, remove_duplicates: bool = f'{version_major}.{version_minor}.{version_revision}) some code may function not as ' f'expected' ) + + if version_combined < 1005000 and detect_rf_use: + warnings.warn('Option detectRFuse is not supported for file format version 1.5.0 and above') + detect_rf_use = False + elif section == '[BLOCKS]': if version_major == 0: raise RuntimeError('Pulseq file MUST include [VERSION] section prior to [BLOCKS] section') @@ -131,7 +137,13 @@ def read(self, path: str, detect_rf_use: bool = False, remove_duplicates: bool = if jemris_generated: self.rf_library = __read_events(input_file, (1, 1, 1, 1, 1), event_library=self.rf_library) else: - if version_combined >= 1004000: # 1.4.x format + if version_combined >= 1005000: # 1.5.x format + self.rf_library = __read_events( + input_file, + (1, 1, 1, 1, 1e-6, 1e-6, 1, 1, 1, 1, np.nan), + event_library=self.rf_library, + ) + elif version_combined >= 1004000: # 1.4.x format self.rf_library = __read_events( input_file, (1, 1, 1, 1, 1e-6, 1, 1), @@ -140,7 +152,9 @@ def read(self, path: str, detect_rf_use: bool = False, remove_duplicates: bool = else: # 1.3.x and below self.rf_library = __read_events(input_file, (1, 1, 1, 1e-6, 1, 1), event_library=self.rf_library) elif section == '[GRADIENTS]': - if version_combined >= 1004000: # 1.4.x format + if version_combined >= 1005000: # 1.5.x format + self.grad_library = __read_events(input_file, (1, 1, 1, 1, 1, 1e-6), 'g', self.grad_library) + elif version_combined >= 1004000: # 1.4.x format self.grad_library = __read_events(input_file, (1, 1, 1, 1e-6), 'g', self.grad_library) else: # 1.3.x and below self.grad_library = __read_events(input_file, (1, 1, 1e-6), 'g', self.grad_library) @@ -150,9 +164,17 @@ def read(self, path: str, detect_rf_use: bool = False, remove_duplicates: bool = else: self.grad_library = __read_events(input_file, (1, 1e-6, 1e-6, 1e-6, 1e-6), 't', self.grad_library) elif section == '[ADC]': - self.adc_library = __read_events( - input_file, (1, 1e-9, 1e-6, 1, 1), event_library=self.adc_library, append=self.system.adc_dead_time - ) + if version_combined >= 1005000: # 1.5.x format + self.adc_library = __read_events( + input_file, + (1, 1e-9, 1e-6, 1, 1, 1, 1, 1), + event_library=self.adc_library, + append=self.system.adc_dead_time, + ) + else: # 1.4.x format and below + self.adc_library = __read_events( + input_file, (1, 1e-9, 1e-6, 1, 1), event_library=self.adc_library, append=self.system.adc_dead_time + ) elif section == '[DELAYS]': if version_combined >= 1004000: raise RuntimeError('Pulseq file revision 1.4.0 and above MUST NOT contain [DELAYS] section') @@ -209,11 +231,19 @@ def l2(s): f'are supported.' ) - # Fix blocks, gradients and RF objects imported from older versions + # Fix blocks, gradients and RF objects imported from older versions (< v1.4.0) if version_combined < 1004000: # Scan through RF objects + self.rf_library.type = dict.fromkeys(self.rf_library.type.keys(), 'u') for i in self.rf_library.data: - self.rf_library.update(i, None, (*self.rf_library.data[i][:3], 0, *self.rf_library.data[i][3:])) + d = self.rf_library.data[i] + rf = self.rf_from_lib_data((d[:3], 0, 0, d[3], 0, 0, d[4:6], 'u')).__delattr__('center') + center = calc_rf_center(rf) + self.rf_library.update( + i, + None, + (d[:3], 0, center, d[3], 0, 0, d[4:6], 'u'), + ) # 0 between [3] and [4:6] are the freq_ppm and phase_ppm # Scan through the gradient objects and update 't'-s (trapezoids) und 'g'-s (free-shape gradients) for i in self.grad_library.data: @@ -263,78 +293,112 @@ def l2(s): # Calculate actual block duration self.block_durations[block_counter] = calc_duration(block) - # TODO: Is it possible to avoid expensive get_block calls here? - grad_channels = ['gx', 'gy', 'gz'] - grad_prev_last = np.zeros(len(grad_channels)) - for block_counter in self.block_events: - block = self.get_block(block_counter) - block_duration = block.block_duration - # We also need to keep track of the event IDs because some PyPulseq files written by external software may contain - # repeated entries so searching by content will fail - event_idx = self.block_events[block_counter] - # Update the objects by filling in the 'first' and 'last' attributes not yet contained in the Pulseq file - for j in range(len(grad_channels)): - grad = getattr(block, grad_channels[j]) - if grad is None: - grad_prev_last[j] = 0 - continue - - if grad.type == 'grad': - if grad.delay > 0: - grad_prev_last[j] = 0 + elif version_combined < 1005000: + # Port from v1.4.x : RF, ADC and GRAD objects need to be updated + # this needs to be done on the level of the libraries, because get_block will fail - # go to next channel, if grad.first and grad.last are already set - if hasattr(grad, 'first') and hasattr(grad, 'last'): - grad_prev_last[j] = grad.last - continue + # Scan though the RFs and add center, freq_ppm, phase_ppm and use fields + self.rf_library.type = dict.fromkeys(self.rf_library.type.keys(), 'u') + for i in self.rf_library.data: + # Use goes into the type field, and this is done separately + d = self.rf_library.data[i] + rf = self.rf_from_lib_data((d[:4], 0, d[4], 0, 0, d[5:7], 'u')).__delattr__('center') + center = calc_rf_center(rf) + self.rf_library.update( + i, + None, + (d[:4], center, d[4], 0, 0, d[5:7], 'u'), + ) # 0 between [4] and [5:7] are the freqPPM and phasePPM + + # Scan through the gradient objects and update 'g'-s (free-shape gradients) + for i in self.grad_library.data: + if self.grad_library.type[i] == 'g': + self.grad_library.update( + i, + None, + ( + self.grad_library.data[i][0], + None, + None, + self.grad_library.data[i][1:4], + ), + self.grad_library.type[i], + ) # We use None to label the non-initialized first/last fields. These will be restored in the code below - # get grad.first and grad.last attributes from the grad_library if they have been set for the current amplitude_ID before - amplitude_ID = event_idx[j + 2] - if amplitude_ID in event_idx[2 : (j + 2)]: - if self.use_block_cache: - grad.first = self.grad_library.data[amplitude_ID][4] - grad.last = self.grad_library.data[amplitude_ID][5] + # Another run through for all older versions + if version_combined < 1005000: + # TODO: Is it possible to avoid expensive get_block calls here? + grad_channels = ['gx', 'gy', 'gz'] + grad_prev_last = np.zeros(len(grad_channels)) + for block_counter in self.block_events: + block = self.get_block(block_counter) + block_duration = block.block_duration + # We also need to keep track of the event IDs because some PyPulseq files written by external software may contain + # repeated entries so searching by content will fail + event_idx = self.block_events[block_counter] + # Update the objects by filling in the 'first' and 'last' attributes not yet contained in the Pulseq file + for j in range(len(grad_channels)): + grad = getattr(block, grad_channels[j]) + if grad is None: + grad_prev_last[j] = 0 continue - # get time_id from grad_library - time_id = self.grad_library.data[amplitude_ID][2] - - # if grad.first is not set, set it to the last value of the previous gradient - grad.first = grad_prev_last[j] + if grad.type == 'grad': + if grad.delay > 0: + grad_prev_last[j] = 0 + + # go to next channel, if grad.first and grad.last are already set + if hasattr(grad, 'first') and hasattr(grad, 'last'): + grad_prev_last[j] = grad.last + continue + + # get grad.first and grad.last attributes from the grad_library if they have been set for the current amplitude_ID before + amplitude_ID = event_idx[j + 2] + if amplitude_ID in event_idx[2 : (j + 2)]: + if self.use_block_cache: + grad.first = self.grad_library.data[amplitude_ID][4] + grad.last = self.grad_library.data[amplitude_ID][5] + continue + + # get time_id from grad_library + time_id = self.grad_library.data[amplitude_ID][2] + + # if grad.first is not set, set it to the last value of the previous gradient + grad.first = grad_prev_last[j] + + # extended trapezoid: use last value of the gradient waveform as grad.last + if time_id != 0: + grad.last = grad.waveform[-1] + grad_duration = grad.delay + grad.tt[-1] + # arbitrary gradients: interpolate grad.last from the gradient waveform + else: + # use a linear extrapolation identical to the one used in the make_arbitrary_grad.py file + grad.last = (3 * grad.waveform[-1] - grad.waveform[-2]) * 0.5 + grad_duration = grad.delay + len(grad.waveform) * self.grad_raster_time + + # Set grad_prev_last to 0 if gradient does not end at block boundary + eps = np.finfo(np.float64).eps + if grad_duration + eps < block_duration: + grad_prev_last[j] = 0 + # Update grad_prev_last for the next iteration if gradient ends at block boundary + else: + grad_prev_last[j] = grad.last + + # Update the grad_library with the new grad.first and grad.last values + amplitude = self.grad_library.data[amplitude_ID][0] + shape_id = self.grad_library.data[amplitude_ID][1] + new_data = ( + amplitude, + shape_id, + time_id, + grad.delay, + grad.first, + grad.last, + ) + self.grad_library.update_data(amplitude_ID, None, new_data, 'g') - # extended trapezoid: use last value of the gradient waveform as grad.last - if time_id != 0: - grad.last = grad.waveform[-1] - grad_duration = grad.delay + grad.tt[-1] - # arbitrary gradients: interpolate grad.last from the gradient waveform else: - # use a linear extrapolation identical to the one used in the make_arbitrary_grad.py file - grad.last = (3 * grad.waveform[-1] - grad.waveform[-2]) * 0.5 - grad_duration = grad.delay + len(grad.waveform) * self.grad_raster_time - - # Set grad_prev_last to 0 if gradient does not end at block boundary - eps = np.finfo(np.float64).eps - if grad_duration + eps < block_duration: grad_prev_last[j] = 0 - # Update grad_prev_last for the next iteration if gradient ends at block boundary - else: - grad_prev_last[j] = grad.last - - # Update the grad_library with the new grad.first and grad.last values - amplitude = self.grad_library.data[amplitude_ID][0] - shape_id = self.grad_library.data[amplitude_ID][1] - new_data = ( - amplitude, - shape_id, - time_id, - grad.delay, - grad.first, - grad.last, - ) - self.grad_library.update_data(amplitude_ID, None, new_data, 'g') - - else: - grad_prev_last[j] = 0 if detect_rf_use: # Find the RF pulses, list flip angles, and work around the current (rev 1.2.0) Pulseq file format limitation @@ -473,7 +537,7 @@ def __read_blocks( def __read_events( - input_file, scale: tuple = (1,), event_type: str = str(), event_library: EventLibrary = None, append=None + input_file, scale: Optional[Tuple] = None, event_type: str = str(), event_library: EventLibrary = None, append=None ) -> EventLibrary: """ Read an event section of a sequence file and return a library of events. @@ -496,18 +560,45 @@ def __read_events( """ if event_library is None: event_library = EventLibrary() - line = __strip_line(input_file) + + # New in v1.5.0 : generate format string; NaN labels character param(s) for 'use' attribute + line, scale, format_spec = __read_format(input_file, scale) + data_mask = np.isfinite(scale) + type_idx = np.where(np.logical_not(data_mask))[0] + + if len(type_idx) > 1: + raise ValueError('Only one type field (marked as NaN) can be provided in scale.') + if len(type_idx) == 0: + type_idx = None + data_mask = None + else: + type_idx = type_idx.item() while line != '' and line != '#': - data = np.fromstring(line, dtype=float, sep=' ') + data = __fromstring(line, format_spec) event_id = data[0] - data = tuple(data[1:] * scale) + + if type_idx is not None: + event_type = data[type_idx + 1] # Need +1 because of the event_id in the first position + + data = data[1:] + data = tuple(data[n] * scale[n] if isinstance(data[n], str) is False else data[n] for n in range(len(data))) + if append is not None: data = (*data, append) - if event_type == '': - event_library.insert(key_id=event_id, new_data=data) + + if event_type == '' and type_idx is None: + if data_mask is None: + event_library.insert(key_id=event_id, new_data=data) + else: + event_library.insert(key_id=event_id, new_data=data[data_mask]) else: - event_library.insert(key_id=event_id, new_data=data, data_type=event_type) + if data_mask is None: + event_library.insert(key_id=event_id, new_data=data, data_type=event_type) + else: + data = tuple(np.asarray(data, dtype=object)[data_mask]) + event_library.insert(key_id=event_id, new_data=data, data_type=event_type) + line = __strip_line(input_file) return event_library @@ -641,3 +732,76 @@ def __strip_line(input_file) -> str: """ line = input_file.readline() # If line is an empty string, end of the file has been reached return line.strip() if line != '' else -1 + + +def __read_format(input_file, scale: Tuple) -> Tuple[List, Tuple, List]: + """ + Generate a format specifier list based on the scale vector. + + '%f' for numeric values, '%s' for string (NaN in scale). + + Parameters + ---------- + input_file: file + Input text + scale : list, default=(1,) + Scale elements according to column vector scale. + + Returns + ------- + line : str + First line in input_file after spaces and newline whitespaces have been removed. Note: File pointer is + remembered, and hence successive calls work as expected. Returns -1 for eof. + scale : tuple[float] + Scaling factor for each element in input line. + Defaults to one for each numeric element. + format : list[str] + List of format tokens for each field (excluding the event ID). + + """ + line = __strip_line(input_file) + + if scale is None: + tok = line.strip().split() + scale = (len(tok) - 1) * (1,) + + scale = np.array(scale, dtype=float) + is_num = np.isfinite(scale) + format_spec = ['%f' if value else '%s' for value in is_num] + + return line, scale, format_spec + + +def __fromstring(line, format_spec: List) -> List: + """ + Parse a line of text using a dynamic format specification. + + Parameters + ---------- + line : str + Input line (e.g., ['23', '1.0', '2.0', '3.0', 'u']) + format_spec : list[str] + Format list like ['%f', '%f', '%f', '%f', '%s'] + + Returns + ------- + data : list + Parsed data. + + """ + tok = line.strip().split() + if len(tok) != len(format_spec) + 1: + raise ValueError('Mismatch between number of tokens and format spec') + + format_spec = ['%f', *format_spec] + data = [] + + for i, fmt in enumerate(format_spec): + if fmt == '%f': + data.append(float(tok[i])) + elif fmt == '%s': + data.append(tok[i]) + else: + raise ValueError(f'Unsupported format: {fmt}') + + return data diff --git a/src/pypulseq/Sequence/sequence.py b/src/pypulseq/Sequence/sequence.py index 2491e3a..003e42d 100644 --- a/src/pypulseq/Sequence/sequence.py +++ b/src/pypulseq/Sequence/sequence.py @@ -1,4 +1,3 @@ -import itertools import math from collections import OrderedDict from copy import deepcopy @@ -13,9 +12,7 @@ Self = TypeVar('Self', bound='Sequence') -import matplotlib as mpl import numpy as np -from matplotlib import pyplot as plt from scipy.interpolate import PPoly from pypulseq import __version__, eps @@ -25,7 +22,7 @@ from pypulseq.decompress_shape import decompress_shape from pypulseq.event_lib import EventLibrary from pypulseq.opts import Opts -from pypulseq.Sequence import block, parula +from pypulseq.Sequence import block from pypulseq.Sequence.calc_grad_spectrum import calculate_gradient_spectrum from pypulseq.Sequence.calc_pns import calc_pns from pypulseq.Sequence.ext_test_report import ext_test_report @@ -33,8 +30,9 @@ from pypulseq.Sequence.read_seq import read from pypulseq.Sequence.write_seq import write as write_seq from pypulseq.Sequence.write_seq import write_v141 as write_seq_v141 -from pypulseq.supported_labels_rf_use import get_supported_labels from pypulseq.utils.cumsum import cumsum +from pypulseq.utils.paper_plot import paper_plot as ext_paper_plot +from pypulseq.utils.seq_plot import SeqPlot from pypulseq.utils.tracing import format_trace, trace, trace_enabled major, minor, revision = __version__.split('.')[:3] @@ -94,6 +92,9 @@ def __init__(self, system: Union[Opts, None] = None, use_block_cache: bool = Tru self.signature_type = '' self.signature_file = '' self.signature_value = '' + self.rf_id_to_name_map = {} + self.adc_id_to_name_map = {} + self.grad_id_to_name_map = {} self.block_durations = {} self.extension_numeric_idx = [] @@ -615,6 +616,33 @@ def evaluate_labels(self, init: Union[dict, None] = None, evolution: str = 'none return labels + import numpy as np + + def find_block_by_time(self, t: float) -> int: + """ + Find the index of the block containing time `t`. + + Parameters + ---------- + t : float + Time (in seconds) to locate within the sequence. + + Returns + ------- + int or None + Index of the block that contains the given time, or None if out of range. + """ + cumsum_durations = np.cumsum(list(self.block_durations.values())) + block_index = np.searchsorted(cumsum_durations, t, side='right').item() + + if block_index >= len(self.block_durations): + return None + + if self.block_durations[block_index] <= 0: + raise ValueError('Block duration cannot be negative') + + return block_index + def flip_grad_axis(self, axis: str) -> None: """ Invert all gradients along the corresponding axis/channel. The function acts on all gradient objects already @@ -627,6 +655,28 @@ def flip_grad_axis(self, axis: str) -> None: """ self.mod_grad_axis(axis, modifier=-1) + def get_raw_block_content_IDs(self, block_index: int) -> SimpleNamespace: + """ + Returns PyPulseq block content IDs at `block_index` position in `self.block_events`. + + No block events are created, only the IDs of the objects are returned. + + See Also + -------- + - `pypulseq.Sequence.sequence.Sequence.get_block()`. + + Parameters + ---------- + block_index : int + Index of block to be retrieved from `Sequence`. + + Returns + ------- + SimpleNamespace + PyPulseq block content IDs at 'block_index' position in `self.block_events`. + """ + return block.get_raw_block_content_IDs(self, block_index) + def get_block(self, block_index: int) -> SimpleNamespace: """ Return a block of the sequence specified by the index. The block is created from the sequence data with all @@ -942,6 +992,43 @@ def mod_grad_axis(self, axis: str, modifier: float) -> None: if self.use_block_cache: self.block_cache.clear() + def paper_plot( + self, + time_range: Tuple[float] = (0, np.inf), + line_width: float = 1.2, + axes_color: Tuple[float] = (0.5, 0.5, 0.5), + rf_color: str = 'black', + gx_color: str = 'blue', + gy_color: str = 'red', + gz_color: Tuple[float] = (0, 0.5, 0.3), + rf_plot: str = 'abs', + ): + """ + Plot sequence using paper-style formatting (minimalist, high-contrast layout). + + Parameters + ---------- + time_range : iterable, default=(0, np.inf) + Time range (x-axis limits) for plotting the sequence. + Default is 0 to infinity (entire sequence). + line_width : float, default=1.2 + Line width used in plots. + axes_color : color, default=(0.5, 0.5, 0.5) + Color of horizontal zero axes (e.g., gray). + rf_color : color, default='black' + Color for RF and ADC events. + gx_color : color, default='blue' + Color for gradient X waveform. + gy_color : color, default='red' + Color for gradient Y waveform. + gz_color : color, default=(0, 0.5, 0.3) + Color for gradient Z waveform. + rf_plot : {'abs', 'real', 'imag'}, default='abs' + Determines how to plot RF waveforms (magnitude, real or imaginary part). + + """ + ext_paper_plot(self, time_range, line_width, axes_color, rf_color, gx_color, gy_color, gz_color, rf_plot) + def plot( self, label: str = str(), @@ -951,10 +1038,7 @@ def plot( time_disp: str = 's', grad_disp: str = 'kHz/m', plot_now: bool = True, - clear: bool = True, - fig1: Union[plt.Figure, None] = None, - fig2: Union[plt.Figure, None] = None, - ) -> Tuple[plt.Figure, Tuple[plt.Axes, plt.Axes, plt.Axes], plt.Figure, Tuple[plt.Axes, plt.Axes, plt.Axes]]: + ) -> SeqPlot: """ Plot `Sequence`. @@ -975,288 +1059,17 @@ def plot( grad_disp : str, default='kHz/m' Gradient display unit, must be one of 'kHz/m' or 'mT/m'. plot_now : bool, default=True - If True, function immediately shows the plots, blocking the rest of the code until plots are exited. - If False, plots are shown when plt.show() is called. Useful if plots are to be modified. - clear : bool, default=True - If True, clear existing figures before plotting (default behavior). - If False, overlay on existing figures 1 and 2 for sequence comparison. - fig1 : Optional[plt.Figure], default=None - Existing figure to plot RF/ADC events on. If None, a new figure is created. - fig2 : Optional[plt.Figure], default=None - Existing figure to plot gradients on. If None, a new figure is created. + If true, function immediately shows the plots, blocking the rest of the code until plots are exited. + If false, plots are shown when plt.show() is called. Useful if plots are to be modified. + plot_type : str, default='Gradient' + Gradients display type, must be one of either 'Gradient' or 'Kspace'. Returns ------- - Tuple[plt.Figure, Tuple[plt.Axes, plt.Axes, plt.Axes], plt.Figure, Tuple[plt.Axes, plt.Axes, plt.Axes]] - Returns (fig1, (sp11, sp12, sp13), fig2, (sp21, sp22, sp23)) for plot customization. - Always returns figures and axes regardless of plot_now setting. + SeqPlot + SeqPlot handle. """ - mpl.rcParams['lines.linewidth'] = 0.75 # Set default Matplotlib linewidth - - valid_time_units = ['s', 'ms', 'us'] - valid_grad_units = ['kHz/m', 'mT/m'] - valid_labels = get_supported_labels() - if not all(isinstance(x, (int, float)) for x in time_range) or len(time_range) != 2: - raise ValueError('Invalid time range') - if time_disp not in valid_time_units: - raise ValueError('Unsupported time unit') - - if grad_disp not in valid_grad_units: - raise ValueError('Unsupported gradient unit. Supported gradient units are: ' + str(valid_grad_units)) - - # Create the two figures (#1 for RF/ADC, #2 for gradients in x, y, z) or reuse existing figures - fig1 = plt.figure() if fig1 is None else fig1 - fig2 = plt.figure() if fig2 is None else fig2 - - # Clear existing figures if clear=True - if clear: - fig1.clear() - fig2.clear() - - # Create or reuse subplots of fig1 - fig1_axes = fig1.get_axes() - if not fig1_axes or clear: - sp11 = fig1.add_subplot(311) - sp12 = fig1.add_subplot(312, sharex=sp11) - sp13 = fig1.add_subplot(313, sharex=sp11) - else: - sp11, sp12, sp13 = fig1_axes[:3] - - # Create or reuse subplots of fig2 - fig2_axes = fig2.get_axes() - if not fig2_axes or clear: - sp21 = fig2.add_subplot(311, sharex=sp11) - sp22 = fig2.add_subplot(312, sharex=sp11) - sp23 = fig2.add_subplot(313, sharex=sp11) - else: - sp21, sp22, sp23 = fig2_axes[:3] - - t_factor_list = [1, 1e3, 1e6] - t_factor = t_factor_list[valid_time_units.index(time_disp)] - - g_factor_list = [1e-3, 1e3 / self.system.gamma] - g_factor = g_factor_list[valid_grad_units.index(grad_disp)] - - t0 = 0 - label_defined = False - label_idx_to_plot = [] - label_legend_to_plot = [] - label_store = {} - for i in range(len(valid_labels)): - label_store[valid_labels[i]] = 0 - if valid_labels[i] in label.upper(): - label_idx_to_plot.append(i) - label_legend_to_plot.append(valid_labels[i]) - - if len(label_idx_to_plot) != 0: - p = parula.main(len(label_idx_to_plot) + 1) - label_colors_to_plot = p(np.arange(len(label_idx_to_plot))) - cycler = mpl.cycler(color=label_colors_to_plot) - sp11.set_prop_cycle(cycler) - - # Block timings - block_edges = np.cumsum([0] + [x[1] for x in sorted(self.block_durations.items())]) - block_edges_in_range = block_edges[(block_edges >= time_range[0]) * (block_edges <= time_range[1])] - if show_blocks: - for sp in [sp11, sp12, sp13, sp21, sp22, sp23]: - sp.set_xticks(t_factor * block_edges_in_range) - sp.set_xticklabels(sp.get_xticklabels(), rotation=90) - - for block_counter in self.block_events: - block = self.get_block(block_counter) - is_valid = time_range[0] <= t0 + self.block_durations[block_counter] and t0 <= time_range[1] - if is_valid: - if getattr(block, 'label', None) is not None: - for i in range(len(block.label)): - if block.label[i].type == 'labelinc': - label_store[block.label[i].label] += block.label[i].value - else: - label_store[block.label[i].label] = block.label[i].value - label_defined = True - - if getattr(block, 'adc', None) is not None: # ADC - adc = block.adc - # From Pulseq: According to the information from Klaus Scheffler and indirectly from Siemens this - # is the present convention - the samples are shifted by 0.5 dwell - t = adc.delay + (np.arange(int(adc.num_samples)) + 0.5) * adc.dwell - sp11.plot(t_factor * (t0 + t), np.zeros(len(t)), 'rx') - sp13.plot( - t_factor * (t0 + t), - np.angle(np.exp(1j * adc.phase_offset) * np.exp(1j * 2 * np.pi * t * adc.freq_offset)), - 'b.', - markersize=0.25, - ) - - if label_defined and len(label_idx_to_plot) != 0: - arr_label_store = list(label_store.values()) - lbl_vals = np.take(arr_label_store, label_idx_to_plot) - t = t0 + adc.delay + (adc.num_samples - 1) / 2 * adc.dwell - _t = [t_factor * t] * len(lbl_vals) - # Plot each label individually to retrieve each corresponding Line2D object - p = itertools.chain.from_iterable( - [sp11.plot(__t, _lbl_vals, '.') for __t, _lbl_vals in zip(_t, lbl_vals, strict=False)] - ) - if len(label_legend_to_plot) != 0: - sp11.legend(list(p), label_legend_to_plot, loc='upper left') - label_legend_to_plot = [] - - if getattr(block, 'rf', None) is not None: # RF - rf = block.rf - time_center, index_center = calc_rf_center(rf) - time = rf.t - signal = rf.signal - - if signal.shape[0] == 2 and rf.freq_offset != 0: - num_samples = min(int(abs(rf.freq_offset)), 256) - time = np.linspace(time[0], time[-1], num_samples) - signal = np.linspace(signal[0], signal[-1], num_samples) - - if abs(signal[0]) != 0: - signal = np.concatenate(([0], signal)) - time = np.concatenate(([time[0]], time)) - index_center += 1 - - if abs(signal[-1]) != 0: - signal = np.concatenate((signal, [0])) - time = np.concatenate((time, [time[-1]])) - - signal_is_real = max(np.abs(np.imag(signal))) / max(np.abs(np.real(signal))) < 1e-6 - - # Compute time vector with delay applied - time_with_delay = t_factor * (t0 + time + rf.delay) - time_center_with_delay = t_factor * (t0 + time_center + rf.delay) - - # Choose plot behavior based on realness of signal - if signal_is_real: - # Plot real part of signal - sp12.plot(time_with_delay, np.real(signal)) - - # Include sign(real(signal)) factor like MATLAB - phase_corrected = ( - signal - * np.sign(np.real(signal)) - * np.exp(1j * rf.phase_offset) - * np.exp(1j * 2 * math.pi * time * rf.freq_offset) - ) - sc_corrected = ( - signal[index_center] - * np.exp(1j * rf.phase_offset) - * np.exp(1j * 2 * math.pi * time[index_center] * rf.freq_offset) - ) - - sp13.plot( - time_with_delay, - np.angle(phase_corrected), - time_center_with_delay, - np.angle(sc_corrected), - 'xb', - ) - else: - # Plot magnitude of complex signal - sp12.plot(time_with_delay, np.abs(signal)) - - # Plot angle of complex signal - phase_corrected = ( - signal * np.exp(1j * rf.phase_offset) * np.exp(1j * 2 * math.pi * time * rf.freq_offset) - ) - sc_corrected = ( - signal[index_center] - * np.exp(1j * rf.phase_offset) - * np.exp(1j * 2 * math.pi * time[index_center] * rf.freq_offset) - ) - - sp13.plot( - time_with_delay, - np.angle(phase_corrected), - time_center_with_delay, - np.angle(sc_corrected), - 'xb', - ) - - grad_channels = ['gx', 'gy', 'gz'] - for x in range(len(grad_channels)): # Gradients - if getattr(block, grad_channels[x], None) is not None: - grad = getattr(block, grad_channels[x]) - if grad.type == 'grad': - # We extend the shape by adding the first and the last points in an effort of making the - # display a bit less confusing... - time = grad.delay + np.array([0, *grad.tt, grad.shape_dur]) - waveform = g_factor * np.array((grad.first, *grad.waveform, grad.last)) - else: - time = np.array( - cumsum( - 0, - grad.delay, - grad.rise_time, - grad.flat_time, - grad.fall_time, - ) - ) - waveform = g_factor * grad.amplitude * np.array([0, 0, 1, 1, 0]) - [sp21, sp22, sp23][x].plot(t_factor * (t0 + time), waveform) - - # Soft delays - plot as shaded regions with annotations - if getattr(block, 'soft_delay', None) is not None: - soft_delay = block.soft_delay - block_duration = self.block_durations[block_counter] - t_mid = t0 + block_duration / 2 # Middle of the block - - # Add shaded region spanning the soft delay block duration on RF phase subplot - sp13.axvspan(t_factor * t0, t_factor * (t0 + block_duration), alpha=0.2, color='orange') - sp12.axvspan(t_factor * t0, t_factor * (t0 + block_duration), alpha=0.2, color='orange') - sp11.axvspan(t_factor * t0, t_factor * (t0 + block_duration), alpha=0.2, color='orange') - - for sp2x in [sp21, sp22, sp23]: - sp2x.axvspan(t_factor * t0, t_factor * (t0 + block_duration), alpha=0.2, color='orange') - - # Add text annotation with soft delay hint - y_lim = sp13.get_ylim() - y_range = y_lim[1] - y_lim[0] - y_pos = y_lim[0] + 0.1 * y_range - y_text = y_lim[0] + 0.3 * y_range - - sp13.annotate( - f'{soft_delay.hint}', - xy=(t_factor * t_mid, y_pos), - xytext=(t_factor * t_mid, y_text), - ha='center', - va='bottom', - fontsize=8, - bbox={'boxstyle': 'round,pad=0.3', 'facecolor': 'orange', 'alpha': 0.7}, - ) - - t0 += self.block_durations[block_counter] - - # Set axis labels - sp11.set_ylabel('ADC') - sp12.set_ylabel('RF mag (Hz)') - sp13.set_ylabel('RF/ADC phase (rad)') - sp13.set_xlabel(f't ({time_disp})') - sp21.set_ylabel(f'Gx ({grad_disp})') - sp22.set_ylabel(f'Gy ({grad_disp})') - sp23.set_ylabel(f'Gz ({grad_disp})') - sp23.set_xlabel(f't ({time_disp})') - - # Set display limits for all subplots - disp_range = t_factor * np.array([time_range[0], min(t0, time_range[1])]) - for sp in [sp11, sp12, sp13, sp21, sp22, sp23]: - sp.set_xlim(disp_range) - - # Enable grid on all subplots (explicitly set to True, don't toggle) - for sp in [sp11, sp12, sp13, sp21, sp22, sp23]: - sp.grid(True) - - fig1.tight_layout() - fig2.tight_layout() - if save: - fig1.savefig('seq_plot1.jpg') - fig2.savefig('seq_plot2.jpg') - - if plot_now: - plt.show() - - # Always return figures and axes for customization - return fig1, (sp11, sp12, sp13), fig2, (sp21, sp22, sp23) + return SeqPlot(self, label, show_blocks, save, time_range, time_disp, grad_disp, plot_now) def read(self, file_path: str, detect_rf_use: bool = False, remove_duplicates: bool = True) -> None: """ @@ -1326,7 +1139,7 @@ def remove_duplicates(self, in_place: bool = False) -> 'Sequence': for grad_id in seq_copy.grad_library.data: if seq_copy.grad_library.type[grad_id] == 'g': data = seq_copy.grad_library.data[grad_id] - new_data = (data[0], mapping[data[1]], mapping[data[2]], *data[3:]) + new_data = (*data[0:3], mapping[data[3]], mapping[data[4]], data[5]) if data != new_data: seq_copy.grad_library.update(grad_id, None, new_data) @@ -1347,14 +1160,14 @@ def remove_duplicates(self, in_place: bool = False) -> 'Sequence': seq_copy.block_events[block_id][4] = mapping[seq_copy.block_events[block_id][4]] # Filter duplicates in RF library - seq_copy.rf_library, mapping = seq_copy.rf_library.remove_duplicates((6, 0, 0, 0, 6, 6, 6)) + seq_copy.rf_library, mapping = seq_copy.rf_library.remove_duplicates((6, 0, 0, 0, 6, 6, 6, 6, 6, 6)) # Remap RF event IDs for block_id in seq_copy.block_events: seq_copy.block_events[block_id][1] = mapping[seq_copy.block_events[block_id][1]] # Filter duplicates in ADC library - seq_copy.adc_library, mapping = seq_copy.adc_library.remove_duplicates((0, -9, -6, 6, 6, 6)) + seq_copy.adc_library, mapping = seq_copy.adc_library.remove_duplicates((0, -9, -6, 6, 6, 6, 6, 6, 6)) # Remap ADC event IDs for block_id in seq_copy.block_events: @@ -1391,7 +1204,7 @@ def rf_from_lib_data(self, lib_data: list, use: str = str()) -> SimpleNamespace: compressed.num_samples = shape_data[0] compressed.data = shape_data[1:] phase = decompress_shape(compressed) - rf.signal = amplitude * mag * np.exp(1j * 2 * np.pi * phase) + rf.signal = amplitude * mag * np.exp(1j * 2 * math.pi * phase) time_shape = lib_data[3] if time_shape > 0: shape_data = self.shape_library.data[time_shape] @@ -1403,22 +1216,26 @@ def rf_from_lib_data(self, lib_data: list, use: str = str()) -> SimpleNamespace: rf.t = (np.arange(1, len(rf.signal) + 1) - 0.5) * self.rf_raster_time rf.shape_dur = len(rf.signal) * self.rf_raster_time - rf.delay = lib_data[4] - rf.freq_offset = lib_data[5] - rf.phase_offset = lib_data[6] + rf.center = lib_data[4] # new in v150 + rf.delay = lib_data[5] # changed in v150 + rf.freq_ppm = lib_data[6] # new in v150 + rf.phase_ppm = lib_data[7] # new in v150 + rf.freq_offset = lib_data[8] # changed in v150 + rf.phase_offset = lib_data[9] # changed in v150 rf.dead_time = self.system.rf_dead_time rf.ringdown_time = self.system.rf_ringdown_time - if use != '': - use_cases = { - 'e': 'excitation', - 'r': 'refocusing', - 'i': 'inversion', - 's': 'saturation', - 'p': 'preparation', - } - rf.use = use_cases.get(use, 'undefined') + # TODO: fixme : use map built from pp.get_supported_rf_uses() + use_cases = { + 'e': 'excitation', + 'r': 'refocusing', + 'i': 'inversion', + 's': 'saturation', + 'p': 'preparation', + 'o': 'other', + } + rf.use = use_cases.get(use, 'undefined') return rf @@ -1469,16 +1286,23 @@ def rf_times( if block.rf is not None: rf = block.rf - t = rf.delay + calc_rf_center(rf)[0] + + tc = calc_rf_center(rf)[0] + t = rf.delay + tc + + full_freq_offset = rf.freq_offset + rf.freq_ppm * 1e-6 * self.system.gamma * self.system.B0 + full_phase_offset = rf.phase_offset + rf.phase_ppm * 1e-6 * self.system.gamma * self.system.B0 + full_phase_offset = full_phase_offset + 2 * math.pi * full_freq_offset * tc + if not hasattr(rf, 'use') or block.rf.use in [ 'excitation', 'undefined', ]: t_excitation.append(curr_dur + t) - fp_excitation.append([block.rf.freq_offset, block.rf.phase_offset]) + fp_excitation.append([full_freq_offset, full_phase_offset]) elif block.rf.use == 'refocusing': t_refocusing.append(curr_dur + t) - fp_refocusing.append([block.rf.freq_offset, block.rf.phase_offset]) + fp_refocusing.append([full_freq_offset, full_phase_offset]) curr_dur += self.block_durations[block_counter] @@ -1798,11 +1622,13 @@ def waveforms(self, append_RF: bool = False, time_range: Union[List[float], None if block.rf is not None: # RF rf = block.rf + full_freq_offset = rf.freq_offset + rf.freq_ppm * 1e-6 * self.system.gamma * self.system.B0 + full_phase_offset = rf.phase_offset + rf.phase_ppm * 1e-6 * self.system.gamma * self.system.B0 if append_RF: rf_piece = np.array( [ curr_dur + rf.delay + rf.t, - rf.signal * np.exp(1j * (rf.phase_offset + 2 * np.pi * rf.freq_offset * rf.t)), + rf.signal * np.exp(1j * (full_phase_offset + 2 * math.pi * full_freq_offset * rf.t)), ] ) out_len[-1] += len(rf.t) @@ -1939,16 +1765,32 @@ def waveforms_export(self, time_range=(0, np.inf)) -> dict: # is the present convention - the samples are shifted by 0.5 dwell t = adc.delay + (np.arange(int(adc.num_samples)) + 0.5) * adc.dwell adc_t = t0 + t - adc_signal = np.exp(1j * adc.phase_offset) * np.exp(1j * 2 * np.pi * t * adc.freq_offset) + + if adc.phase_modulation is None or len(adc.phase_modulation) == 0: + phase_modulation = 0 + else: + phase_modulation = adc.phase_modulation + + full_freq_offset = np.atleast_1d(adc.freq_offset + adc.freq_ppm * 1e-6 * self.system.B0) + full_phase_offset = np.atleast_1d( + adc.phase_offset + adc.phase_offset * 1e-6 * self.system.B0 + phase_modulation + ) + + adc_signal = np.exp(1j * full_phase_offset) * np.exp(1j * 2 * math.pi * t * full_phase_offset) adc_t_all = np.concatenate((adc_t_all, adc_t)) adc_signal_all = np.concatenate((adc_signal_all, adc_signal)) if block.rf is not None: rf = block.rf + tc, ic = calc_rf_center(rf) t = rf.t + rf.delay tc = tc + rf.delay + full_freq_offset = rf.freq_offset + rf.freq_ppm * 1e-6 * self.system.gamma * self.system.B0 + full_phase_offset = rf.phase_offset + rf.phase_ppm * 1e-6 * self.system.gamma * self.system.B0 + full_phase_offset = full_phase_offset + 2 * math.pi * full_freq_offset * tc + # Debug - visualize # sp12.plot(t_factor * (t0 + t), np.abs(rf.signal)) # sp13.plot(t_factor * (t0 + t), np.angle(rf.signal * np.exp(1j * rf.phase_offset) @@ -1958,7 +1800,7 @@ def waveforms_export(self, time_range=(0, np.inf)) -> dict: # 'xb') rf_t = t0 + t - rf = rf.signal * np.exp(1j * rf.phase_offset) * np.exp(1j * 2 * math.pi * rf.t * rf.freq_offset) + rf = rf.signal * np.exp(1j * (full_phase_offset + 2 * math.pi * full_freq_offset * rf.t)) rf_t_all = np.concatenate((rf_t_all, rf_t)) rf_signal_all = np.concatenate((rf_signal_all, rf)) rf_t_centers = np.concatenate((rf_t_centers, [rf_t[ic]])) diff --git a/src/pypulseq/Sequence/write_seq.py b/src/pypulseq/Sequence/write_seq.py index 6d7e19e..97d59ba 100644 --- a/src/pypulseq/Sequence/write_seq.py +++ b/src/pypulseq/Sequence/write_seq.py @@ -5,7 +5,10 @@ import numpy as np -from pypulseq.supported_labels_rf_use import get_supported_labels +from pypulseq import __version__ +from pypulseq.supported_labels_rf_use import get_supported_labels, get_supported_rf_uses + +version_major, version_minor, version_revision = __version__.split('.')[:3] def write(self, file_name: Union[str, Path], create_signature, remove_duplicates=True) -> Union[str, None]: @@ -100,15 +103,19 @@ def write(self, file_name: Union[str, Path], create_signature, remove_duplicates if len(self.rf_library.data) != 0: output_file.write('# Format of RF events:\n') - output_file.write('# id amplitude mag_id phase_id time_shape_id delay freq phase\n') - output_file.write('# .. Hz .... .... .... us Hz rad\n') + output_file.write('# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use\n') + output_file.write('# .. Hz .. .. .. us us ppm rad/MHz Hz rad ..\n') + output_file.write(f'# Field "use" is the initial of: {" ".join(get_supported_rf_uses()).strip()}\n') output_file.write('[RF]\n') - id_format_str = '{:.0f} {:12g} {:.0f} {:.0f} {:.0f} {:g} {:g} {:g}\n' # Refer lines 20-21 + id_format_str = ( + '{:.0f} {:12g} {:.0f} {:.0f} {:.0f} {:g} {:g} {:g} {:g} {:g} {:g} {:s}\n' # Refer lines 20-21 + ) for k in self.rf_library.data: lib_data1 = self.rf_library.data[k][0:4] - lib_data2 = self.rf_library.data[k][5:7] - delay = round(self.rf_library.data[k][4] / self.rf_raster_time) * self.rf_raster_time * 1e6 - s = id_format_str.format(k, *lib_data1, delay, *lib_data2) + lib_data2 = self.rf_library.data[k][6:10] + center = self.rf_library.data[k][4] * 1e6 # us + delay = round(self.rf_library.data[k][5] / self.rf_raster_time) * self.rf_raster_time * 1e6 + s = id_format_str.format(k, *lib_data1, center, delay, *lib_data2, self.rf_library.type[k]) output_file.write(s) output_file.write('\n') @@ -121,16 +128,16 @@ def write(self, file_name: Union[str, Path], create_signature, remove_duplicates output_file.write( '# time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)\n' ) - output_file.write('# id amplitude amp_shape_id time_shape_id delay\n') - output_file.write('# .. Hz/m .. .. us\n') + output_file.write('# id amplitude first last amp_shape_id time_shape_id delay\n') + output_file.write('# .. Hz/m Hz/m Hz/m .. .. us\n') output_file.write('[GRADIENTS]\n') - id_format_str = '{:.0f} {:12g} {:.0f} {:.0f} {:.0f}\n' # Refer lines 20-21 + id_format_str = '{:.0f} {:12g} {:12g} {:12g} {:.0f} {:.0f} {:.0f}\n' # Refer lines 20-21 keys = np.array(list(self.grad_library.data.keys())) for k in keys[arb_grad_mask]: s = id_format_str.format( k, - *self.grad_library.data[k][:3], - round(self.grad_library.data[k][3] * 1e6), + *self.grad_library.data[k][:5], + round(self.grad_library.data[k][5] * 1e6), ) output_file.write(s) output_file.write('\n') @@ -156,12 +163,12 @@ def write(self, file_name: Union[str, Path], create_signature, remove_duplicates if len(self.adc_library.data) != 0: output_file.write('# Format of ADC events:\n') - output_file.write('# id num dwell delay freq phase\n') - output_file.write('# .. .. ns us Hz rad\n') + output_file.write('# id num dwell delay freqPPM phasePPM freq phase phase_id\n') + output_file.write('# .. .. ns us ppm rad/MHz Hz rad ..\n') output_file.write('[ADC]\n') - id_format_str = '{:.0f} {:.0f} {:.0f} {:.0f} {:g} {:g}\n' # Refer lines 20-21 + id_format_str = '{:.0f} {:.0f} {:.0f} {:.0f} {:g} {:g} {:g} {:g} {:.0f}\n' # Refer lines 20-21 for k in self.adc_library.data: - data = np.multiply(self.adc_library.data[k][0:5], [1, 1e9, 1e6, 1, 1]) + data = np.multiply(self.adc_library.data[k][0:8], [1, 1e9, 1e6, 1, 1, 1, 1, 1]) s = id_format_str.format(k, *data) output_file.write(s) output_file.write('\n') diff --git a/src/pypulseq/add_gradients.py b/src/pypulseq/add_gradients.py index c4e2aac..13bfef0 100644 --- a/src/pypulseq/add_gradients.py +++ b/src/pypulseq/add_gradients.py @@ -86,24 +86,29 @@ def add_gradients( return grad # Find out the general delay of all gradients and other statistics - delays, firsts, lasts, durs, is_trap, is_arb = [], [], [], [], [], [] + delays, firsts, lasts, durs, is_trap, is_arb, is_osa = [], [], [], [], [], [], [] for ii in range(len(grads)): if grads[ii].channel != channel: raise ValueError('Cannot add gradients on different channels.') delays.append(grads[ii].delay) - firsts.append(grads[ii].first) - lasts.append(grads[ii].last) durs.append(calc_duration(grads[ii])) is_trap.append(grads[ii].type == 'trap') if is_trap[-1]: is_arb.append(False) + is_osa.append(False) + firsts.append(0.0) + lasts.append(0.0) else: - tt_rast = grads[ii].tt / system.grad_raster_time - 0.5 - is_arb.append(np.all(np.abs(tt_rast - np.arange(len(tt_rast)))) < eps) + tt_rast = grads[ii].tt / system.grad_raster_time + is_arb.append(np.all(np.abs(tt_rast + 0.5 - np.arange(1, len(tt_rast) + 1))) < eps) + is_osa.append(np.all(np.abs(tt_rast - 0.5 * np.arange(1, len(tt_rast) + 1)) < eps)) + firsts.append(grads[ii].first) + lasts.append(grads[ii].last) # Check if we only have arbitrary grads on irregular time samplings, optionally mixed with trapezoids - if np.all(np.logical_or(is_trap, np.logical_not(is_arb))): + is_etrap = np.logical_and.reduce((np.logical_not(is_trap), np.logical_not(is_arb), np.logical_not(is_osa))) + if np.all(np.logical_or(is_trap, is_etrap)): # Keep shapes still rather simple times = [] for ii in range(len(grads)): @@ -161,21 +166,35 @@ def add_gradients( # Convert to numpy.ndarray for fancy-indexing later on firsts, lasts = np.array(firsts), np.array(lasts) common_delay = np.min(delays) + total_duration = np.max(durs) durs = np.array(durs) # Convert everything to a regularly-sampled waveform waveforms = {} max_length = 0 + + if np.any(is_osa): + target_raster = 0.5 * system.grad_raster_time + else: + target_raster = system.grad_raster_time + for ii in range(len(grads)): g = grads[ii] if g.type == 'grad': - if is_arb[ii]: - waveforms[ii] = g.waveform + if is_arb[ii] or is_osa[ii]: + if np.any(is_osa) and is_arb[ii]: # Porting MATLAB here, maybe a bit ugly + # Interpolate missing samples + idx = np.arange(0, len(g.waveform) - 0.5 + eps, 0.5) + wf = g.waveform + interp_waveform = 0.5 * (wf[np.floor(idx).astype(int)] + wf[np.ceil(idx).astype(int)]) + waveforms[ii] = interp_waveform + else: + waveforms[ii] = g.waveform else: waveforms[ii] = points_to_waveform( amplitudes=g.waveform, times=g.tt, - grad_raster_time=system.grad_raster_time, + grad_raster_time=target_raster, ) elif g.type == 'trap': if g.flat_time > 0: # Triangle or trapezoid @@ -200,7 +219,7 @@ def add_gradients( waveforms[ii] = points_to_waveform( amplitudes=amplitudes, times=times, - grad_raster_time=system.grad_raster_time, + grad_raster_time=target_raster, ) else: raise ValueError('Unknown gradient type') @@ -228,6 +247,9 @@ def add_gradients( max_slew=max_slew, max_grad=max_grad, delay=common_delay, + oversampling=np.any(is_osa), + first=np.sum(firsts[delays == common_delay]), + last=np.sum(lasts[durs == total_duration]), ) # Fix the first and the last values # First is defined by the sum of firsts with the minimal delay (common_delay) diff --git a/src/pypulseq/add_ramps.py b/src/pypulseq/add_ramps.py index c1c3c8d..0bffa6b 100644 --- a/src/pypulseq/add_ramps.py +++ b/src/pypulseq/add_ramps.py @@ -14,6 +14,7 @@ def add_ramps( max_slew: int = 0, rf: Union[SimpleNamespace, None] = None, system: Union[Opts, None] = None, + oversampling: bool = False, ) -> List[np.ndarray]: """ Add segments to the trajectory to ramp to and from the given trajectory. @@ -32,6 +33,8 @@ def add_ramps( Maximum gradient amplitude. max_slew : int, default=0 Maximum slew rate. + oversampling : bool, default=False + Boolean flag to indicate if gradient is oversampled by a factor of 2. Returns ------- @@ -64,8 +67,8 @@ def add_ramps( num_channels = k.shape[0] k = np.vstack((k, np.zeros((3 - num_channels, k.shape[1])))) # Pad with zeros if needed - k_up, ok1 = calc_ramp(k0=np.zeros((3, 2)), k_end=k[:, :2], system=system) - k_down, ok2 = calc_ramp(k0=k[:, -2:], k_end=np.zeros((3, 2)), system=system) + k_up, ok1 = calc_ramp(k0=np.zeros((3, 2)), k_end=k[:, :2], system=system, oversampling=oversampling) + k_down, ok2 = calc_ramp(k0=k[:, -2:], k_end=np.zeros((3, 2)), system=system, oversampling=oversampling) if not (ok1 and ok2): raise RuntimeError('Failed to calculate gradient ramps') diff --git a/src/pypulseq/calc_ramp.py b/src/pypulseq/calc_ramp.py index fd3342e..0b3bbd3 100644 --- a/src/pypulseq/calc_ramp.py +++ b/src/pypulseq/calc_ramp.py @@ -12,6 +12,7 @@ def calc_ramp( max_points: int = 500, max_slew: Union[np.ndarray, None] = None, system: Union[Opts, None] = None, + oversampling: bool = False, ) -> Tuple[np.ndarray, bool]: """ Join the points `k0` and `k_end` in three-dimensional k-space in minimal time, observing the gradient and slew @@ -34,6 +35,8 @@ def calc_ramp( Maximum total slew rate. Either a single value or one value for each coordinate, of shape `[3, 1]`. system : Opts, default=Opts() System limits. + oversampling : bool, default=False + Boolean flag to indicate if gradient is oversampled by a factor of 2. Returns ------- @@ -321,7 +324,10 @@ def __joinright1(k0, k_end, use_points, G0, G_end): if np.all(np.where(max_slew <= 0)): max_slew = [system.max_slew] - grad_raster = system.grad_raster_time + if oversampling: + grad_raster = 0.5 * system.grad_raster_time + else: + grad_raster = system.grad_raster_time if len(max_grad) == 1 and len(max_slew) == 1: mode = 0 diff --git a/src/pypulseq/calc_rf_bandwidth.py b/src/pypulseq/calc_rf_bandwidth.py index 0a0a7bf..cf72148 100644 --- a/src/pypulseq/calc_rf_bandwidth.py +++ b/src/pypulseq/calc_rf_bandwidth.py @@ -1,10 +1,12 @@ import math +import warnings from types import SimpleNamespace from typing import Tuple, Union import numpy as np from pypulseq.calc_rf_center import calc_rf_center +from pypulseq.opts import Opts def calc_rf_bandwidth( @@ -12,6 +14,8 @@ def calc_rf_bandwidth( cutoff: float = 0.5, return_axis: bool = False, return_spectrum: bool = False, + dw: float = 10, + dt: float = 1e-6, ) -> Union[float, Tuple[float, np.ndarray], Tuple[float, np.ndarray, float]]: """ Calculate the spectrum of the RF pulse. Returns the bandwidth of the pulse (calculated by a simple FFT, e.g. @@ -27,6 +31,10 @@ def calc_rf_bandwidth( Boolean flag to indicate if frequency axis of RF pulse will be returned. return_spectrum : bool, default=False Boolean flag to indicate if spectrum of RF pulse will be returned. + dw : float, default=10 + Spectral resolution in (Hz). + dt : float, default=1e-6 + Sampling time in (s). Returns ------- @@ -36,13 +44,22 @@ def calc_rf_bandwidth( """ time_center, _ = calc_rf_center(rf) + if abs(rf.freq_ppm) > np.finfo(float).eps: + warnings.warn( + 'calc_rf_bandwidth(): relying on the system properties, like B0 and gamma, ' + 'stored in the global environment by calling pypulseq.Opts()' + ) + sys = Opts() + full_freq_offset = rf.freq_offset + rf.freq_ppm * 1e-6 * sys.gamma * sys.B0 + else: + full_freq_offset = rf.freq_offset + # Resample the pulse to a reasonable time array - dw = 10 # Hz - dt = 1e-6 # For now, 1 MHz nn = round(1 / dw / dt) tt = np.arange(-math.floor(nn / 2), math.ceil(nn / 2) - 1) * dt - rfs = np.interp(xp=rf.t - time_center, fp=rf.signal, x=tt) + rf_signal = rf.signal * np.exp(1j * rf.phase_offset + 2 * math.pi * full_freq_offset * rf.t) + rfs = np.interp(xp=rf.t - time_center, fp=rf_signal, x=tt) spectrum = np.fft.fftshift(np.fft.fft(np.fft.fftshift(rfs))) w = np.arange(-math.floor(nn / 2), math.ceil(nn / 2) - 1) * dw @@ -53,7 +70,9 @@ def calc_rf_bandwidth( if return_spectrum and not return_axis: return bw, spectrum - if return_axis: + elif return_axis and not return_spectrum: + return bw, w + elif return_spectrum and return_axis: return bw, spectrum, w return bw diff --git a/src/pypulseq/calc_rf_center.py b/src/pypulseq/calc_rf_center.py index bce48bf..cfe2cc1 100644 --- a/src/pypulseq/calc_rf_center.py +++ b/src/pypulseq/calc_rf_center.py @@ -22,6 +22,9 @@ def calc_rf_center(rf: SimpleNamespace) -> Tuple[float, float]: id_center : float Corresponding position of `time_center` in the radio-frequency pulse's envelope. """ + if hasattr(rf, 'center'): + return rf.center, np.argmin(abs(rf.t - rf.center)).item() + # Detect the excitation peak; if i is a plateau take its center rf_max = np.max(np.abs(rf.signal)) i_peak = np.where(np.abs(rf.signal) >= rf_max * 0.99999)[0] diff --git a/src/pypulseq/make_adc.py b/src/pypulseq/make_adc.py index a7bb670..231386b 100644 --- a/src/pypulseq/make_adc.py +++ b/src/pypulseq/make_adc.py @@ -4,6 +4,8 @@ from typing import List, Optional, Tuple, Union from warnings import warn +import numpy as np + from pypulseq.opts import Opts from pypulseq.utils.tracing import trace, trace_enabled @@ -16,6 +18,9 @@ def make_adc( freq_offset: float = 0, phase_offset: float = 0, system: Union[Opts, None] = None, + freq_ppm: float = 0, + phase_ppm: float = 0, + phase_modulation: np.ndarray = None, ) -> SimpleNamespace: """ Create an ADC readout event. @@ -36,6 +41,13 @@ def make_adc( Frequency offset of ADC readout event. phase_offset : float, default=0 Phase offset of ADC readout event. + freq_ppm : float, default=0 + PPM frequency offset of ADC readout event. + phase_ppm : float, default=0 + PPM phase offset of ADC readout event. + phase_modulation : numpy.ndarray, default=None + Phase modulation array for FOV shifting. + If provided, it must have `num_samples` number of samples. Returns ------- @@ -57,11 +69,17 @@ def make_adc( adc.delay = delay adc.freq_offset = freq_offset adc.phase_offset = phase_offset + adc.freq_ppm = freq_ppm + adc.phase_ppm = phase_ppm adc.dead_time = system.adc_dead_time if (dwell == 0 and duration == 0) or (dwell > 0 and duration > 0): raise ValueError('Either dwell or duration must be defined') + if phase_modulation is not None and len(phase_modulation) != num_samples: + raise ValueError('ADC Phase modulation vector must have the same length as the number of samples') + adc.phase_modulation = phase_modulation + if duration > 0: adc.dwell = duration / num_samples diff --git a/src/pypulseq/make_adiabatic_pulse.py b/src/pypulseq/make_adiabatic_pulse.py index 0eaea43..a64bd5c 100644 --- a/src/pypulseq/make_adiabatic_pulse.py +++ b/src/pypulseq/make_adiabatic_pulse.py @@ -30,7 +30,9 @@ def make_adiabatic_pulse( return_gz: bool = False, slice_thickness: float = 0, system: Union[Opts, None] = None, - use: str = str(), + use: str = 'inversion', + freq_ppm: float = 0, + phase_ppm: float = 0, ) -> Union[ SimpleNamespace, Tuple[SimpleNamespace, SimpleNamespace, SimpleNamespace], @@ -115,6 +117,10 @@ def make_adiabatic_pulse( System limits. use : str Whether it is a 'refocusing' pulse (for k-space calculation). + freq_ppm : float, default=0 + PPM frequency offset. + phase_ppm : float, default=0 + PPM phase offset. Returns ------- @@ -209,10 +215,13 @@ def make_adiabatic_pulse( rf.shape_dur = n_samples * dwell rf.freq_offset = freq_offset rf.phase_offset = phase_offset + rf.freq_ppm = freq_ppm + rf.phase_ppm = phase_ppm rf.dead_time = system.rf_dead_time rf.ringdown_time = system.rf_ringdown_time rf.delay = delay - rf.use = use if use != '' else 'inversion' + rf.center, _ = calc_rf_center(rf) + rf.use = use if rf.dead_time > rf.delay: warn( f'Specified RF delay {rf.delay * 1e6:.2f} us is less than the dead time {rf.dead_time * 1e6:.0f} us. Delay was increased to the dead time.', diff --git a/src/pypulseq/make_arbitrary_grad.py b/src/pypulseq/make_arbitrary_grad.py index 720e94c..75ba236 100644 --- a/src/pypulseq/make_arbitrary_grad.py +++ b/src/pypulseq/make_arbitrary_grad.py @@ -1,3 +1,4 @@ +# import warnings from types import SimpleNamespace from typing import Union @@ -17,6 +18,7 @@ def make_arbitrary_grad( max_grad: Union[float, None] = None, max_slew: Union[float, None] = None, system: Union[Opts, None] = None, + oversampling: bool = False, ) -> SimpleNamespace: """ Creates a gradient event from an arbitrary waveform. @@ -51,6 +53,8 @@ def make_arbitrary_grad( Will default to `system.max_slew` if not provided. delay : float, default=0 Delay in seconds (s). + oversampling : bool, default=False + Boolean flag to indicate if gradient is oversampled by a factor of 2. Returns ------- @@ -76,28 +80,59 @@ def make_arbitrary_grad( if channel not in ['x', 'y', 'z']: raise ValueError(f'Invalid channel. Must be one of x, y or z. Passed: {channel}') - slew_rate = np.diff(waveform) / system.grad_raster_time + if first is None: + # warnings.warn( + # 'it will be compulsory to provide the first point of the gradient shape in the future releases; finding the first by extrapolation for now...', + # FutureWarning, + # ) + if oversampling: + first = 2 * waveform[0] - waveform[1] # extrapolate by 1 gradient raster + else: + first = (3 * waveform[0] - waveform[1]) * 0.5 # extrapolate by 1/2 gradient of the raster + + if last is None: + # warnings.warn( + # 'it will be compulsory to provide the last point of the gradient shape in the future releases; finding the last by extrapolation for now...', + # FutureWarning, + # ) + if oversampling: + last = 2 * waveform[-1] - waveform[-2] # extrapolate by 1 gradient raster + else: + last = (3 * waveform[-1] - waveform[-2]) * 0.5 # extrapolate by 1/2 gradient of the raster + + # Slew rate calculation + if oversampling: + slew_rate = np.concatenate([[first - waveform[0]], np.diff(waveform), [last - waveform[-1]]]) / ( + system.grad_raster_time * 2 + ) + else: + slew_rate = ( + np.concatenate([[2 * (first - waveform[0])], np.diff(waveform), [2 * (waveform[-1] - last)]]) + / system.grad_raster_time + ) + if max(abs(slew_rate)) > max_slew * (1 + eps): raise ValueError(f'Slew rate violation {max(abs(slew_rate)) / max_slew * 100}') if max(abs(waveform)) > max_grad + eps: raise ValueError(f'Gradient amplitude violation {max(abs(waveform)) / max_grad * 100}') - if first is None: - first = (3 * waveform[0] - waveform[1]) * 0.5 # linear extrapolation - - if last is None: - last = (3 * waveform[-1] - waveform[-2]) * 0.5 # linear extrapolation - grad = SimpleNamespace() grad.type = 'grad' grad.channel = channel grad.waveform = waveform grad.delay = delay - grad.tt = (np.arange(len(waveform)) + 0.5) * system.grad_raster_time - grad.shape_dur = len(waveform) * system.grad_raster_time + if oversampling: + if len(waveform) % 2 == 0: + raise ValueError('When oversampling is active, waveform must have an odd number of samples') + grad.area = (waveform[::2] * system.grad_raster_time).sum() + grad.tt = np.arange(1, len(waveform) + 1) * 0.5 * system.grad_raster_time + grad.shape_dur = (len(waveform) + 1) * 0.5 * system.grad_raster_time + else: + grad.area = (waveform * system.grad_raster_time).sum() + grad.tt = (np.arange(len(waveform)) + 0.5) * system.grad_raster_time + grad.shape_dur = len(waveform) * system.grad_raster_time grad.first = first grad.last = last - grad.area = (waveform * system.grad_raster_time).sum() if trace_enabled(): grad.trace = trace() diff --git a/src/pypulseq/make_arbitrary_rf.py b/src/pypulseq/make_arbitrary_rf.py index 8c42ce2..1c424f3 100644 --- a/src/pypulseq/make_arbitrary_rf.py +++ b/src/pypulseq/make_arbitrary_rf.py @@ -6,6 +6,7 @@ import numpy as np +from pypulseq.calc_rf_center import calc_rf_center from pypulseq.make_trapezoid import make_trapezoid from pypulseq.opts import Opts from pypulseq.supported_labels_rf_use import get_supported_rf_uses @@ -27,7 +28,10 @@ def make_arbitrary_rf( slice_thickness: float = 0, system: Union[Opts, None] = None, time_bw_product: float = 0, - use: str = str(), + use: str = 'undefined', + freq_ppm: float = 0, + phase_ppm: float = 0, + center: Union[float, None] = None, ) -> Union[SimpleNamespace, Tuple[SimpleNamespace, SimpleNamespace]]: """ Create an RF pulse with the given pulse shape. @@ -63,8 +67,14 @@ def make_arbitrary_rf( System limits. time_bw_product : float, default=4 Time-bandwidth product. - use : str, default=str() + use : str, default='undefined' Use of arbitrary radio-frequency pulse event. Must be one of 'excitation', 'refocusing' or 'inversion'. + freq_ppm : float, default=0 + PPM frequency offset. + phase_ppm : float, default=0 + PPM phase offset. + center : float, default=None + RF center (s). Returns ------- @@ -110,12 +120,12 @@ def make_arbitrary_rf( rf.shape_dur = duration rf.freq_offset = freq_offset rf.phase_offset = phase_offset + rf.freq_ppm = freq_ppm + rf.phase_ppm = phase_ppm rf.dead_time = system.rf_dead_time rf.ringdown_time = system.rf_ringdown_time rf.delay = delay - - if use != '': - rf.use = use + rf.use = use if rf.dead_time > rf.delay: warn( @@ -124,6 +134,15 @@ def make_arbitrary_rf( ) rf.delay = rf.dead_time + if center is not None: + rf.center = center + if rf.center < 0: + rf.center = 0 + if rf.center > rf.shape_dur: + rf.center = rf.shape_dur + else: + rf.center, _ = calc_rf_center(rf) + if return_gz: if slice_thickness <= 0: raise ValueError('Slice thickness must be provided.') diff --git a/src/pypulseq/make_block_pulse.py b/src/pypulseq/make_block_pulse.py index 8fad36c..0359547 100644 --- a/src/pypulseq/make_block_pulse.py +++ b/src/pypulseq/make_block_pulse.py @@ -18,7 +18,9 @@ def make_block_pulse( freq_offset: float = 0, phase_offset: float = 0, system: Union[Opts, None] = None, - use: str = str(), + use: str = 'undefined', + freq_ppm: float = 0, + phase_ppm: float = 0, ) -> SimpleNamespace: """ Create a block (RECT or hard) pulse. @@ -46,8 +48,12 @@ def make_block_pulse( Phase offset Hertz (Hz). system : Opts, default=Opts() System limits. - use : str, default=str() + use : str, default='undefined' Use of radio-frequency block pulse event. + freq_ppm : float, default=0 + PPM frequency offset. + phase_ppm : float, default=0 + PPM phase offset. Returns ------- @@ -101,12 +107,13 @@ def make_block_pulse( rf.shape_dur = t[-1] rf.freq_offset = freq_offset rf.phase_offset = phase_offset + rf.freq_ppm = freq_ppm + rf.phase_ppm = phase_ppm rf.dead_time = system.rf_dead_time rf.ringdown_time = system.rf_ringdown_time rf.delay = delay - - if use != '': - rf.use = use + rf.center = rf.shape_dur / 2 + rf.use = use if rf.dead_time > rf.delay: warn( diff --git a/src/pypulseq/make_gauss_pulse.py b/src/pypulseq/make_gauss_pulse.py index d549525..f7c80b3 100644 --- a/src/pypulseq/make_gauss_pulse.py +++ b/src/pypulseq/make_gauss_pulse.py @@ -28,7 +28,9 @@ def make_gauss_pulse( slice_thickness: float = 0, system: Union[Opts, None] = None, time_bw_product: float = 4, - use: str = str(), + use: str = 'undefined', + freq_ppm: float = 0, + phase_ppm: float = 0, ) -> Union[ SimpleNamespace, Tuple[SimpleNamespace, SimpleNamespace, SimpleNamespace], @@ -70,8 +72,12 @@ def make_gauss_pulse( System limits. time_bw_product : int, default=4 Time-bandwidth product. - use : str, default=str() + use : str, default='undefined' Use of radio-frequency gauss pulse event. Must be one defined in pypulseq.supported_labels_rf_use.get_supported_rf_uses. + freq_ppm : float, default=0 + PPM frequency offset. + phase_ppm : float, default=0 + PPM phase offset. Returns ------- @@ -117,9 +123,12 @@ def make_gauss_pulse( rf.shape_dur = n_samples * dwell rf.freq_offset = freq_offset rf.phase_offset = phase_offset + rf.freq_ppm = freq_ppm + rf.phase_ppm = phase_ppm rf.dead_time = system.rf_dead_time rf.ringdown_time = system.rf_ringdown_time rf.delay = delay + rf.center = duration * center_pos if use != '': rf.use = use diff --git a/src/pypulseq/make_sigpy_pulse.py b/src/pypulseq/make_sigpy_pulse.py index e5fd033..715a444 100644 --- a/src/pypulseq/make_sigpy_pulse.py +++ b/src/pypulseq/make_sigpy_pulse.py @@ -17,6 +17,7 @@ from pypulseq.make_trapezoid import make_trapezoid from pypulseq.opts import Opts from pypulseq.sigpy_pulse_opts import SigpyPulseOpts +from pypulseq.supported_labels_rf_use import get_supported_rf_uses from pypulseq.utils.tracing import trace, trace_enabled @@ -34,8 +35,10 @@ def sigpy_n_seq( system: Union[Opts, None] = None, time_bw_product: float = 4, pulse_cfg: Union[SigpyPulseOpts, None] = None, - use: str = str(), + use: str = 'undefined', plot: bool = True, + freq_ppm: float = 0, + phase_ppm: float = 0, ) -> Union[SimpleNamespace, Tuple[SimpleNamespace, SimpleNamespace, SimpleNamespace]]: """ Creates a radio-frequency sinc pulse event using the sigpy rf pulse library and optionally accompanying slice select, slice select rephasing @@ -88,10 +91,14 @@ def sigpy_n_seq( Band separation. SMS only. - phs_0_pt: str, optional, default='None' Phase 0 point. SMS only. - use : str, optional, default=str() + use : str, default='undefined' Use of radio-frequency sinc pulse. Must be one of 'excitation', 'refocusing' or 'inversion'. plot: bool, optional, default=True Show sigpy plot outputs + freq_ppm : float, default=0 + PPM frequency offset. + phase_ppm : float, default=0 + PPM phase offset. Returns ------- @@ -114,11 +121,9 @@ def sigpy_n_seq( if pulse_cfg is None: pulse_cfg = SigpyPulseOpts() - valid_use_pulses = ['excitation', 'refocusing', 'inversion'] - if use != '' and use not in valid_use_pulses: - raise ValueError( - f"Invalid use parameter. Must be one of 'excitation', 'refocusing' or 'inversion'. Passed: {use}" - ) + valid_pulse_uses = get_supported_rf_uses() + if use != '' and use not in valid_pulse_uses: + raise ValueError(f'Invalid use parameter. Must be one of {valid_pulse_uses}. Passed: {use}') if pulse_cfg.pulse_type == 'slr': [signal, t, _] = make_slr( @@ -146,12 +151,13 @@ def sigpy_n_seq( rfp.shape_dur = t[-1] rfp.freq_offset = freq_offset rfp.phase_offset = phase_offset + rfp.freq_ppm = freq_ppm + rfp.phase_ppm = phase_ppm rfp.dead_time = system.rf_dead_time rfp.ringdown_time = system.rf_ringdown_time rfp.delay = delay - - if use != '': - rfp.use = use + rfp.center = center_pos * rfp.shape_dur + rfp.use = use if rfp.dead_time > rfp.delay: warn( diff --git a/src/pypulseq/make_sinc_pulse.py b/src/pypulseq/make_sinc_pulse.py index c791b00..fdb4c1c 100644 --- a/src/pypulseq/make_sinc_pulse.py +++ b/src/pypulseq/make_sinc_pulse.py @@ -27,7 +27,9 @@ def make_sinc_pulse( slice_thickness: float = 0, system: Union[Opts, None] = None, time_bw_product: float = 4, - use: str = str(), + use: str = 'undefined', + freq_ppm: float = 0, + phase_ppm: float = 0, ) -> Union[ SimpleNamespace, Tuple[SimpleNamespace, SimpleNamespace, SimpleNamespace], @@ -66,8 +68,12 @@ def make_sinc_pulse( System limits. Default is a system limits object initialized to default values. time_bw_product : float, default=4 Time-bandwidth product. - use : str, default=str() + use : str, default='undefined' Use of radio-frequency sinc pulse. Must be one of 'excitation', 'refocusing' or 'inversion'. + freq_ppm : float, default=0 + PPM frequency offset. + phase_ppm : float, default=0 + PPM phase offset. See also `pypulseq.Sequence.sequence.Sequence.add_block()`. @@ -116,12 +122,13 @@ def make_sinc_pulse( rf.shape_dur = n_samples * dwell rf.freq_offset = freq_offset rf.phase_offset = phase_offset + rf.freq_ppm = freq_ppm + rf.phase_ppm = phase_ppm rf.dead_time = system.rf_dead_time rf.ringdown_time = system.rf_ringdown_time rf.delay = delay - - if use != str(): - rf.use = use + rf.center = duration * center_pos + rf.use = use if rf.dead_time > rf.delay: warn( diff --git a/src/pypulseq/scale_grad.py b/src/pypulseq/scale_grad.py index 37af6e2..3245879 100644 --- a/src/pypulseq/scale_grad.py +++ b/src/pypulseq/scale_grad.py @@ -1,8 +1,14 @@ from copy import copy from types import SimpleNamespace +from typing import Union +import numpy as np -def scale_grad(grad: SimpleNamespace, scale: float) -> SimpleNamespace: +from pypulseq import eps +from pypulseq.opts import Opts + + +def scale_grad(grad: SimpleNamespace, scale: float, system: Union[Opts, None] = None) -> SimpleNamespace: """ Scales the gradient with the scalar. @@ -12,6 +18,8 @@ def scale_grad(grad: SimpleNamespace, scale: float) -> SimpleNamespace: Gradient event to be scaled. scale : float Scaling factor. + system : Opts, default=Opts() + System limits. Returns ------- @@ -22,12 +30,36 @@ def scale_grad(grad: SimpleNamespace, scale: float) -> SimpleNamespace: scaled_grad = copy(grad) if scaled_grad.type == 'trap': scaled_grad.amplitude = scaled_grad.amplitude * scale - scaled_grad.area = scaled_grad.area * scale scaled_grad.flat_area = scaled_grad.flat_area * scale + if system is not None: + if abs(scaled_grad.amplitude) > system.max_grad: + raise ValueError( + f'scale_grad: maximum amplitude exceeded {100 * abs(scaled_grad.amplitude) / system.max_grad} %' + ) + if ( + abs(grad.amplitude) > eps + and abs(scaled_grad.amplitude) / min(scaled_grad.rise_time, scaled_grad.fall_time) > system.max_slew + ): + raise ValueError( + 'mr.scale_grad: maximum slew rate exceeded {100 * abs(scaled_grad.amplitude) / min(scaled_grad.rise_time, scaled_grad.fall_time) / system.max_slew} %' + ) + else: scaled_grad.waveform = scaled_grad.waveform * scale scaled_grad.first = scaled_grad.first * scale scaled_grad.last = scaled_grad.last * scale + if system is not None: + if max(abs(scaled_grad.waveform)) > system.max_grad: + raise ValueError( + f'scale_grad: maximum amplitude exceeded {100 * max(abs(scaled_grad.waveform)) / system.max_grad} %' + ) + if max(abs(scaled_grad.waveform)) > eps: + scaled_grad_max_abs_slew = max(abs(np.diff(scaled_grad.waveform) / np.diff(grad.tt))) + if scaled_grad_max_abs_slew > system.max_slew: + raise ValueError( + f'scale_grad: maximum slew rate exceeded {100 * scaled_grad_max_abs_slew / system.max_slew} %' + ) + scaled_grad.area = scaled_grad.area * scale if hasattr(scaled_grad, 'id'): delattr(scaled_grad, 'id') diff --git a/src/pypulseq/supported_labels_rf_use.py b/src/pypulseq/supported_labels_rf_use.py index eb62668..13aa2cb 100644 --- a/src/pypulseq/supported_labels_rf_use.py +++ b/src/pypulseq/supported_labels_rf_use.py @@ -34,11 +34,11 @@ def get_supported_labels() -> Tuple[str, str, str, str, str, str, str, str, str, ) -def get_supported_rf_uses() -> Tuple[str, str, str, str, str]: +def get_supported_rf_uses() -> Tuple[str, str, str, str, str, str, str]: """ Returns ------- tuple Supported RF use labels. """ - return 'excitation', 'refocusing', 'inversion', 'saturation', 'preparation' + return 'excitation', 'refocusing', 'inversion', 'saturation', 'preparation', 'other', 'undefined' diff --git a/src/pypulseq/utils/paper_plot.py b/src/pypulseq/utils/paper_plot.py new file mode 100644 index 0000000..a83aa8c --- /dev/null +++ b/src/pypulseq/utils/paper_plot.py @@ -0,0 +1,121 @@ +from typing import Tuple + +import matplotlib.pyplot as plt +import numpy as np + + +def paper_plot( + seq, + time_range: Tuple[float] = (0, np.inf), + line_width: float = 1.2, + axes_color: Tuple[float] = (0.5, 0.5, 0.5), + rf_color: str = 'black', + gx_color: str = 'blue', + gy_color: str = 'red', + gz_color: Tuple[float] = (0, 0.5, 0.3), + rf_plot: str = 'abs', +): + """ + Plot sequence using paper-style formatting (minimalist, high-contrast layout). + + Parameters + ---------- + seq : Sequence + The Pulseq sequence object to plot. + time_range : iterable, default=(0, np.inf) + Time range (x-axis limits) for plotting the sequence. + Default is 0 to infinity (entire sequence). + line_width : float, default=1.2 + Line width used in plots. + axes_color : color, default=(0.5, 0.5, 0.5) + Color of horizontal zero axes (e.g., gray). + rf_color : color, default='black' + Color for RF and ADC events. + gx_color : color, default='blue' + Color for gradient X waveform. + gy_color : color, default='red' + Color for gradient Y waveform. + gz_color : color, default=(0, 0.5, 0.3) + Color for gradient Z waveform. + rf_plot : {'abs', 'real', 'imag'}, default='abs' + Determines how to plot RF waveforms (magnitude, real or imaginary part). + + """ + # Get waveform data + wave_data, _, _, t_adc, _ = seq.waveforms_and_times(append_RF=True, time_range=time_range) + + # Max amplitudes for scaling + gwm = np.max(np.abs(np.concatenate(wave_data[:3], axis=1)), axis=1) + gwm[0] = max(gwm[0], t_adc[-1]) + rfm = np.max(np.abs(wave_data[3]), axis=1) + + # Handle complex RF + if rf_plot == 'real': + rf_waveform = np.real(wave_data[3][1]) + elif rf_plot == 'imag': + rf_waveform = np.imag(wave_data[3][1]) + else: + rf_waveform = np.abs(wave_data[3][1]) + wave_data[3] = np.stack((wave_data[3][0].real, rf_waveform), axis=0) + + # Clean waveforms by inserting NaNs between zero plateaus + for i in range(4): + data = wave_data[i] + j = data.shape[1] - 1 + while j > 0: + if data[1, j] == 0 and data[1, j - 1] == 0: + midpoint = 0.5 * (data[0, j] + data[0, j - 1]) + data = np.hstack([data[:, :j], np.array([[midpoint], [np.nan]]), data[:, j:]]) + wave_data[i] = data + j -= 1 + + # Create figure + fig = plt.figure(figsize=(12, 10), constrained_layout=True) + fig.patch.set_facecolor('white') + spec = fig.add_gridspec(nrows=4, ncols=1, hspace=0.0) + axes = [] + + def format_axis(ax, xlim, ylim): + ax.set_xlim(xlim) + ax.set_ylim(ylim) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_facecolor('white') + ax.spines[:].set_visible(False) + + # ADC + ax = fig.add_subplot(spec[0]) + ax.vlines(t_adc, ymin=0, ymax=rfm[1] / 5, color=rf_color, lw=line_width / 4, zorder=2.5) + + # RF + ax.plot([-0.01 * gwm[0], 1.01 * gwm[0]], [0, 0], color=axes_color, lw=line_width / 5) + ax.plot(wave_data[3][0], wave_data[3][1], color=rf_color, lw=line_width) + + # Format RF + ADC + format_axis(ax, [-0.03 * gwm[0], 1.03 * gwm[0]], [-1.03 * rfm[1], 1.03 * rfm[1]]) + axes.append(ax) + + # Gradient Z + ax = fig.add_subplot(spec[1]) + ax.plot([-0.01 * gwm[0], 1.01 * gwm[0]], [0, 0], color=axes_color, lw=line_width / 5) + ax.plot(wave_data[2][0], wave_data[2][1], color=gz_color, lw=line_width) + format_axis(ax, [-0.03 * gwm[0], 1.03 * gwm[0]], [-1.03 * gwm[1], 1.03 * gwm[1]]) + axes.append(ax) + + # Gradient Y + ax = fig.add_subplot(spec[2]) + ax.plot([-0.01 * gwm[0], 1.01 * gwm[0]], [0, 0], color=axes_color, lw=line_width / 5) + ax.plot(wave_data[1][0], wave_data[1][1], color=gy_color, lw=line_width) + format_axis(ax, [-0.03 * gwm[0], 1.03 * gwm[0]], [-1.03 * gwm[1], 1.03 * gwm[1]]) + axes.append(ax) + + # Gradient X + ax = fig.add_subplot(spec[3]) + ax.plot([-0.01 * gwm[0], 1.01 * gwm[0]], [0, 0], color=axes_color, lw=line_width / 5) + ax.plot(wave_data[0][0], wave_data[0][1], color=gx_color, lw=line_width) + format_axis(ax, [-0.03 * gwm[0], 1.03 * gwm[0]], [-1.03 * gwm[1], 1.03 * gwm[1]]) + axes.append(ax) + + # Link X-axes (time axis) + for ax in axes[1:]: + ax.sharex(axes[0]) diff --git a/src/pypulseq/utils/seq_plot.py b/src/pypulseq/utils/seq_plot.py new file mode 100644 index 0000000..7d2e7e6 --- /dev/null +++ b/src/pypulseq/utils/seq_plot.py @@ -0,0 +1,399 @@ +import itertools +import math + +import matplotlib as mpl +import numpy as np +from matplotlib import pyplot as plt + +from pypulseq.calc_rf_center import calc_rf_center +from pypulseq.Sequence import parula +from pypulseq.supported_labels_rf_use import get_supported_labels +from pypulseq.utils.cumsum import cumsum + +try: + import mplcursors + + __MPLCURSORS_AVAILABLE__ = True +except ImportError: + __MPLCURSORS_AVAILABLE__ = False + + +class SeqPlot: + """ + Interactive plotter for a Pulseq `Sequence` object. + + Parameters + ---------- + seq : Sequence + The Pulseq sequence object to plot. + label : str, default=str() + Plot label values for ADC events: in this example for LIN and REP labels; other valid labes are accepted as + a comma-separated list. + save : bool, default=False + Boolean flag indicating if plots should be saved. The two figures will be saved as JPG with numerical + suffixes to the filename 'seq_plot'. + show_blocks : bool, default=False + Boolean flag to indicate if grid and tick labels at the block boundaries are to be plotted. + time_range : iterable, default=(0, np.inf) + Time range (x-axis limits) for plotting the sequence. Default is 0 to infinity (entire sequence). + time_disp : str, default='s' + Time display type, must be one of `s`, `ms` or `us`. + grad_disp : str, default='s' + Gradient display unit, must be one of `kHz/m` or `mT/m`. + plot_now : bool, default=True + If true, function immediately shows the plots, blocking the rest of the code until plots are exited. + If false, plots are shown when plt.show() is called. Useful if plots are to be modified. + plot_type : str, default='Gradient' + Gradients display type, must be one of either 'Gradient' or 'Kspace'. + + Attributes + ---------- + fig1 : matplotlib.figure.Figure + Figure containing RF and ADC channels. + fig2 : matplotlib.figure.Figure + Figure containing Gradient or K-space channels. + ax1 : matplotlib.axes.Axes + Axes for fig1. + ax2 : matplotlib.axes.Axes + Axes for fig2. + """ + + def __init__( + self, + seq, + label: str = str(), + show_blocks: bool = False, + save: bool = False, + time_range=(0, np.inf), + time_disp: str = 's', + grad_disp: str = 'kHz/m', + plot_now: bool = True, + ): + self.seq = seq + self._cursors = [] + + self.fig1, self.fig2 = _seq_plot( + seq, + label=label, + save=save, + show_blocks=show_blocks, + time_range=time_range, + time_disp=time_disp, + grad_disp=grad_disp, + ) + self.ax1 = self.fig1.axes[0] + self.ax2 = self.fig2.axes[0] + + if __MPLCURSORS_AVAILABLE__: + self._setup_cursor(self.fig1) + self._setup_cursor(self.fig2) + + if plot_now: + self.show() + + def show(self): + plt.show() + + def _setup_cursor(self, fig): + for ax in fig.axes: + lines = ax.get_lines() + for line in lines: + cursor = mplcursors.cursor(line, multiple=True) + cursor.connect('add', lambda sel: self._on_datatip(sel)) + self._cursors.append(cursor) + + def _on_datatip(self, sel): + artist = sel.artist + ax = artist.axes + x, y = sel.target + ylabel = ax.get_ylabel().lower() + + if ylabel.startswith('adc') or ( + ylabel.startswith('rf/adc') and artist.get_linestyle() == 'none' and artist.get_marker() == '.' + ): + field = 'adc' + else: + field = ylabel[:2] + + t0 = artist.get_xdata()[0] if hasattr(artist, 'get_xdata') else x + block_index = self.seq.find_block_by_time(t0) + rb = self.seq.get_raw_block_content_IDs(block_index) + + lines_txt = [f't: {x:.3f}', f'Y: {y:.3f}'] + + val = getattr(rb, field, None) + if val is not None: + try: + if field[0] == 'a': + name = self.seq.adc_id2name_map[val] + elif field[0] == 'r': + name = self.seq.rf_id2name_map[val] + else: + name = self.seq.grad_id2name_map[val] + + lines_txt.append(f"blk: {block_index} {field}_id: {val} '{name}'") + except Exception: + lines_txt.append(f'blk: {block_index} {field}_id: {val}') + else: + lines_txt.append(f'blk: {block_index}') + + sel.annotation.set_text('\n'.join(lines_txt)) + self._update_guides() + + def _update_guides(self): + for ax in (self.ax1, self.ax2): + ax.relim() + ax.autoscale_view() + + for fig in (self.fig1, self.fig2): + fig.canvas.draw_idle() + + +def _seq_plot( + seq, + label, + show_blocks, + save, + time_range, + time_disp, + grad_disp, +): + mpl.rcParams['lines.linewidth'] = 0.75 # Set default Matplotlib linewidth + + valid_time_units = ['s', 'ms', 'us'] + valid_grad_units = ['kHz/m', 'mT/m'] + valid_labels = get_supported_labels() + if not all(isinstance(x, (int, float)) for x in time_range) or len(time_range) != 2: + raise ValueError('Invalid time range') + if time_disp not in valid_time_units: + raise ValueError('Unsupported time unit') + + if grad_disp not in valid_grad_units: + raise ValueError('Unsupported gradient unit. Supported gradient units are: ' + str(valid_grad_units)) + + fig1, fig2 = plt.figure(), plt.figure() + sp11 = fig1.add_subplot(311) + sp12 = fig1.add_subplot(312, sharex=sp11) + sp13 = fig1.add_subplot(313, sharex=sp11) + fig2_subplots = [ + fig2.add_subplot(311, sharex=sp11), + fig2.add_subplot(312, sharex=sp11), + fig2.add_subplot(313, sharex=sp11), + ] + + t_factor_list = [1, 1e3, 1e6] + t_factor = t_factor_list[valid_time_units.index(time_disp)] + + g_factor_list = [1e-3, 1e3 / seq.system.gamma] + g_factor = g_factor_list[valid_grad_units.index(grad_disp)] + + t0 = 0 + label_defined = False + label_idx_to_plot = [] + label_legend_to_plot = [] + label_store = {} + for i in range(len(valid_labels)): + label_store[valid_labels[i]] = 0 + if valid_labels[i] in label.upper(): + label_idx_to_plot.append(i) + label_legend_to_plot.append(valid_labels[i]) + + if len(label_idx_to_plot) != 0: + p = parula.main(len(label_idx_to_plot) + 1) + label_colors_to_plot = p(np.arange(len(label_idx_to_plot))) + cycler = mpl.cycler(color=label_colors_to_plot) + sp11.set_prop_cycle(cycler) + + # Block timings + block_edges = np.cumsum([0] + [x[1] for x in sorted(seq.block_durations.items())]) + block_edges_in_range = block_edges[(block_edges >= time_range[0]) * (block_edges <= time_range[1])] + if show_blocks: + for sp in [sp11, sp12, sp13, *fig2_subplots]: + sp.set_xticks(t_factor * block_edges_in_range) + sp.set_xticklabels(sp.get_xticklabels(), rotation=90) + + for block_counter in seq.block_events: + block = seq.get_block(block_counter) + is_valid = time_range[0] <= t0 + seq.block_durations[block_counter] and t0 <= time_range[1] + if is_valid: + if getattr(block, 'label', None) is not None: + for i in range(len(block.label)): + if block.label[i].type == 'labelinc': + label_store[block.label[i].label] += block.label[i].value + else: + label_store[block.label[i].label] = block.label[i].value + label_defined = True + + if getattr(block, 'adc', None) is not None: # ADC + adc = block.adc + # From Pulseq: According to the information from Klaus Scheffler and indirectly from Siemens this + # is the present convention - the samples are shifted by 0.5 dwell + t = adc.delay + (np.arange(int(adc.num_samples)) + 0.5) * adc.dwell + sp11.plot(t_factor * (t0 + t), np.zeros(len(t)), 'rx') + + if adc.phase_modulation is None or len(adc.phase_modulation) == 0: + phase_modulation = 0 + else: + phase_modulation = adc.phase_modulation + + full_freq_offset = np.atleast_1d(adc.freq_offset + adc.freq_ppm * 1e-6 * seq.system.B0) + full_phase_offset = np.atleast_1d( + adc.phase_offset + adc.phase_offset * 1e-6 * seq.system.B0 + phase_modulation + ) + + sp13.plot( + t_factor * (t0 + t), + np.angle(np.exp(1j * full_phase_offset) * np.exp(1j * 2 * math.pi * t * full_freq_offset)), + 'b.', + markersize=0.25, + ) + + if label_defined and len(label_idx_to_plot) != 0: + arr_label_store = list(label_store.values()) + lbl_vals = np.take(arr_label_store, label_idx_to_plot) + t = t0 + adc.delay + (adc.num_samples - 1) / 2 * adc.dwell + _t = [t_factor * t] * len(lbl_vals) + # Plot each label individually to retrieve each corresponding Line2D object + p = itertools.chain.from_iterable( + [sp11.plot(__t, _lbl_vals, '.') for __t, _lbl_vals in zip(_t, lbl_vals, strict=True)] + ) + if len(label_legend_to_plot) != 0: + sp11.legend(list(p), label_legend_to_plot, loc='upper left') + label_legend_to_plot = [] + + if getattr(block, 'rf', None) is not None: # RF + rf = block.rf + time_center, index_center = calc_rf_center(rf) + time = rf.t + signal = rf.signal + + if signal.shape[0] == 2 and rf.freq_offset != 0: + num_samples = min(int(abs(rf.freq_offset)), 256) + time = np.linspace(time[0], time[-1], num_samples) + signal = np.linspace(signal[0], signal[-1], num_samples) + + if abs(signal[0]) != 0: + signal = np.concatenate(([0], signal)) + time = np.concatenate(([time[0]], time)) + index_center += 1 + + if abs(signal[-1]) != 0: + signal = np.concatenate((signal, [0])) + time = np.concatenate((time, [time[-1]])) + + signal_is_real = max(np.abs(np.imag(signal))) / max(np.abs(np.real(signal))) < 1e-6 + + full_freq_offset = rf.freq_offset + rf.freq_ppm * 1e-6 * seq.system.B0 + full_phase_offset = rf.phase_offset + rf.phase_ppm * 1e-6 * seq.system.B0 + + # If off-resonant and rectangular (2 samples), interpolate the pulse + if len(signal) == 2 and full_freq_offset != 0: + num_interp = min(int(abs(full_freq_offset)), 256) + time = np.linspace(time[0], time[-1], num_interp) + signal = np.linspace(signal[0], signal[-1], num_interp) + if abs(signal[0]) != 0: # fix strangely looking phase / amplitude in the beginning + signal = np.concatenate([[0], signal]) + time = np.concatenate([[time[0]], time]) + if abs(signal[-1]) != 0: # fix strangely looking phase / amplitude at the end + signal = np.concatenate([signal, [0]]) + time = np.concatenate([time, [time[-1]]]) + + # Compute time vector with delay applied + time_with_delay = t_factor * (t0 + time + rf.delay) + time_center_with_delay = t_factor * (t0 + time_center + rf.delay) + + # Choose plot behavior based on realness of signal + if signal_is_real: + # Plot real part of signal + sp12.plot(time_with_delay, np.real(signal)) + + # Include sign(real(signal)) factor like MATLAB + phase_corrected = ( + signal + * np.sign(np.real(signal)) + * np.exp(1j * full_phase_offset) + * np.exp(1j * 2 * math.pi * time * full_freq_offset) + ) + sc_corrected = ( + signal[index_center] + * np.exp(1j * full_phase_offset) + * np.exp(1j * 2 * math.pi * time[index_center] * full_freq_offset) + ) + + sp13.plot( + time_with_delay, + np.angle(phase_corrected), + time_center_with_delay, + np.angle(sc_corrected), + 'xb', + ) + else: + # Plot magnitude of complex signal + sp12.plot(time_with_delay, np.abs(signal)) + + # Plot angle of complex signal + phase_corrected = ( + signal * np.exp(1j * full_phase_offset) * np.exp(1j * 2 * math.pi * time * full_freq_offset) + ) + sc_corrected = ( + signal[index_center] + * np.exp(1j * full_phase_offset) + * np.exp(1j * 2 * math.pi * time[index_center] * full_freq_offset) + ) + + sp13.plot( + time_with_delay, + np.angle(phase_corrected), + time_center_with_delay, + np.angle(sc_corrected), + 'xb', + ) + + grad_channels = ['gx', 'gy', 'gz'] + for x in range(len(grad_channels)): # Gradients + if getattr(block, grad_channels[x], None) is not None: + grad = getattr(block, grad_channels[x]) + if grad.type == 'grad': + # We extend the shape by adding the first and the last points in an effort of making the + # display a bit less confusing... + time = grad.delay + np.array([0, *grad.tt, grad.shape_dur]) + waveform = g_factor * np.array((grad.first, *grad.waveform, grad.last)) + else: + time = np.array( + cumsum( + 0, + grad.delay, + grad.rise_time, + grad.flat_time, + grad.fall_time, + ) + ) + waveform = g_factor * grad.amplitude * np.array([0, 0, 1, 1, 0]) + fig2_subplots[x].plot(t_factor * (t0 + time), waveform) + t0 += seq.block_durations[block_counter] + + grad_plot_labels = ['x', 'y', 'z'] + sp11.set_ylabel('ADC') + sp12.set_ylabel('RF mag (Hz)') + sp13.set_ylabel('RF/ADC phase (rad)') + sp13.set_xlabel(f't ({time_disp})') + for x in range(3): + _label = grad_plot_labels[x] + fig2_subplots[x].set_ylabel(f'G{_label} ({grad_disp})') + fig2_subplots[-1].set_xlabel(f't ({time_disp})') + + # Setting display limits + disp_range = t_factor * np.array([time_range[0], min(t0, time_range[1])]) + [x.set_xlim(disp_range) for x in [sp11, sp12, sp13, *fig2_subplots]] + + # Grid on + for sp in [sp11, sp12, sp13, *fig2_subplots]: + sp.grid() + + fig1.tight_layout() + fig2.tight_layout() + if save: + fig1.savefig('seq_plot1.jpg') + fig2.savefig('seq_plot2.jpg') + + return fig1, fig2 diff --git a/tests/expected_output/seq1.seq b/tests/expected_output/seq1.seq index 158922f..4651f35 100644 --- a/tests/expected_output/seq1.seq +++ b/tests/expected_output/seq1.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -25,10 +25,11 @@ TotalDuration 0.0051 7 83 0 4 0 1 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 125 1 2 3 0 0 0 +1 125 1 2 3 500 0 0 0 0 0 u # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -64,4 +65,4 @@ num_samples 2 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 39ca49fddd6e69e6ebea24a67968a298 +Hash eeb227bcd331d4a060a40836ffcc638a diff --git a/tests/expected_output/seq2.seq b/tests/expected_output/seq2.seq index 34598fa..f5dd762 100644 --- a/tests/expected_output/seq2.seq +++ b/tests/expected_output/seq2.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -24,11 +24,12 @@ TotalDuration 0.0142 6 1000 0 4 0 0 1 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 250 1 2 3 0 0 0 -2 500 1 2 3 0 0 0 +1 250 1 2 3 500 0 0 0 0 0 u +2 500 1 2 3 500 0 0 0 0 0 u # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -40,10 +41,10 @@ TotalDuration 0.0142 4 102459 240 9520 240 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 100 100000 0 0 0 +1 100 100000 0 0 0 0 0 0 # Sequence Shapes [SHAPES] @@ -69,4 +70,4 @@ num_samples 2 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 208eb3220a345215abbb506cf8a3d495 +Hash 84be8c860fdbebf68977a6779ae84cb4 diff --git a/tests/expected_output/seq3.seq b/tests/expected_output/seq3.seq index d227e52..7db723a 100644 --- a/tests/expected_output/seq3.seq +++ b/tests/expected_output/seq3.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -68,10 +68,11 @@ TotalDuration 0.12722 50 1000 0 3 0 0 1 1 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 62.5 1 2 3 0 0 0 +1 62.5 1 2 3 500 0 0 0 0 0 u # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -91,10 +92,10 @@ TotalDuration 0.12722 12 1.66667e+06 240 0 240 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 100 100000 0 0 0 +1 100 100000 0 0 0 0 0 0 # Format of extension lists: # id type ref next_id @@ -132,4 +133,4 @@ num_samples 2 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 54a863dda55463e26b0cc8a849a0111b +Hash 890a2a51ec2754a0007f5e1bda944975 diff --git a/tests/expected_output/seq4.seq b/tests/expected_output/seq4.seq index 8af4668..95b432d 100644 --- a/tests/expected_output/seq4.seq +++ b/tests/expected_output/seq4.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -68,10 +68,11 @@ TotalDuration 0.12722 50 1000 0 3 0 0 1 10 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 62.5 1 2 3 0 0 0 +1 62.5 1 2 3 500 0 0 0 0 0 u # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -91,10 +92,10 @@ TotalDuration 0.12722 12 1.66667e+06 240 0 240 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 100 100000 0 0 0 +1 100 100000 0 0 0 0 0 0 # Format of extension lists: # id type ref next_id @@ -150,4 +151,4 @@ num_samples 2 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 51b6b4fac6d80f3a8ef6476f58e840ee +Hash 74a5b1a80817326ed54758f9cbe8c446 diff --git a/tests/expected_output/seq_make_block_pulses.seq b/tests/expected_output/seq_make_block_pulses.seq index 6d0e40e..6851d0f 100644 --- a/tests/expected_output/seq_make_block_pulses.seq +++ b/tests/expected_output/seq_make_block_pulses.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -31,15 +31,16 @@ TotalDuration 6.018 13 100 4 0 0 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 39.7887 1 2 3 0 0 0 -2 39.7887 1 2 3 1000 0 0 -3 62.5 1 2 3 0 0 0 -4 250 1 2 4 0 0 0 -5 125 1 2 5 0 0 1.5708 -6 250 1 2 4 0 1000 1.5708 +1 39.7887 1 2 3 2000 0 0 0 0 0 u +2 39.7887 1 2 3 2000 1000 0 0 0 0 u +3 62.5 1 2 3 2000 0 0 0 0 0 u +4 250 1 2 4 500 0 0 0 0 0 u +5 125 1 2 5 1000 0 0 0 0 1.5708 u +6 250 1 2 4 500 0 0 0 1000 1.5708 u # Sequence Shapes [SHAPES] @@ -75,4 +76,4 @@ num_samples 2 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 295982a782d6eaa751b5d3e83bcbeebe +Hash c4e19a13c587d72daa11521193fdd313 diff --git a/tests/expected_output/seq_make_gauss_pulses.seq b/tests/expected_output/seq_make_gauss_pulses.seq index 6dfa7ad..18ae7f9 100644 --- a/tests/expected_output/seq_make_gauss_pulses.seq +++ b/tests/expected_output/seq_make_gauss_pulses.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -33,17 +33,18 @@ TotalDuration 7.019 15 100 8 0 0 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 159.155 1 2 0 0 0 0 -2 159.155 1 2 0 1000 0 0 -3 250 1 2 0 0 0 0 -4 999.988 3 4 0 0 0 0 -5 499.999 5 6 0 0 0 1.5708 -6 999.988 3 4 0 0 1000 1.5708 -7 316.492 7 4 0 0 0 0 -8 1018.14 8 4 0 0 0 0 +1 159.155 1 2 0 2000 0 0 0 0 0 u +2 159.155 1 2 0 2000 1000 0 0 0 0 u +3 250 1 2 0 2000 0 0 0 0 0 u +4 999.988 3 4 0 500 0 0 0 0 0 u +5 499.999 5 6 0 1000 0 0 0 0 1.5708 u +6 999.988 3 4 0 500 0 0 0 1000 1.5708 u +7 316.492 7 4 0 500 0 0 0 0 0 u +8 1018.14 8 4 0 500 0 0 0 0 0 u # Sequence Shapes [SHAPES] @@ -9051,4 +9052,4 @@ num_samples 1000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash f295d5778ececa314823a3ba290ab1b6 +Hash 8832bf4b6933e396cf89908b426e8f06 diff --git a/tests/expected_output/seq_make_sinc_pulses.seq b/tests/expected_output/seq_make_sinc_pulses.seq index b9cb938..c2267e1 100644 --- a/tests/expected_output/seq_make_sinc_pulses.seq +++ b/tests/expected_output/seq_make_sinc_pulses.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -33,17 +33,18 @@ TotalDuration 7.019 15 100 8 0 0 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 176.286 1 2 0 0 0 0 -2 176.286 1 2 0 1000 0 0 -3 276.909 1 2 0 0 0 0 -4 1107.63 3 4 0 0 0 0 -5 553.817 5 6 0 0 0 1.5708 -6 1107.63 3 4 0 0 1000 1.5708 -7 286.482 7 8 0 0 0 0 -8 1081.31 9 4 0 0 0 0 +1 176.286 1 2 0 2000 0 0 0 0 0 u +2 176.286 1 2 0 2000 1000 0 0 0 0 u +3 276.909 1 2 0 2000 0 0 0 0 0 u +4 1107.63 3 4 0 500 0 0 0 0 0 u +5 553.817 5 6 0 1000 0 0 0 0 1.5708 u +6 1107.63 3 4 0 500 0 0 0 1000 1.5708 u +7 286.482 7 8 0 500 0 0 0 0 0 u +8 1081.31 9 4 0 500 0 0 0 0 0 u # Sequence Shapes [SHAPES] @@ -9120,4 +9121,4 @@ num_samples 1000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash ec62133c1b1a636648c233d62ef54ca2 +Hash 9019cfa5ffd98baa3847d1b085527790 diff --git a/tests/expected_output/write_epi.seq b/tests/expected_output/write_epi.seq index 78eac0c..c0dcd1a 100644 --- a/tests/expected_output/write_epi.seq +++ b/tests/expected_output/write_epi.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -408,12 +408,13 @@ TotalDuration 0.15405 390 6 0 0 6 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 329.152 1 2 0 100 -1333.33 0 -2 329.152 1 2 0 100 0 0 -3 329.152 1 2 0 100 1333.33 0 +1 329.152 1 2 0 1500 100 0 0 -1333.33 0 e +2 329.152 1 2 0 1500 100 0 0 0 0 e +3 329.152 1 2 0 1500 100 0 0 1333.33 0 e # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -428,10 +429,10 @@ TotalDuration 0.15405 7 -1.13636e+06 210 260 210 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 4000 214 0 0 +1 64 4000 214 0 0 0 0 0 # Sequence Shapes [SHAPES] @@ -3460,4 +3461,4 @@ num_samples 3000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 5492b1d869053f516a0feda6747fa5dd +Hash e6dba7be0f6647ea5e8ceb85c5afe0ac diff --git a/tests/expected_output/write_epi_label.seq b/tests/expected_output/write_epi_label.seq index a525ce3..1b2ae8a 100644 --- a/tests/expected_output/write_epi_label.seq +++ b/tests/expected_output/write_epi_label.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -5654,16 +5654,17 @@ TotalDuration 4.32868 5636 0 0 0 0 0 0 18 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 329.152 1 2 0 100 -4000 6 -2 329.152 1 2 0 100 -2666.67 4 -3 329.152 1 2 0 100 -1333.33 2 -4 329.152 1 2 0 100 0 -0 -5 329.152 1 2 0 100 1333.33 -2 -6 329.152 1 2 0 100 2666.67 -4 -7 329.152 1 2 0 100 4000 -6 +1 329.152 1 2 0 1500 100 0 0 -4000 6 e +2 329.152 1 2 0 1500 100 0 0 -2666.67 4 e +3 329.152 1 2 0 1500 100 0 0 -1333.33 2 e +4 329.152 1 2 0 1500 100 0 0 0 -0 e +5 329.152 1 2 0 1500 100 0 0 1333.33 -2 e +6 329.152 1 2 0 1500 100 0 0 2666.67 -4 e +7 329.152 1 2 0 1500 100 0 0 4000 -6 e # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -5679,10 +5680,10 @@ TotalDuration 4.32868 8 1.35307e+06 250 610 250 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 4000 214 0 0 +1 64 4000 214 0 0 0 0 0 # Format of extension lists: # id type ref next_id @@ -8762,4 +8763,4 @@ num_samples 3000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 49597aa6ec446d9a34de2d40e5b1eb25 +Hash b2950251a04185feafc64b2bb9101dca diff --git a/tests/expected_output/write_epi_se.seq b/tests/expected_output/write_epi_se.seq index f715238..e4d2632 100644 --- a/tests/expected_output/write_epi_se.seq +++ b/tests/expected_output/write_epi_se.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -154,11 +154,12 @@ TotalDuration 0.08315 136 10 0 0 0 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 329.152 1 2 0 100 0 0 -2 1000 3 4 5 100 0 0 +1 329.152 1 2 0 1500 100 0 0 0 0 e +2 1000 3 4 5 250 100 0 0 0 0 r # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -174,10 +175,10 @@ TotalDuration 0.08315 8 -781250 150 320 150 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 5000 150 0 0 +1 64 5000 150 0 0 0 0 0 # Sequence Shapes [SHAPES] @@ -3221,4 +3222,4 @@ num_samples 2 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash e3079d0c5f49e461b83b55b58a2a4aeb +Hash fa0b012993045142ae20c0ec821d5ad7 diff --git a/tests/expected_output/write_epi_se_rs.seq b/tests/expected_output/write_epi_se_rs.seq index 96e357b..ebcdd78 100644 --- a/tests/expected_output/write_epi_se_rs.seq +++ b/tests/expected_output/write_epi_se_rs.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -198,26 +198,27 @@ TotalDuration 0.21735 180 48 0 8 10 0 1 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 129.712 1 2 0 100 -424.504 0 -2 493.727 3 4 0 130 -2000 0 -3 987.454 3 4 0 1740 -2000 1.5708 -4 493.727 3 4 0 130 0 0 -5 987.454 3 4 0 1740 0 1.5708 -6 493.727 3 4 0 130 2000 0 -7 987.454 3 4 0 1740 2000 1.5708 +1 129.712 1 2 0 4000 100 0 0 -424.504 0 s +2 493.727 3 4 0 1000 130 0 0 -2000 0 e +3 987.454 3 4 0 1000 1740 0 0 -2000 1.5708 r +4 493.727 3 4 0 1000 130 0 0 0 0 e +5 987.454 3 4 0 1000 1740 0 0 0 1.5708 r +6 493.727 3 4 0 1000 130 0 0 2000 0 e +7 987.454 3 4 0 1000 1740 0 0 2000 1.5708 r # Format of arbitrary gradients: # time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster) -# id amplitude amp_shape_id time_shape_id delay -# .. Hz/m .. .. us +# id amplitude first last amp_shape_id time_shape_id delay +# .. Hz/m Hz/m Hz/m .. .. us [GRADIENTS] -3 1.34624e+06 5 6 0 -7 -133333 7 8 450 -9 -133333 9 10 0 -10 -133333 11 8 0 +3 1.34624e+06 0 0 5 6 0 +7 -133333 0 -133333 7 8 450 +9 -133333 -133333 -133333 9 10 0 +10 -133333 -133333 0 11 8 0 # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -231,10 +232,10 @@ TotalDuration 0.21735 8 -1.00036e+06 220 40 220 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 104 3900 37 0 0 +1 104 3900 37 0 0 0 0 0 # Format of extension lists: # id type ref next_id @@ -10335,4 +10336,4 @@ num_samples 2 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash a34c8de5a75831a58c38f69e0b1815fc +Hash 7754f3ee61085517dc7839965fb2b3fd diff --git a/tests/expected_output/write_gre.seq b/tests/expected_output/write_gre.seq index 0d6688a..a26c5c5 100644 --- a/tests/expected_output/write_gre.seq +++ b/tests/expected_output/write_gre.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -338,33 +338,34 @@ TotalDuration 0.768 320 378 0 6 9 8 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 37.2185 1 2 0 100 0 0 -2 37.2185 1 2 0 100 0 2.04204 -3 37.2185 1 2 0 100 0 6.12611 -4 37.2185 1 2 0 100 0 5.96903 -5 37.2185 1 2 0 100 0 1.5708 -6 37.2185 1 2 0 100 0 5.49779 -7 37.2185 1 2 0 100 0 5.18363 -8 37.2185 1 2 0 100 0 0.628319 -9 37.2185 1 2 0 100 0 4.39823 -10 37.2185 1 2 0 100 0 3.92699 -11 37.2185 1 2 0 100 0 2.82743 -12 37.2185 1 2 0 100 0 2.19911 -13 37.2185 1 2 0 100 0 3.61283 -14 37.2185 1 2 0 100 0 0.785398 -15 37.2185 1 2 0 100 0 1.25664 -16 37.2185 1 2 0 100 0 4.55531 -17 37.2185 1 2 0 100 0 4.71239 -18 37.2185 1 2 0 100 0 0.471239 -19 37.2185 1 2 0 100 0 1.41372 -20 37.2185 1 2 0 100 0 3.14159 -21 37.2185 1 2 0 100 0 5.34071 -22 37.2185 1 2 0 100 0 2.35619 -23 37.2185 1 2 0 100 0 3.76991 -24 37.2185 1 2 0 100 0 2.98451 +1 37.2185 1 2 0 1500 100 0 0 0 0 e +2 37.2185 1 2 0 1500 100 0 0 0 2.04204 e +3 37.2185 1 2 0 1500 100 0 0 0 6.12611 e +4 37.2185 1 2 0 1500 100 0 0 0 5.96903 e +5 37.2185 1 2 0 1500 100 0 0 0 1.5708 e +6 37.2185 1 2 0 1500 100 0 0 0 5.49779 e +7 37.2185 1 2 0 1500 100 0 0 0 5.18363 e +8 37.2185 1 2 0 1500 100 0 0 0 0.628319 e +9 37.2185 1 2 0 1500 100 0 0 0 4.39823 e +10 37.2185 1 2 0 1500 100 0 0 0 3.92699 e +11 37.2185 1 2 0 1500 100 0 0 0 2.82743 e +12 37.2185 1 2 0 1500 100 0 0 0 2.19911 e +13 37.2185 1 2 0 1500 100 0 0 0 3.61283 e +14 37.2185 1 2 0 1500 100 0 0 0 0.785398 e +15 37.2185 1 2 0 1500 100 0 0 0 1.25664 e +16 37.2185 1 2 0 1500 100 0 0 0 4.55531 e +17 37.2185 1 2 0 1500 100 0 0 0 4.71239 e +18 37.2185 1 2 0 1500 100 0 0 0 0.471239 e +19 37.2185 1 2 0 1500 100 0 0 0 1.41372 e +20 37.2185 1 2 0 1500 100 0 0 0 3.14159 e +21 37.2185 1 2 0 1500 100 0 0 0 5.34071 e +22 37.2185 1 2 0 1500 100 0 0 0 2.35619 e +23 37.2185 1 2 0 1500 100 0 0 0 3.76991 e +24 37.2185 1 2 0 1500 100 0 0 0 2.98451 e # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -443,33 +444,33 @@ TotalDuration 0.768 71 0 10 980 10 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 50000 20 0 0 -2 64 50000 20 0 2.04204 -3 64 50000 20 0 6.12611 -4 64 50000 20 0 5.96903 -5 64 50000 20 0 1.5708 -6 64 50000 20 0 5.49779 -7 64 50000 20 0 5.18363 -8 64 50000 20 0 0.628319 -9 64 50000 20 0 4.39823 -10 64 50000 20 0 3.92699 -11 64 50000 20 0 2.82743 -12 64 50000 20 0 2.19911 -13 64 50000 20 0 3.61283 -14 64 50000 20 0 0.785398 -15 64 50000 20 0 1.25664 -16 64 50000 20 0 4.55531 -17 64 50000 20 0 4.71239 -18 64 50000 20 0 0.471239 -19 64 50000 20 0 1.41372 -20 64 50000 20 0 3.14159 -21 64 50000 20 0 5.34071 -22 64 50000 20 0 2.35619 -23 64 50000 20 0 3.76991 -24 64 50000 20 0 2.98451 +1 64 50000 20 0 0 0 0 0 +2 64 50000 20 0 0 0 2.04204 0 +3 64 50000 20 0 0 0 6.12611 0 +4 64 50000 20 0 0 0 5.96903 0 +5 64 50000 20 0 0 0 1.5708 0 +6 64 50000 20 0 0 0 5.49779 0 +7 64 50000 20 0 0 0 5.18363 0 +8 64 50000 20 0 0 0 0.628319 0 +9 64 50000 20 0 0 0 4.39823 0 +10 64 50000 20 0 0 0 3.92699 0 +11 64 50000 20 0 0 0 2.82743 0 +12 64 50000 20 0 0 0 2.19911 0 +13 64 50000 20 0 0 0 3.61283 0 +14 64 50000 20 0 0 0 0.785398 0 +15 64 50000 20 0 0 0 1.25664 0 +16 64 50000 20 0 0 0 4.55531 0 +17 64 50000 20 0 0 0 4.71239 0 +18 64 50000 20 0 0 0 0.471239 0 +19 64 50000 20 0 0 0 1.41372 0 +20 64 50000 20 0 0 0 3.14159 0 +21 64 50000 20 0 0 0 5.34071 0 +22 64 50000 20 0 0 0 2.35619 0 +23 64 50000 20 0 0 0 3.76991 0 +24 64 50000 20 0 0 0 2.98451 0 # Sequence Shapes [SHAPES] @@ -3498,4 +3499,4 @@ num_samples 3000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 67684cb72935bc35f3efc93cd509afb0 +Hash 7160d56755113d789f083de4e9d7a9a1 diff --git a/tests/expected_output/write_gre_label.seq b/tests/expected_output/write_gre_label.seq index ec2bc82..8a33988 100644 --- a/tests/expected_output/write_gre_label.seq +++ b/tests/expected_output/write_gre_label.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -339,33 +339,34 @@ TotalDuration 0.64 321 248 0 6 9 8 0 4 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 25.6007 1 2 0 100 0 0 -2 25.6007 1 2 0 100 0 2.04204 -3 25.6007 1 2 0 100 0 6.12611 -4 25.6007 1 2 0 100 0 5.96903 -5 25.6007 1 2 0 100 0 1.5708 -6 25.6007 1 2 0 100 0 5.49779 -7 25.6007 1 2 0 100 0 5.18363 -8 25.6007 1 2 0 100 0 0.628319 -9 25.6007 1 2 0 100 0 4.39823 -10 25.6007 1 2 0 100 0 3.92699 -11 25.6007 1 2 0 100 0 2.82743 -12 25.6007 1 2 0 100 0 2.19911 -13 25.6007 1 2 0 100 0 3.61283 -14 25.6007 1 2 0 100 0 0.785398 -15 25.6007 1 2 0 100 0 1.25664 -16 25.6007 1 2 0 100 0 4.55531 -17 25.6007 1 2 0 100 0 4.71239 -18 25.6007 1 2 0 100 0 0.471239 -19 25.6007 1 2 0 100 0 1.41372 -20 25.6007 1 2 0 100 0 3.14159 -21 25.6007 1 2 0 100 0 5.34071 -22 25.6007 1 2 0 100 0 2.35619 -23 25.6007 1 2 0 100 0 3.76991 -24 25.6007 1 2 0 100 0 2.98451 +1 25.6007 1 2 0 1500 100 0 0 0 0 e +2 25.6007 1 2 0 1500 100 0 0 0 2.04204 e +3 25.6007 1 2 0 1500 100 0 0 0 6.12611 e +4 25.6007 1 2 0 1500 100 0 0 0 5.96903 e +5 25.6007 1 2 0 1500 100 0 0 0 1.5708 e +6 25.6007 1 2 0 1500 100 0 0 0 5.49779 e +7 25.6007 1 2 0 1500 100 0 0 0 5.18363 e +8 25.6007 1 2 0 1500 100 0 0 0 0.628319 e +9 25.6007 1 2 0 1500 100 0 0 0 4.39823 e +10 25.6007 1 2 0 1500 100 0 0 0 3.92699 e +11 25.6007 1 2 0 1500 100 0 0 0 2.82743 e +12 25.6007 1 2 0 1500 100 0 0 0 2.19911 e +13 25.6007 1 2 0 1500 100 0 0 0 3.61283 e +14 25.6007 1 2 0 1500 100 0 0 0 0.785398 e +15 25.6007 1 2 0 1500 100 0 0 0 1.25664 e +16 25.6007 1 2 0 1500 100 0 0 0 4.55531 e +17 25.6007 1 2 0 1500 100 0 0 0 4.71239 e +18 25.6007 1 2 0 1500 100 0 0 0 0.471239 e +19 25.6007 1 2 0 1500 100 0 0 0 1.41372 e +20 25.6007 1 2 0 1500 100 0 0 0 3.14159 e +21 25.6007 1 2 0 1500 100 0 0 0 5.34071 e +22 25.6007 1 2 0 1500 100 0 0 0 2.35619 e +23 25.6007 1 2 0 1500 100 0 0 0 3.76991 e +24 25.6007 1 2 0 1500 100 0 0 0 2.98451 e # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -444,33 +445,33 @@ TotalDuration 0.64 71 -0 10 980 10 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 50000 20 0 0 -2 64 50000 20 0 2.04204 -3 64 50000 20 0 6.12611 -4 64 50000 20 0 5.96903 -5 64 50000 20 0 1.5708 -6 64 50000 20 0 5.49779 -7 64 50000 20 0 5.18363 -8 64 50000 20 0 0.628319 -9 64 50000 20 0 4.39823 -10 64 50000 20 0 3.92699 -11 64 50000 20 0 2.82743 -12 64 50000 20 0 2.19911 -13 64 50000 20 0 3.61283 -14 64 50000 20 0 0.785398 -15 64 50000 20 0 1.25664 -16 64 50000 20 0 4.55531 -17 64 50000 20 0 4.71239 -18 64 50000 20 0 0.471239 -19 64 50000 20 0 1.41372 -20 64 50000 20 0 3.14159 -21 64 50000 20 0 5.34071 -22 64 50000 20 0 2.35619 -23 64 50000 20 0 3.76991 -24 64 50000 20 0 2.98451 +1 64 50000 20 0 0 0 0 0 +2 64 50000 20 0 0 0 2.04204 0 +3 64 50000 20 0 0 0 6.12611 0 +4 64 50000 20 0 0 0 5.96903 0 +5 64 50000 20 0 0 0 1.5708 0 +6 64 50000 20 0 0 0 5.49779 0 +7 64 50000 20 0 0 0 5.18363 0 +8 64 50000 20 0 0 0 0.628319 0 +9 64 50000 20 0 0 0 4.39823 0 +10 64 50000 20 0 0 0 3.92699 0 +11 64 50000 20 0 0 0 2.82743 0 +12 64 50000 20 0 0 0 2.19911 0 +13 64 50000 20 0 0 0 3.61283 0 +14 64 50000 20 0 0 0 0.785398 0 +15 64 50000 20 0 0 0 1.25664 0 +16 64 50000 20 0 0 0 4.55531 0 +17 64 50000 20 0 0 0 4.71239 0 +18 64 50000 20 0 0 0 0.471239 0 +19 64 50000 20 0 0 0 1.41372 0 +20 64 50000 20 0 0 0 3.14159 0 +21 64 50000 20 0 0 0 5.34071 0 +22 64 50000 20 0 0 0 2.35619 0 +23 64 50000 20 0 0 0 3.76991 0 +24 64 50000 20 0 0 0 2.98451 0 # Format of extension lists: # id type ref next_id @@ -3521,4 +3522,4 @@ num_samples 3000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash f7a94ec301ff8db005b677e444d25922 +Hash 9685a123498aa0542870622e0c1ce861 diff --git a/tests/expected_output/write_haste.seq b/tests/expected_output/write_haste.seq index fe6d6db..27cded8 100644 --- a/tests/expected_output/write_haste.seq +++ b/tests/expected_output/write_haste.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -185,26 +185,27 @@ TotalDuration 7 167 500000 0 0 0 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 394.982 3 4 0 100 0 1.5708 -2 987.454 9 10 0 100 0 0 +1 394.982 3 4 0 1250 100 0 0 0 1.5708 e +2 987.454 9 10 0 1000 100 0 0 0 0 r # Format of arbitrary gradients: # time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster) -# id amplitude amp_shape_id time_shape_id delay -# .. Hz/m .. .. us +# id amplitude first last amp_shape_id time_shape_id delay +# .. Hz/m Hz/m Hz/m .. .. us [GRADIENTS] -1 320000 1 2 0 -2 320000 5 6 0 -3 320000 7 8 0 -5 320000 5 11 0 -6 180372 12 13 0 -8 491667 14 13 0 -9 38940.8 5 15 0 -10 180372 16 13 0 -12 491667 17 13 0 +1 320000 0 320000 1 2 0 +2 320000 320000 320000 5 6 0 +3 320000 320000 320000 7 8 0 +5 320000 320000 320000 5 11 0 +6 180372 0 38940.8 12 13 0 +8 491667 320000 0 14 13 0 +9 38940.8 38940.8 38940.8 5 15 0 +10 180372 38940.8 0 16 13 0 +12 491667 0 320000 17 13 0 # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -276,10 +277,10 @@ TotalDuration 7 73 -84092.9 250 1190 250 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 100000 10 0 0 +1 64 100000 10 0 0 0 0 0 # Sequence Shapes [SHAPES] @@ -4905,4 +4906,4 @@ num_samples 4 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash ae19094c326220b5f3bba8a568014c23 +Hash af032808145852407caf689951a9cd32 diff --git a/tests/expected_output/write_mprage.seq b/tests/expected_output/write_mprage.seq index 3ace553..866ff97 100644 --- a/tests/expected_output/write_mprage.seq +++ b/tests/expected_output/write_mprage.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -5958,42 +5958,43 @@ TotalDuration 150 5940 124807 0 0 0 5 0 3 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 563.681 1 2 3 100 0 0 -2 194.444 4 5 6 610 0 0 -3 194.444 4 5 6 610 0 2.04204 -4 194.444 4 5 6 610 0 6.12611 -5 194.444 4 5 6 610 0 5.96903 -6 194.444 4 5 6 610 0 1.5708 -7 194.444 4 5 6 610 0 5.49779 -8 194.444 4 5 6 610 0 5.18363 -9 194.444 4 5 6 610 0 0.628319 -10 194.444 4 5 6 610 0 4.39823 -11 194.444 4 5 6 610 0 3.92699 -12 194.444 4 5 6 610 0 2.82743 -13 194.444 4 5 6 610 0 2.19911 -14 194.444 4 5 6 610 0 3.61283 -15 194.444 4 5 6 610 0 0.785398 -16 194.444 4 5 6 610 0 1.25664 -17 194.444 4 5 6 610 0 4.55531 -18 194.444 4 5 6 610 0 4.71239 -19 194.444 4 5 6 610 0 0.471239 -20 194.444 4 5 6 610 0 1.41372 -21 194.444 4 5 6 610 0 3.14159 -22 194.444 4 5 6 610 0 5.34071 -23 194.444 4 5 6 610 0 2.35619 -24 194.444 4 5 6 610 0 3.76991 -25 194.444 4 5 6 610 0 2.98451 +1 563.681 1 2 3 5125 100 0 0 0 0 i +2 194.444 4 5 6 50 610 0 0 0 0 e +3 194.444 4 5 6 50 610 0 0 0 2.04204 e +4 194.444 4 5 6 50 610 0 0 0 6.12611 e +5 194.444 4 5 6 50 610 0 0 0 5.96903 e +6 194.444 4 5 6 50 610 0 0 0 1.5708 e +7 194.444 4 5 6 50 610 0 0 0 5.49779 e +8 194.444 4 5 6 50 610 0 0 0 5.18363 e +9 194.444 4 5 6 50 610 0 0 0 0.628319 e +10 194.444 4 5 6 50 610 0 0 0 4.39823 e +11 194.444 4 5 6 50 610 0 0 0 3.92699 e +12 194.444 4 5 6 50 610 0 0 0 2.82743 e +13 194.444 4 5 6 50 610 0 0 0 2.19911 e +14 194.444 4 5 6 50 610 0 0 0 3.61283 e +15 194.444 4 5 6 50 610 0 0 0 0.785398 e +16 194.444 4 5 6 50 610 0 0 0 1.25664 e +17 194.444 4 5 6 50 610 0 0 0 4.55531 e +18 194.444 4 5 6 50 610 0 0 0 4.71239 e +19 194.444 4 5 6 50 610 0 0 0 0.471239 e +20 194.444 4 5 6 50 610 0 0 0 1.41372 e +21 194.444 4 5 6 50 610 0 0 0 3.14159 e +22 194.444 4 5 6 50 610 0 0 0 5.34071 e +23 194.444 4 5 6 50 610 0 0 0 2.35619 e +24 194.444 4 5 6 50 610 0 0 0 3.76991 e +25 194.444 4 5 6 50 610 0 0 0 2.98451 e # Format of arbitrary gradients: # time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster) -# id amplitude amp_shape_id time_shape_id delay -# .. Hz/m .. .. us +# id amplitude first last amp_shape_id time_shape_id delay +# .. Hz/m Hz/m Hz/m .. .. us [GRADIENTS] -4 -708063 7 8 0 -5 984720 9 10 0 +4 -708063 0 49824.6 7 8 0 +5 984720 49824.6 0 9 10 0 # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -6099,33 +6100,33 @@ TotalDuration 150 100 -23148.1 180 0 180 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 78400 380 0 0 -2 64 78400 380 0 2.04204 -3 64 78400 380 0 6.12611 -4 64 78400 380 0 5.96903 -5 64 78400 380 0 1.5708 -6 64 78400 380 0 5.49779 -7 64 78400 380 0 5.18363 -8 64 78400 380 0 0.628319 -9 64 78400 380 0 4.39823 -10 64 78400 380 0 3.92699 -11 64 78400 380 0 2.82743 -12 64 78400 380 0 2.19911 -13 64 78400 380 0 3.61283 -14 64 78400 380 0 0.785398 -15 64 78400 380 0 1.25664 -16 64 78400 380 0 4.55531 -17 64 78400 380 0 4.71239 -18 64 78400 380 0 0.471239 -19 64 78400 380 0 1.41372 -20 64 78400 380 0 3.14159 -21 64 78400 380 0 5.34071 -22 64 78400 380 0 2.35619 -23 64 78400 380 0 3.76991 -24 64 78400 380 0 2.98451 +1 64 78400 380 0 0 0 0 0 +2 64 78400 380 0 0 0 2.04204 0 +3 64 78400 380 0 0 0 6.12611 0 +4 64 78400 380 0 0 0 5.96903 0 +5 64 78400 380 0 0 0 1.5708 0 +6 64 78400 380 0 0 0 5.49779 0 +7 64 78400 380 0 0 0 5.18363 0 +8 64 78400 380 0 0 0 0.628319 0 +9 64 78400 380 0 0 0 4.39823 0 +10 64 78400 380 0 0 0 3.92699 0 +11 64 78400 380 0 0 0 2.82743 0 +12 64 78400 380 0 0 0 2.19911 0 +13 64 78400 380 0 0 0 3.61283 0 +14 64 78400 380 0 0 0 0.785398 0 +15 64 78400 380 0 0 0 1.25664 0 +16 64 78400 380 0 0 0 4.55531 0 +17 64 78400 380 0 0 0 4.71239 0 +18 64 78400 380 0 0 0 0.471239 0 +19 64 78400 380 0 0 0 1.41372 0 +20 64 78400 380 0 0 0 3.14159 0 +21 64 78400 380 0 0 0 5.34071 0 +22 64 78400 380 0 0 0 2.35619 0 +23 64 78400 380 0 0 0 3.76991 0 +24 64 78400 380 0 0 0 2.98451 0 # Format of extension lists: # id type ref next_id @@ -8262,4 +8263,4 @@ num_samples 4 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash 8a2e143cc8c92f04e1d6676612fcac28 +Hash aa9547958f21808d0bc90c8b0ba90fda diff --git a/tests/expected_output/write_radial_gre.seq b/tests/expected_output/write_radial_gre.seq index b801f9c..d1d29e9 100644 --- a/tests/expected_output/write_radial_gre.seq +++ b/tests/expected_output/write_radial_gre.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -423,33 +423,34 @@ TotalDuration 1.62081 405 923 0 180 135 7 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 27.4293 1 2 0 100 0 0 -2 27.4293 1 2 0 100 0 2.04204 -3 27.4293 1 2 0 100 0 6.12611 -4 27.4293 1 2 0 100 0 5.96903 -5 27.4293 1 2 0 100 0 1.5708 -6 27.4293 1 2 0 100 0 5.49779 -7 27.4293 1 2 0 100 0 5.18363 -8 27.4293 1 2 0 100 0 0.628319 -9 27.4293 1 2 0 100 0 4.39823 -10 27.4293 1 2 0 100 0 3.92699 -11 27.4293 1 2 0 100 0 2.82743 -12 27.4293 1 2 0 100 0 2.19911 -13 27.4293 1 2 0 100 0 3.61283 -14 27.4293 1 2 0 100 0 0.785398 -15 27.4293 1 2 0 100 0 1.25664 -16 27.4293 1 2 0 100 0 4.55531 -17 27.4293 1 2 0 100 0 4.71239 -18 27.4293 1 2 0 100 0 0.471239 -19 27.4293 1 2 0 100 0 1.41372 -20 27.4293 1 2 0 100 0 3.14159 -21 27.4293 1 2 0 100 0 5.34071 -22 27.4293 1 2 0 100 0 2.35619 -23 27.4293 1 2 0 100 0 3.76991 -24 27.4293 1 2 0 100 0 2.98451 +1 27.4293 1 2 0 2000 100 0 0 0 0 e +2 27.4293 1 2 0 2000 100 0 0 0 2.04204 e +3 27.4293 1 2 0 2000 100 0 0 0 6.12611 e +4 27.4293 1 2 0 2000 100 0 0 0 5.96903 e +5 27.4293 1 2 0 2000 100 0 0 0 1.5708 e +6 27.4293 1 2 0 2000 100 0 0 0 5.49779 e +7 27.4293 1 2 0 2000 100 0 0 0 5.18363 e +8 27.4293 1 2 0 2000 100 0 0 0 0.628319 e +9 27.4293 1 2 0 2000 100 0 0 0 4.39823 e +10 27.4293 1 2 0 2000 100 0 0 0 3.92699 e +11 27.4293 1 2 0 2000 100 0 0 0 2.82743 e +12 27.4293 1 2 0 2000 100 0 0 0 2.19911 e +13 27.4293 1 2 0 2000 100 0 0 0 3.61283 e +14 27.4293 1 2 0 2000 100 0 0 0 0.785398 e +15 27.4293 1 2 0 2000 100 0 0 0 1.25664 e +16 27.4293 1 2 0 2000 100 0 0 0 4.55531 e +17 27.4293 1 2 0 2000 100 0 0 0 4.71239 e +18 27.4293 1 2 0 2000 100 0 0 0 0.471239 e +19 27.4293 1 2 0 2000 100 0 0 0 1.41372 e +20 27.4293 1 2 0 2000 100 0 0 0 3.14159 e +21 27.4293 1 2 0 2000 100 0 0 0 5.34071 e +22 27.4293 1 2 0 2000 100 0 0 0 2.35619 e +23 27.4293 1 2 0 2000 100 0 0 0 3.76991 e +24 27.4293 1 2 0 2000 100 0 0 0 2.98451 e # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -637,33 +638,33 @@ TotalDuration 1.62081 180 -768177 160 0 160 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 20000 40 0 0.471239 -2 64 20000 40 0 1.41372 -3 64 20000 40 0 4.39823 -4 64 20000 40 0 3.14159 -5 64 20000 40 0 3.92699 -6 64 20000 40 0 5.34071 -7 64 20000 40 0 5.96903 -8 64 20000 40 0 2.35619 -9 64 20000 40 0 0.785398 -10 64 20000 40 0 1.25664 -11 64 20000 40 0 3.76991 -12 64 20000 40 0 2.04204 -13 64 20000 40 0 4.71239 -14 64 20000 40 0 2.82743 -15 64 20000 40 0 2.98451 -16 64 20000 40 0 5.18363 -17 64 20000 40 0 1.5708 -18 64 20000 40 0 3.61283 -19 64 20000 40 0 4.55531 -20 64 20000 40 0 0 -21 64 20000 40 0 2.19911 -22 64 20000 40 0 5.49779 -23 64 20000 40 0 0.628319 -24 64 20000 40 0 6.12611 +1 64 20000 40 0 0 0 0.471239 0 +2 64 20000 40 0 0 0 1.41372 0 +3 64 20000 40 0 0 0 4.39823 0 +4 64 20000 40 0 0 0 3.14159 0 +5 64 20000 40 0 0 0 3.92699 0 +6 64 20000 40 0 0 0 5.34071 0 +7 64 20000 40 0 0 0 5.96903 0 +8 64 20000 40 0 0 0 2.35619 0 +9 64 20000 40 0 0 0 0.785398 0 +10 64 20000 40 0 0 0 1.25664 0 +11 64 20000 40 0 0 0 3.76991 0 +12 64 20000 40 0 0 0 2.04204 0 +13 64 20000 40 0 0 0 4.71239 0 +14 64 20000 40 0 0 0 2.82743 0 +15 64 20000 40 0 0 0 2.98451 0 +16 64 20000 40 0 0 0 5.18363 0 +17 64 20000 40 0 0 0 1.5708 0 +18 64 20000 40 0 0 0 3.61283 0 +19 64 20000 40 0 0 0 4.55531 0 +20 64 20000 40 0 0 0 0 0 +21 64 20000 40 0 0 0 2.19911 0 +22 64 20000 40 0 0 0 5.49779 0 +23 64 20000 40 0 0 0 0.628319 0 +24 64 20000 40 0 0 0 6.12611 0 # Sequence Shapes [SHAPES] @@ -4692,4 +4693,4 @@ num_samples 4000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash bc6a5dca3453696e2a15de3cb78e83a0 +Hash 58cdbaad4e1e8082860b71b9826935db diff --git a/tests/expected_output/write_tse.seq b/tests/expected_output/write_tse.seq index 351b2e6..3704c86 100644 --- a/tests/expected_output/write_tse.seq +++ b/tests/expected_output/write_tse.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -368,26 +368,27 @@ TotalDuration 10 350 179761 0 0 0 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 394.982 3 4 0 100 0 1.5708 -2 987.454 9 10 0 100 0 0 +1 394.982 3 4 0 1250 100 0 0 0 1.5708 e +2 987.454 9 10 0 1000 100 0 0 0 0 r # Format of arbitrary gradients: # time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster) -# id amplitude amp_shape_id time_shape_id delay -# .. Hz/m .. .. us +# id amplitude first last amp_shape_id time_shape_id delay +# .. Hz/m Hz/m Hz/m .. .. us [GRADIENTS] -1 320000 1 2 0 -2 320000 5 6 0 -3 320000 7 8 0 -5 320000 5 11 0 -6 180372 12 13 0 -8 491667 14 13 0 -9 38940.8 5 15 0 -10 180372 16 13 0 -11 491667 17 13 0 +1 320000 0 320000 1 2 0 +2 320000 320000 320000 5 6 0 +3 320000 320000 320000 7 8 0 +5 320000 320000 320000 5 11 0 +6 180372 0 38940.8 12 13 0 +8 491667 320000 0 14 13 0 +9 38940.8 38940.8 38940.8 5 15 0 +10 180372 38940.8 0 16 13 0 +11 491667 0 320000 17 13 0 # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -461,10 +462,10 @@ TotalDuration 10 75 86805.6 250 1190 250 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 64 100000 10 0 0 +1 64 100000 10 0 0 0 0 0 # Sequence Shapes [SHAPES] @@ -5090,4 +5091,4 @@ num_samples 4 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash b17ad152434c4e0e5d3dd6f499fd02e5 +Hash 143ed0af8842a255835442c1ba7df1e7 diff --git a/tests/expected_output/write_ute.seq b/tests/expected_output/write_ute.seq index dc9b911..03bddf2 100644 --- a/tests/expected_output/write_ute.seq +++ b/tests/expected_output/write_ute.seq @@ -3,8 +3,8 @@ [VERSION] major 1 -minor 4 -revision 2 +minor 5 +revision 0 [DEFINITIONS] AdcRasterTime 1e-07 @@ -274,33 +274,34 @@ TotalDuration 0.64 256 584 0 15 37 0 0 0 # Format of RF events: -# id amplitude mag_id phase_id time_shape_id delay freq phase -# .. Hz .... .... .... us Hz rad +# id ampl. mag_id phase_id time_shape_id center delay freqPPm phasePPM freq phase use +# .. Hz .. .. .. us us ppm rad/MHz Hz rad .. +# Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined [RF] -1 161.288 1 2 0 160 0 0 -2 161.288 1 2 0 160 0 2.04204 -3 161.288 1 2 0 160 0 6.12611 -4 161.288 1 2 0 160 0 5.96903 -5 161.288 1 2 0 160 0 1.5708 -6 161.288 1 2 0 160 0 5.49779 -7 161.288 1 2 0 160 0 5.18363 -8 161.288 1 2 0 160 0 0.628319 -9 161.288 1 2 0 160 0 4.39823 -10 161.288 1 2 0 160 0 3.92699 -11 161.288 1 2 0 160 0 2.82743 -12 161.288 1 2 0 160 0 2.19911 -13 161.288 1 2 0 160 0 3.61283 -14 161.288 1 2 0 160 0 0.785398 -15 161.288 1 2 0 160 0 1.25664 -16 161.288 1 2 0 160 0 4.55531 -17 161.288 1 2 0 160 0 4.71239 -18 161.288 1 2 0 160 0 0.471239 -19 161.288 1 2 0 160 0 1.41372 -20 161.288 1 2 0 160 0 3.14159 -21 161.288 1 2 0 160 0 5.34071 -22 161.288 1 2 0 160 0 2.35619 -23 161.288 1 2 0 160 0 3.76991 -24 161.288 1 2 0 160 0 2.98451 +1 161.288 1 2 0 1000 160 0 0 0 0 e +2 161.288 1 2 0 1000 160 0 0 0 2.04204 e +3 161.288 1 2 0 1000 160 0 0 0 6.12611 e +4 161.288 1 2 0 1000 160 0 0 0 5.96903 e +5 161.288 1 2 0 1000 160 0 0 0 1.5708 e +6 161.288 1 2 0 1000 160 0 0 0 5.49779 e +7 161.288 1 2 0 1000 160 0 0 0 5.18363 e +8 161.288 1 2 0 1000 160 0 0 0 0.628319 e +9 161.288 1 2 0 1000 160 0 0 0 4.39823 e +10 161.288 1 2 0 1000 160 0 0 0 3.92699 e +11 161.288 1 2 0 1000 160 0 0 0 2.82743 e +12 161.288 1 2 0 1000 160 0 0 0 2.19911 e +13 161.288 1 2 0 1000 160 0 0 0 3.61283 e +14 161.288 1 2 0 1000 160 0 0 0 0.785398 e +15 161.288 1 2 0 1000 160 0 0 0 1.25664 e +16 161.288 1 2 0 1000 160 0 0 0 4.55531 e +17 161.288 1 2 0 1000 160 0 0 0 4.71239 e +18 161.288 1 2 0 1000 160 0 0 0 0.471239 e +19 161.288 1 2 0 1000 160 0 0 0 1.41372 e +20 161.288 1 2 0 1000 160 0 0 0 3.14159 e +21 161.288 1 2 0 1000 160 0 0 0 5.34071 e +22 161.288 1 2 0 1000 160 0 0 0 2.35619 e +23 161.288 1 2 0 1000 160 0 0 0 3.76991 e +24 161.288 1 2 0 1000 160 0 0 0 2.98451 e # Format of trapezoid gradients: # id amplitude rise flat fall delay @@ -372,33 +373,33 @@ TotalDuration 0.64 64 -5.8783e-11 80 0 80 0 # Format of ADC events: -# id num dwell delay freq phase -# .. .. ns us Hz rad +# id num dwell delay freqPPM phasePPM freq phase phase_id +# .. .. ns us ppm rad/MHz Hz rad .. [ADC] -1 128 20000 20 0 0 -2 128 20000 20 0 2.04204 -3 128 20000 20 0 6.12611 -4 128 20000 20 0 5.96903 -5 128 20000 20 0 1.5708 -6 128 20000 20 0 5.49779 -7 128 20000 20 0 5.18363 -8 128 20000 20 0 0.628319 -9 128 20000 20 0 4.39823 -10 128 20000 20 0 3.92699 -11 128 20000 20 0 2.82743 -12 128 20000 20 0 2.19911 -13 128 20000 20 0 3.61283 -14 128 20000 20 0 0.785398 -15 128 20000 20 0 1.25664 -16 128 20000 20 0 4.55531 -17 128 20000 20 0 4.71239 -18 128 20000 20 0 0.471239 -19 128 20000 20 0 1.41372 -20 128 20000 20 0 3.14159 -21 128 20000 20 0 5.34071 -22 128 20000 20 0 2.35619 -23 128 20000 20 0 3.76991 -24 128 20000 20 0 2.98451 +1 128 20000 20 0 0 0 0 0 +2 128 20000 20 0 0 0 2.04204 0 +3 128 20000 20 0 0 0 6.12611 0 +4 128 20000 20 0 0 0 5.96903 0 +5 128 20000 20 0 0 0 1.5708 0 +6 128 20000 20 0 0 0 5.49779 0 +7 128 20000 20 0 0 0 5.18363 0 +8 128 20000 20 0 0 0 0.628319 0 +9 128 20000 20 0 0 0 4.39823 0 +10 128 20000 20 0 0 0 3.92699 0 +11 128 20000 20 0 0 0 2.82743 0 +12 128 20000 20 0 0 0 2.19911 0 +13 128 20000 20 0 0 0 3.61283 0 +14 128 20000 20 0 0 0 0.785398 0 +15 128 20000 20 0 0 0 1.25664 0 +16 128 20000 20 0 0 0 4.55531 0 +17 128 20000 20 0 0 0 4.71239 0 +18 128 20000 20 0 0 0 0.471239 0 +19 128 20000 20 0 0 0 1.41372 0 +20 128 20000 20 0 0 0 3.14159 0 +21 128 20000 20 0 0 0 5.34071 0 +22 128 20000 20 0 0 0 2.35619 0 +23 128 20000 20 0 0 0 3.76991 0 +24 128 20000 20 0 0 0 2.98451 0 # Sequence Shapes [SHAPES] @@ -1423,4 +1424,4 @@ num_samples 1000 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE] # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification) Type md5 -Hash b46103f50a4aad91cdfe7c3da4935da9 +Hash d068b005279832626cd595b50b51baea diff --git a/tests/test_check_timing.py b/tests/test_check_timing.py index c640b62..f6747d5 100644 --- a/tests/test_check_timing.py +++ b/tests/test_check_timing.py @@ -54,7 +54,7 @@ def test_check_timing(): seq.add_block(adc) # Block 4: RASTER adc = pp.make_adc(num_samples=100, duration=1e-3, system=system_broken) - seq.add_block(adc) # Block 5: ADC_DEAD_TIME, POST_ADC_DEAD_TIME + seq.add_block(adc) # Block 5: ADC_DEAD_TIME, POST_ADC_DEAD_TIME, BLOCK_DURATION_MISMATCH gx = pp.make_trapezoid(channel='x', area=1, duration=1, system=system) seq.add_block(gx) # Block 6: No error @@ -84,6 +84,9 @@ def test_check_timing(): assert exists_in_error_report(error_report, 5, event='adc', field='delay', error_type='ADC_DEAD_TIME') assert exists_in_error_report(error_report, 5, event='adc', field='duration', error_type='POST_ADC_DEAD_TIME') + assert exists_in_error_report( + error_report, 5, event='block', field='duration', error_type='BLOCK_DURATION_MISMATCH' + ) assert exists_in_error_report(error_report, 7, event='block', field='duration', error_type='RASTER') assert exists_in_error_report(error_report, 7, event='gx', field='flat_time', error_type='RASTER') @@ -94,4 +97,4 @@ def test_check_timing(): assert exists_in_error_report(error_report, 9, event='gx', field='delay', error_type='NEGATIVE_DELAY') - assert len(error_report) == 12, 'Total number of timing errors was expected to be 12' + assert len(error_report) == 13, 'Total number of timing errors was expected to be 12' diff --git a/tests/test_scale_grad.py b/tests/test_scale_grad.py new file mode 100644 index 0000000..7d98170 --- /dev/null +++ b/tests/test_scale_grad.py @@ -0,0 +1,132 @@ +import numpy as np +import pypulseq as pp +import pytest +from pypulseq import scale_grad +from pypulseq.opts import Opts + +# Updated gradients with realistic hardware limits +grad_list = [ + pp.make_trapezoid(channel='x', amplitude=10, duration=13, max_grad=30, max_slew=200), + pp.make_trapezoid(channel='y', amplitude=10, duration=13, max_grad=30, max_slew=200), + pp.make_trapezoid(channel='z', amplitude=10, duration=13, max_grad=30, max_slew=200), + pp.make_trapezoid(channel='x', amplitude=20, duration=5, max_grad=25, max_slew=150), + pp.make_trapezoid(channel='y', amplitude=20, duration=5, max_grad=25, max_slew=150), + pp.make_trapezoid(channel='z', amplitude=20, duration=5, max_grad=25, max_slew=150), + pp.make_extended_trapezoid( + 'x', [0, 15, 5, 10], convert_to_arbitrary=True, times=[1, 3, 4, 7], max_grad=40, max_slew=300 + ), + pp.make_extended_trapezoid( + 'y', [0, 15, 5, 10], convert_to_arbitrary=True, times=[1, 3, 4, 7], max_grad=40, max_slew=300 + ), + pp.make_extended_trapezoid( + 'z', [0, 15, 5, 10], convert_to_arbitrary=True, times=[1, 3, 4, 7], max_grad=40, max_slew=300 + ), + pp.make_extended_trapezoid( + 'x', [0, 20, 10, 15], convert_to_arbitrary=False, times=[1, 3, 4, 7], max_grad=25, max_slew=150 + ), + pp.make_extended_trapezoid( + 'y', [0, 20, 10, 15], convert_to_arbitrary=False, times=[1, 3, 4, 7], max_grad=25, max_slew=150 + ), + pp.make_extended_trapezoid( + 'z', [0, 20, 10, 15], convert_to_arbitrary=False, times=[1, 3, 4, 7], max_grad=25, max_slew=150 + ), + pp.make_extended_trapezoid( + 'x', [0, 10, 5, 10], convert_to_arbitrary=False, times=[1, 2, 3, 4], max_grad=15, max_slew=80 + ), + pp.make_extended_trapezoid( + 'y', [0, 10, 5, 10], convert_to_arbitrary=False, times=[1, 2, 3, 4], max_grad=15, max_slew=80 + ), + pp.make_extended_trapezoid( + 'z', [0, 10, 5, 10], convert_to_arbitrary=False, times=[1, 2, 3, 4], max_grad=15, max_slew=80 + ), +] + +# ----------- TEST SCALING IS CORRECT ----------- + + +@pytest.mark.parametrize('grad', grad_list) +def test_scale_grad_correct_scaling(grad): + scale = 0.5 # Safe scale + system = Opts(max_grad=40, max_slew=300) + scaled = scale_grad(grad, scale, system) + + if grad.type == 'trap': + assert np.isclose(scaled.amplitude, grad.amplitude * scale) + assert np.isclose(scaled.flat_area, grad.flat_area * scale) + else: + assert np.allclose(scaled.waveform, grad.waveform * scale) + assert np.isclose(scaled.first, grad.first * scale) + assert np.isclose(scaled.last, grad.last * scale) + + assert np.isclose(scaled.area, grad.area * scale) + + # ID must be removed + assert not hasattr(scaled, 'id') + + +# ----------- TEST AMPLITUDE VIOLATION ----------- + + +def test_scale_grad_violates_amplitude(): + scale = 100.0 + system = Opts(max_grad=40, max_slew=999999999) # make sure we have infinite slew rate + + expected_failures = 0 + actual_failures = 0 + + for grad in grad_list: + if grad.type == 'trap': + should_fail = abs(grad.amplitude) * scale > system.max_grad + else: + should_fail = max(abs(grad.waveform)) * scale > system.max_grad + + if should_fail: + expected_failures += 1 + with pytest.raises(ValueError, match='maximum amplitude exceeded'): + scale_grad(grad, scale, system) + actual_failures += 1 + else: + scale_grad(grad, scale, system) + + assert expected_failures == actual_failures, ( + f'Expected {expected_failures} amplitude violations, but got {actual_failures}.' + ) + + +# ----------- TEST SLEW RATE VIOLATION ----------- + + +def test_scale_grad_violates_slew(): + scale = 100.0 + system = Opts(max_grad=999999999, max_slew=300) # make sure we have infinite gradient amp + + expected_failures = 0 + actual_failures = 0 + + for grad in grad_list: + if grad.type == 'trap': + if abs(grad.amplitude) > 1e-6: + approx_slew = abs(grad.amplitude * scale) / min(grad.rise_time, grad.fall_time) + should_fail = approx_slew > system.max_slew + else: + should_fail = False + else: + if max(abs(grad.waveform)) > 1e-6: + waveform = np.array(grad.waveform) * scale + tt = np.array(grad.tt) + diffs = np.abs(np.diff(waveform) / np.diff(tt)) + should_fail = max(diffs) > system.max_slew + else: + should_fail = False + + if should_fail: + expected_failures += 1 + with pytest.raises(ValueError, match='maximum slew rate exceeded'): + scale_grad(grad, scale, system) + actual_failures += 1 + else: + scale_grad(grad, scale, system) + + assert expected_failures == actual_failures, ( + f'Expected {expected_failures} slew violations, but got {actual_failures}.' + ) diff --git a/tests/test_sequence.py b/tests/test_sequence.py index 2e0c2f4..4824146 100644 --- a/tests/test_sequence.py +++ b/tests/test_sequence.py @@ -265,12 +265,12 @@ def seq6(): 'write_gre_label_softdelay', 'write_haste', 'write_radial_gre', - 'write_tse', + 'write_tse', # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last') 'write_epi', 'write_epi_label', 'write_epi_se', - 'write_epi_se_rs', - 'write_mprage', + 'write_epi_se_rs', # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last') + 'write_mprage', # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last') 'write_ute', ] @@ -365,11 +365,11 @@ def test_writeread(self, seq_func, tmp_path, compare_seq_file): block_orig = seq.get_block(block_counter) block_compare = seq2.get_block(block_counter) - if hasattr(block_orig, 'rf') and hasattr(block_orig.rf, 'use'): - from copy import deepcopy + # if hasattr(block_orig, 'rf') and hasattr(block_orig.rf, 'use'): + # from copy import deepcopy - block_orig = deepcopy(block_orig) - block_orig.rf.use = 'undefined' + # block_orig = deepcopy(block_orig) + # block_orig.rf.use = 'undefined' assert block_compare == Approx(block_orig, abs=1e-5, rel=1e-5), f'Block {block_counter} does not match' @@ -393,9 +393,9 @@ def test_writeread(self, seq_func, tmp_path, compare_seq_file): # Restore RF use for k-space calculation for block_counter in seq.block_events: block_orig = seq.get_block(block_counter) - if hasattr(block_orig, 'rf') and hasattr(block_orig.rf, 'use'): - block_compare = seq2.get_block(block_counter) - block_compare.rf.use = block_orig.rf.use + # if hasattr(block_orig, 'rf') and hasattr(block_orig.rf, 'use'): + # block_compare = seq2.get_block(block_counter) + # block_compare.rf.use = block_orig.rf.use # Test for approximate equality of kspace calculation assert seq2.calculate_kspace() == Approx(seq.calculate_kspace(), abs=1e-1, nan_ok=True)