diff --git a/files.txt b/files.txt new file mode 100644 index 0000000..dfe9733 --- /dev/null +++ b/files.txt @@ -0,0 +1,84 @@ +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/EGG-INFO/dependency_links.txt +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/EGG-INFO/PKG-INFO +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/EGG-INFO/top_level.txt +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/EGG-INFO/SOURCES.txt +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/EGG-INFO/not-zip-safe +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/EGG-INFO/requires.txt +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/split_gradient.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_gauss_pulse.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_block_pulse.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/add_ramps.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_extended_trapezoid.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/calc_duration.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/opts.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_delay.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/calc_rf_center.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/traj_to_grad.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/convert.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_trap_pulse.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/calc_ramp.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/add_gradients.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__init__.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/check_timing.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/align.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/event_lib.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/split_gradient_at.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/decompress_shape.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_sinc_pulse.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/compress_shape.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_adc.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_arbitrary_grad.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/make_arbitrary_rf.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/points_to_waveform.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/block.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/write_seq.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/test_report.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/sequence.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/__init__.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/read_seq.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/__pycache__/sequence.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/__pycache__/read_seq.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/__pycache__/block.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/__pycache__/__init__.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/__pycache__/test_report.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/Sequence/__pycache__/write_seq.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_trap_pulse.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_delay.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/check_timing.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_arbitrary_rf.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_extended_trapezoid.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_block_pulse.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/compress_shape.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/calc_ramp.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/convert.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_gauss_pulse.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/__init__.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_arbitrary_grad.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/calc_duration.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/decompress_shape.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/align.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/opts.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/points_to_waveform.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/event_lib.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/traj_to_grad.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_adc.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/split_gradient_at.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/add_ramps.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/split_gradient.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/add_gradients.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/make_sinc_pulse.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/__pycache__/calc_rf_center.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/vds_2d.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/k2g.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/__init__.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/SAR/__init__.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/SAR/QGlobal.mat +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/SAR/SAR_calc.py +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/SAR/__pycache__/SAR_calc.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/SAR/__pycache__/__init__.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/__pycache__/__init__.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/__pycache__/vds_2d.cpython-37.pyc +/home/veldmannm/miniconda3/lib/python3.7/site-packages/pypulseq-1.2.0.post3-py3.7.egg/pypulseq/utils/__pycache__/k2g.cpython-37.pyc +/home/veldmannm/miniconda3/bin/f2py +/home/veldmannm/miniconda3/bin/f2py3 +/home/veldmannm/miniconda3/bin/f2py3.7 diff --git a/pypulseq/.vscode/settings.json b/pypulseq/.vscode/settings.json new file mode 100644 index 0000000..e72dbb0 --- /dev/null +++ b/pypulseq/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "python.pythonPath": "/home/veldmannm/miniconda3/envs/spiral/bin/python" +} \ No newline at end of file diff --git a/pypulseq/Sequence/block.py b/pypulseq/Sequence/block.py index 8268870..141d8c3 100755 --- a/pypulseq/Sequence/block.py +++ b/pypulseq/Sequence/block.py @@ -1,224 +1,242 @@ from types import SimpleNamespace import numpy as np from pypulseq.calc_duration import calc_duration from pypulseq.compress_shape import compress_shape from pypulseq.decompress_shape import decompress_shape def add_block(self, block_index: int, *args): """ Inserts pulse sequence events into `self.block_events` at position `block_index`. Also performs gradient checks. Parameters ---------- block_index : int Index at which `SimpleNamespace` events have to be inserted into `self.block_events`. args : list List of `SimpleNamespace` pulse sequence events to be added to `self.block_events`. """ block_duration = calc_duration(*args) - self.block_events[block_index] = np.zeros(6, dtype=np.int) + self.block_events[block_index] = np.zeros(7, dtype=np.int) # changed to 7 entries for trigger support by mveldmann duration = 0 check_g = {} for event in args: if event.type == 'rf': mag = np.abs(event.signal) amplitude = np.max(mag) mag = np.divide(mag, amplitude) # Following line of code is a workaround for numpy's divide functions returning NaN when mathematical # edge cases are encountered (eg. divide by 0) mag[np.isnan(mag)] = 0 phase = np.angle(event.signal) phase[phase < 0] += 2 * np.pi phase /= 2 * np.pi mag_shape = compress_shape(mag) data = np.insert(mag_shape.data, 0, mag_shape.num_samples) mag_id, found = self.shape_library.find(data) if not found: self.shape_library.insert(mag_id, data) phase_shape = compress_shape(phase) data = np.insert(phase_shape.data, 0, phase_shape.num_samples) phase_id, found = self.shape_library.find(data) if not found: self.shape_library.insert(phase_id, data) use = 0 use_cases = {'excitation': 1, 'refocusing': 2, 'inversion': 3} if hasattr(event, 'use'): use = use_cases[event.use] data = np.array((amplitude, mag_id, phase_id, event.delay, event.freq_offset, event.phase_offset, event.dead_time, event.ringdown_time, use)) data_id, found = self.rf_library.find(data) if not found: self.rf_library.insert(data_id, data) self.block_events[block_index][1] = data_id duration = max(duration, max( mag.shape) * self.rf_raster_time + event.dead_time + event.ringdown_time + event.delay) elif event.type == 'grad': channel_num = ['x', 'y', 'z'].index(event.channel) idx = 2 + channel_num check_g[channel_num] = SimpleNamespace() check_g[channel_num].idx = idx check_g[channel_num].start = np.array((event.delay + min(event.t), event.first)) check_g[channel_num].stop = np.array( (event.delay + max(event.t) + self.system.grad_raster_time, event.last)) amplitude = max(abs(event.waveform)) if amplitude > 0: g = event.waveform / amplitude else: g = event.waveform shape = compress_shape(g) data = np.insert(shape.data, 0, shape.num_samples) shape_id, found = self.shape_library.find(data) if not found: self.shape_library.insert(shape_id, data) data = np.array([amplitude, shape_id, event.delay, event.first, event.last]) grad_id, found = self.grad_library.find(data) if not found: self.grad_library.insert(grad_id, data, 'g') idx = 2 + channel_num self.block_events[block_index][idx] = grad_id duration = max(duration, len(g) * self.grad_raster_time) elif event.type == 'trap': channel_num = ['x', 'y', 'z'].index(event.channel) idx = 2 + channel_num check_g[channel_num] = SimpleNamespace() check_g[channel_num].idx = idx check_g[channel_num].start = np.array((0, 0)) check_g[channel_num].stop = np.array((event.delay + event.rise_time + event.fall_time + event.flat_time, 0)) data = np.array([event.amplitude, event.rise_time, event.flat_time, event.fall_time, event.delay]) trap_id, found = self.grad_library.find(data) if not found: self.grad_library.insert(trap_id, data, 't') self.block_events[block_index][idx] = trap_id duration = max(duration, event.delay + event.rise_time + event.flat_time + event.fall_time) elif event.type == 'adc': data = np.array( [event.num_samples, event.dwell, max(event.delay, event.dead_time), event.freq_offset, event.phase_offset, event.dead_time]) adc_id, found = self.adc_library.find(data) if not found: self.adc_library.insert(adc_id, data) self.block_events[block_index][5] = adc_id duration = max(duration, event.delay + event.num_samples * event.dwell + event.dead_time) elif event.type == 'delay': data = np.array([event.delay]) delay_id, found = self.delay_library.find(data) if not found: self.delay_library.insert(delay_id, data) self.block_events[block_index][0] = delay_id duration = max(duration, event.delay) + # inserted for trigger support by mveldmann + elif event.type == 'trigger': + data = np.array([event.delay, event.duration]) + extension_id, found = self.extension_library.find(data) + if not found: + self.extension_library.insert(extension_id, data) + self.block_events[block_index][6] = extension_id + duration = max(duration, event.delay + event.duration) + # ========= # PERFORM GRADIENT CHECKS # ========= for cg_temp in check_g.keys(): cg = check_g[cg_temp] if abs(cg.start[1]) > self.system.max_slew * self.system.grad_raster_time: if cg.start[0] != 0: raise ValueError('No delay allowed for gradients which start with a non-zero amplitude') if block_index > 1: prev_id = self.block_events[block_index - 1][cg.idx] if prev_id != 0: prev_lib = self.grad_library.get(prev_id) prev_dat = prev_lib['data'] prev_type = prev_lib['type'] if prev_type == 't': raise Exception( 'Two consecutive gradients need to have the same amplitude at the connection point') elif prev_type == 'g': last = prev_dat[4] if abs(last - cg.start[1]) > self.system.max_slew * self.system.grad_raster_time: raise Exception( 'Two consecutive gradients need to have the same amplitude at the connection point') else: raise Exception('First gradient in the the first block has to start at 0.') if cg.stop[1] > self.system.max_slew * self.system.grad_raster_time and abs(cg.stop[0] - block_duration) > 1e-7: raise Exception('A gradient that doesnt end at zero needs to be aligned to the block boundary.') def get_block(self, block_index: int) -> SimpleNamespace: """ Returns Pulseq block at `block_index` position in `self.block_events`. Parameters ---------- block_index : int Index of Pulseq block to be retrieved from `self.block_events`. Returns ------- block : SimpleNamespace Pulseq block at 'block_index' position in `self.block_events`. """ block = SimpleNamespace() event_ind = self.block_events[block_index] if event_ind[0] > 0: delay = SimpleNamespace() delay.type = 'delay' delay.delay = self.delay_library.data[event_ind[0]][0] block.delay = delay elif event_ind[1] > 0: block.rf = self.rf_from_lib_data(self.rf_library.data[event_ind[1]]) grad_channels = ['gx', 'gy', 'gz'] for i in range(1, len(grad_channels) + 1): if event_ind[2 + (i - 1)] > 0: grad, compressed = SimpleNamespace(), SimpleNamespace() type = self.grad_library.type[event_ind[2 + (i - 1)]] lib_data = self.grad_library.data[event_ind[2 + (i - 1)]] grad.type = 'trap' if type == 't' else 'grad' grad.channel = grad_channels[i - 1][1] if grad.type == 'grad': amplitude = lib_data[0] shape_id = lib_data[1] delay = lib_data[2] shape_data = self.shape_library.data[shape_id] compressed.num_samples = shape_data[0] compressed.data = shape_data[1:] g = decompress_shape(compressed) grad.waveform = amplitude * g grad.t = np.arange(g.size) * self.grad_raster_time grad.delay = delay if len(lib_data) > 4: grad.first = lib_data[3] grad.last = lib_data[4] else: grad.first = grad.waveform[0] grad.last = grad.waveform[-1] else: if max(lib_data.shape) < 5: # added by GT grad.amplitude, grad.rise_time, grad.flat_time, grad.fall_time = [lib_data[x] for x in range(4)] grad.delay = 0 else: grad.amplitude, grad.rise_time, grad.flat_time, grad.fall_time, grad.delay = [lib_data[x] for x in range(5)] grad.area = grad.amplitude * (grad.flat_time + grad.rise_time / 2 + grad.fall_time / 2) grad.flat_area = grad.amplitude * grad.flat_time setattr(block, grad_channels[i - 1], grad) if event_ind[5] > 0: lib_data = self.adc_library.data[event_ind[5]] if max(lib_data.shape) < 6: lib_data = np.append(lib_data, 0) 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 = int(adc.num_samples) adc.type = 'adc' block.adc = adc + + # inserted for trigger support by mveldmann + if event_ind[6] > 0: + trig = SimpleNamespace() + trig.type = 'trigger' + trig.delay = self.extension_library.data[event_ind[6]][0] + trig.duration = self.extension_library.data[event_ind[6]][1] + block.trig = trig + return block diff --git a/pypulseq/Sequence/read_seq.py b/pypulseq/Sequence/read_seq.py index 83c928a..a070009 100755 --- a/pypulseq/Sequence/read_seq.py +++ b/pypulseq/Sequence/read_seq.py @@ -1,279 +1,323 @@ from pathlib import Path import numpy as np from pypulseq.event_lib import EventLibrary def read(self, path: str, **kwargs): """ Reads a `.seq` file from `path`. Parameters ---------- path : Path Path of .seq file to be read. """ detect_rf_use = True if 'detect_rf_use' in kwargs else False input_file = open(path, 'r') self.shape_library = EventLibrary() self.rf_library = EventLibrary() self.grad_library = EventLibrary() self.adc_library = EventLibrary() self.delay_library = EventLibrary() self.block_events = {} self.rf_raster_time = self.system.rf_raster_time self.grad_raster_time = self.system.grad_raster_time self.definitions = {} jemris_generated = False while True: section = __skip_comments(input_file) if section == -1: break if section == '[DEFINITIONS]': self.definitions = __read_definitions(input_file) elif section == '[JEMRIS]': # Added by GT jemris_generated = True elif section == '[VERSION]': version_major, version_minor, version_revision = __read_version(input_file) if version_major != self.version_major or version_minor != self.version_minor or version_revision != self.version_revision: raise Exception( f'Unsupported version: {version_major, version_minor, version_revision}. ' f'Expected: {self.version_major, self.version_minor, self.version_revision}') elif section == '[BLOCKS]': self.block_events = __read_blocks(input_file) elif section == '[RF]': if jemris_generated: self.rf_library = __read_events(input_file, [1, 1, 1, 1, 1]) else: self.rf_library = __read_events(input_file, [1, 1, 1, 1e-6, 1, 1]) elif section == '[GRADIENTS]': self.grad_library = __read_events(input_file, [1, 1, 1e-6], 'g', self.grad_library) elif section == '[TRAP]': if jemris_generated: self.grad_library = __read_events(input_file, [1, 1e-6, 1e-6, 1e-6], 't', self.grad_library) 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]) elif section == '[DELAYS]': self.delay_library = __read_events(input_file, 1e-6) elif section == '[SHAPES]': self.shape_library = __read_shapes(input_file) + # inserted for trigger support by mveldmann + elif section == '[EXTENSIONS]': + line = __strip_line(input_file) + elif section == 'extension TRIGGERS 1': + self.extension_library = __read_triggers(input_file, [1e-6, 1e-6], event_library=self.extension_library) + ### else: raise ValueError(f'Unknown section code: {section}') if detect_rf_use: for k in self.rf_library.keys(): lib_data = self.rf_library.data[k] rf = self.rf_from_lib_data(lib_data) flip_deg = abs(np.sum(rf.signal)) * rf.t[0] * 360 if len(lib_data) < 9: if flip_deg <= 90: lib_data[8] = 0 else: lib_data[8] = 2 self.rf_library.data[k] = lib_data def __read_definitions(input_file) -> dict: """ Read definitions from .seq file. Parameters ---------- input_file : file object .seq file Returns ------- definitions : dict Dict object containing key value pairs of definitions. """ definitions = dict() line = __strip_line(input_file) while line != '' and line[0] != '#': tok = line.split(' ') try: # Try converting every element into a float [float(x) for x in tok[1:]] definitions[tok[0]] = np.array(tok[1:], dtype=float) except ValueError: # Try clause did not work! definitions[tok[0]] = tok[1:] line = __strip_line(input_file) return definitions def __read_version(input_file) -> tuple: """ Read version from .seq file. Parameters ---------- input_file : file object .seq file Returns ------- tuple Tuple of major, minor and revision number. """ line = __strip_line(input_file) while line != '' and line[0] != '#': tok = line.split(' ') if tok[0] == 'major': major = np.array(tok[1:], dtype=float) elif tok[0] == 'minor': minor = np.array(tok[1:], dtype=float) elif tok[0] == 'revision': revision = np.array(tok[1:], dtype=float) line = __strip_line(input_file) return major, minor, revision def __read_blocks(input_file) -> dict: """ Read Pulseq blocks from .seq file. Parameters ---------- input_file : file .seq file Returns ------- event_table : dict Dict object containing key value pairs of Pulseq block ID and block definition. """ line = __strip_line(input_file) event_table = dict() while line != '' and line != '#': block_events = np.fromstring(line, dtype=int, sep=' ') event_table[block_events[0]] = block_events[1:] line = __strip_line(input_file) return event_table def __read_events(input_file, scale, type: str = None, event_library: EventLibrary = None) -> EventLibrary: """ Read Pulseq events from .seq file. Parameters ---------- input_file : file object .seq file scale : int Scaling factor. type : str Type of Pulseq event. event_library : EventLibrary EventLibrary Returns ------- definitions : dict `EventLibrary` object containing Pulseq event definitions. """ scale = 1 if scale is None else scale event_library = event_library if event_library is not None else EventLibrary() line = __strip_line(input_file) while line != '' and line != '#': data = np.fromstring(line, dtype=float, sep=' ') id = data[0] data = np.round(data[1:] * scale, decimals=6) if type is None: event_library.insert(key_id=id, new_data=data) else: event_library.insert(key_id=id, new_data=data, data_type=type) line = __strip_line(input_file) return event_library +### inserted for trigger support by mveldmann +def __read_triggers(input_file, scale, type: str = None, event_library: EventLibrary = None) -> EventLibrary: + """ + Read Pulseq triggers from .seq file. + + Parameters + ---------- + input_file : file object + .seq file + scale : int + Scaling factor. + type : str + Type of Pulseq event. + event_library : EventLibrary + EventLibrary + + Returns + ------- + definitions : dict + `EventLibrary` object containing Pulseq trigger definitions. + """ + scale = 1 if scale is None else scale + event_library = event_library if event_library is not None else EventLibrary() + + line = __strip_line(input_file) + + while line != '' and line != '#': + data = np.fromstring(line, dtype=float, sep=' ') + id = data[0] + data = np.round(data[3:] * scale, decimals=6) + if type is None: + event_library.insert(key_id=id, new_data=data) + else: + event_library.insert(key_id=id, new_data=data, data_type=type) + line = __strip_line(input_file) + + return event_library +### def __read_shapes(input_file) -> EventLibrary: """ Read Pulseq shapes from .seq file. Parameters ---------- input_file : file .seq file Returns ------- shape_library : EventLibrary `EventLibrary` object containing shape definitions. """ shape_library = EventLibrary() line = __skip_comments(input_file) while line != -1 and (line != '' or line[0:8] == 'shape_id'): tok = line.split(' ') id = int(tok[1]) line = __skip_comments(input_file) tok = line.split(' ') num_samples = int(tok[1]) data = [] line = __skip_comments(input_file) while line != '' and line != '#': data.append(float(line)) line = __strip_line(input_file) line = __skip_comments(input_file) data.insert(0, num_samples) data = np.asarray(data) shape_library.insert(key_id=id, new_data=data) return shape_library def __skip_comments(input_file) -> str: """ Skip one '#' comment in .seq file. Parameters ---------- input_file : file .seq file Returns ------- line : str First line in `input_file` after skipping one '#' comment block. Note: File pointer is remembered, so successive calls work as expected. """ line = __strip_line(input_file) while line != -1 and (line == '' or line[0] == '#'): line = __strip_line(input_file) return line def __strip_line(input_file): """ Removes spaces and newline whitespaces. Parameters ---------- input_file : file .seq file Returns ------- line : str First line in input_file after removing spaces and newline whitespaces. Note: File pointer is remembered, so successive calls work as expected. Returns -1 for eof. """ line = input_file.readline() # If line is an empty string, end of the file has been reached return line.strip() if line != '' else -1 diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py index 9572926..a4b9a7b 100755 --- a/pypulseq/Sequence/sequence.py +++ b/pypulseq/Sequence/sequence.py @@ -1,480 +1,482 @@ import math from types import SimpleNamespace import matplotlib as mpl import numpy as np mpl.use('TkAgg') from matplotlib import pyplot as plt from pypulseq.Sequence.test_report import test_report as ext_test_report from pypulseq.check_timing import check_timing as ext_check_timing from pypulseq.Sequence import block from pypulseq.Sequence.read_seq import read from pypulseq.Sequence.write_seq import write from pypulseq.calc_duration import calc_duration from pypulseq.calc_rf_center import calc_rf_center from pypulseq.decompress_shape import decompress_shape from pypulseq.event_lib import EventLibrary from pypulseq.opts import Opts from pypulseq.points_to_waveform import points_to_waveform class Sequence: def __init__(self, system=Opts()): self.version_major = 1 self.version_minor = 2 self.version_revision = 0 self.system = system self.definitions = dict() self.grad_library = EventLibrary() self.shape_library = EventLibrary() self.rf_library = EventLibrary() self.adc_library = EventLibrary() self.delay_library = EventLibrary() + self.extension_library = EventLibrary() # inserted for trigger support by mveldmann self.block_events = {} self.rf_raster_time = self.system.rf_raster_time self.grad_raster_time = self.system.grad_raster_time def __str__(self): s = "Sequence:" s += "\nshape_library: " + str(self.shape_library) s += "\nrf_library: " + str(self.rf_library) s += "\ngrad_library: " + str(self.grad_library) s += "\nadc_library: " + str(self.adc_library) s += "\ndelay_library: " + str(self.delay_library) + s += "\nextension_library: " + str(self.extension_library) # inserted for trigger support by mveldmann s += "\nrf_raster_time: " + str(self.rf_raster_time) s += "\ngrad_raster_time: " + str(self.grad_raster_time) s += "\nblock_events: " + str(len(self.block_events)) return s def duration(self): """ Get duration of this sequence. Returns ------- duration : float Duration of this sequence in millis. num_blocks : int Number of blocks in this sequence. event_count : int Number of events in this sequence. """ num_blocks = len(self.block_events) event_count = np.zeros(len(self.block_events[1])) duration = 0 for i in range(num_blocks): block = self.get_block(i + 1) event_count += self.block_events[i + 1] > 0 duration += calc_duration(block) return duration, num_blocks, event_count def check_timing(self): """ Check timing of the events in each block based on grad raster time system limit. Returns ------- is_ok : bool Boolean flag indicating timing errors. error_report : str Error report in case of timing errors. """ num_blocks = len(self.block_events) is_ok = True error_report = [] for i in range(num_blocks): block = self.get_block(i + 1) event_names = ['rf', 'gx', 'gy', 'gz', 'adc', 'delay'] ind = [hasattr(block, 'rf'), hasattr(block, 'gx'), hasattr(block, 'gy'), hasattr(block, 'gz'), hasattr(block, 'adc'), hasattr(block, 'delay')] ev = [getattr(block, event_names[i]) for i in range(len(event_names)) if ind[i] == 1] res, rep = ext_check_timing(self, *ev) is_ok = is_ok and res if len(rep) != 0: error_report.append(f'Block: {i} - {rep}\n') return is_ok, error_report def test_report(self): """ Analyze the sequence and return a text report. """ return ext_test_report(self) def set_definition(self, key: str, val: str): """ Sets custom definition to the `Sequence`. Parameters ---------- key : str Definition key. val : str Definition value. """ self.definitions[key] = val def get_definition(self, key: str) -> str: """ Retrieves definition identified by `key` from `Sequence`. Returns `None` if no matching definition is found. Parameters ---------- key : str Key of definition to retrieve. Returns ------- str Definition identified by `key` if found, else returns `None`. """ if key in self.definitions: return self.definitions[key] else: return None def add_block(self, *args): """ Adds event(s) as a block to `Sequence`. Parameters ---------- args Event or list of events to be added as a block to `Sequence`. """ block.add_block(self, len(self.block_events) + 1, *args) def get_block(self, block_index: int) -> SimpleNamespace: """ Retrieves block of events identified by `block_index` from `Sequence`. Parameters ---------- block_index : int Index of block to be retrieved from `Sequence`. Returns ------- SimpleNamespace Block identified by `block_index`. """ return block.get_block(self, block_index) def read(self, file_path: str): """ Read `.seq` file from `file_path`. Parameters ---------- file_path : str Path to `.seq` file to be read. """ read(self, file_path) def write(self, name: str): """ Writes the calling `Sequence` object as a `.seq` file with filename `name`. Parameters ---------- name :str Filename of `.seq` file to be written to disk. """ write(self, name) def rf_from_lib_data(self, lib_data): """ Construct RF object from `lib_data`. Parameters ---------- lib_data : list RF envelope. Returns ------- rf : SimpleNamespace RF object constructed from lib_data. """ rf = SimpleNamespace() rf.type = 'rf' amplitude, mag_shape, phase_shape = lib_data[0], lib_data[1], lib_data[2] shape_data = self.shape_library.data[mag_shape] compressed = SimpleNamespace() compressed.num_samples = shape_data[0] compressed.data = shape_data[1:] mag = decompress_shape(compressed) shape_data = self.shape_library.data[phase_shape] compressed.num_samples = shape_data[0] compressed.data = shape_data[1:] phase = decompress_shape(compressed) rf.signal = 1j * 2 * np.pi * phase rf.signal = amplitude * mag * np.exp(rf.signal) rf.t = np.arange(1, max(mag.shape) + 1) * self.rf_raster_time if max(lib_data.shape) < 6: rf.delay = 0 rf.freq_offset = lib_data[3] rf.phase_offset = lib_data[4] lib_data = np.append(lib_data, 0) else: rf.delay = lib_data[3] rf.freq_offset = lib_data[4] rf.phase_offset = lib_data[5] if max(lib_data.shape) < 7: lib_data = np.append(lib_data, 0) rf.dead_time = lib_data[6] if max(lib_data.shape) < 8: lib_data = np.append(lib_data, 0) rf.ringdown_time = lib_data[7] if max(lib_data.shape) < 9: lib_data = np.append(lib_data, 0) use_cases = {1: 'excitation', 2: 'refocusing', 3: 'inversion'} if lib_data[8] in use_cases: rf.use = use_cases[lib_data[8]] return rf def calculate_kspace(self, trajectory_delay: int = 0): """ Calculates the k-space trajectory of the entire pulse sequence. Parameters ---------- trajectory_delay : int Compensation factor in millis to align ADC and gradients in the reconstruction. Returns ------- k_traj_adc : numpy.ndarray K-space trajectory sampled at `t_adc` timepoints. k_traj : numpy.ndarray K-space trajectory of the entire pulse sequence. t_excitation : numpy.ndarray Excitation timepoints. t_refocusing : numpy.ndarray Refocusing timepoints. t_adc : numpy.ndarray Sampling timepoints. """ c_excitation = 0 c_refocusing = 0 c_adc_samples = 0 for i in range(len(self.block_events)): block = self.get_block(i + 1) if hasattr(block, 'rf'): if not hasattr(block.rf, 'use') or block.rf.use != 'refocusing': c_excitation += 1 else: c_refocusing += 1 if hasattr(block, 'adc'): c_adc_samples += int(block.adc.num_samples) t_excitation = np.zeros(c_excitation) t_refocusing = np.zeros(c_refocusing) k_time = np.zeros(c_adc_samples) current_dur = 0 c_excitation = 0 c_refocusing = 0 k_counter = 0 traj_recon_delay = trajectory_delay for i in range(len(self.block_events)): block = self.get_block(i + 1) if hasattr(block, 'rf'): rf = block.rf rf_center, _ = calc_rf_center(rf) t = rf.delay + rf_center if not hasattr(block.rf, 'use') or block.rf.use != 'refocusing': t_excitation[c_excitation] = current_dur + t c_excitation += 1 else: t_refocusing[c_refocusing] = current_dur + t c_refocusing += 1 if hasattr(block, 'adc'): k_time[k_counter:k_counter + block.adc.num_samples] = np.arange( block.adc.num_samples) * block.adc.dwell + block.adc.delay + current_dur + traj_recon_delay k_counter += block.adc.num_samples current_dur += calc_duration(block) gw = self.gradient_waveforms() i_excitation = np.round(t_excitation / self.grad_raster_time) i_refocusing = np.round(t_refocusing / self.grad_raster_time) k_traj = np.zeros(gw.shape) k = [0, 0, 0] for i in range(gw.shape[1]): k += gw[:, i] * self.grad_raster_time k_traj[:, i] = k if len(np.where(i_excitation == i + 1)[0]) >= 1: k = 0 k_traj[:, i] = np.nan if len(np.where(i_refocusing == i + 1)[0]) >= 1: k = -k k_traj_adc = [] for i in range(k_traj.shape[0]): k_traj_adc.append(np.interp(k_time, np.arange(1, k_traj.shape[1] + 1) * self.grad_raster_time, k_traj[i])) k_traj_adc = np.asarray(k_traj_adc) t_adc = k_time return k_traj_adc, k_traj, t_excitation, t_refocusing, t_adc def gradient_waveforms(self) -> np.ndarray: """ Decompress the entire gradient waveform. Returns an array of shape `gradient_axesxtimepoints`. `gradient_axes` is typically 3. Returns ------- grad_waveforms : numpy.ndarray Decompressed gradient waveform. """ duration, num_blocks, _ = self.duration() wave_length = math.ceil(duration / self.grad_raster_time) grad_channels = 3 grad_waveforms = np.zeros((grad_channels, wave_length)) grad_channels = ['gx', 'gy', 'gz'] t0 = 0 t0_n = 0 for i in range(num_blocks): block = self.get_block(i + 1) for j in range(len(grad_channels)): if hasattr(block, grad_channels[j]): grad = getattr(block, grad_channels[j]) if grad.type == 'grad': nt_start = round((grad.delay + grad.t[0]) / self.grad_raster_time) waveform = grad.waveform else: nt_start = round(grad.delay / self.grad_raster_time) if abs(grad.flat_time) > np.finfo(float).eps: t = np.cumsum([0, grad.rise_time, grad.flat_time, grad.fall_time]) trap_form = np.multiply([0, 1, 1, 0], grad.amplitude) else: t = np.cumsum([0, grad.rise_time, grad.fall_time]) trap_form = np.multiply([0, 1, 0], grad.amplitude) tn = math.floor(t[-1] / self.grad_raster_time) t = np.append(t, t[-1] + self.grad_raster_time) trap_form = np.append(trap_form, 0) if abs(grad.amplitude) > np.finfo(float).eps: waveform = points_to_waveform(t, trap_form, self.grad_raster_time) else: waveform = np.zeros(tn + 1) if waveform.size != np.sum(np.isfinite(waveform)): raise Warning('Not all elements of the generated waveform are finite') """ Matlab dynamically resizes arrays during slice assignment operation if assignment is out of bounds Numpy does not Following lines are a workaround """ l1, l2 = int(t0_n + nt_start), int(t0_n + nt_start + max(waveform.shape)) if l2 > grad_waveforms.shape[1]: grad_waveforms.resize((len(grad_channels), l2)) grad_waveforms[j, l1:l2] = waveform t0 += calc_duration(block) t0_n = round(t0 / self.grad_raster_time) return grad_waveforms def plot(self, type: str = 'Gradient', time_range=(0, np.inf), time_disp: str = 's', save: bool=False): """ Plot `Sequence`. Parameters ---------- type : str Gradients display type, must be one of either 'Gradient' or 'Kspace'. time_range : List Time range (x-axis limits) for plotting the sequence. Default is 0 to infinity (entire sequence). time_disp : str Time display type, must be one of `s`, `ms` or `us`. save_as : bool Boolean flag indicating if plots should be saved. The two figures will be saved as JPG with numerical suffixes to the filename 'seq_plot'. """ valid_plot_types = ['Gradient', 'Kspace'] valid_time_units = ['s', 'ms', 'us'] if type not in valid_plot_types: raise Exception() if time_disp not in valid_time_units: raise Exception() fig1, fig2 = plt.figure(1), plt.figure(2) sp11 = fig1.add_subplot(311) sp12, sp13 = fig1.add_subplot(312, sharex=sp11), fig1.add_subplot(313, sharex=sp11) fig2_sp_list = [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)] t0 = 0 for iB in range(1, len(self.block_events) + 1): block = self.get_block(iB) is_valid = time_range[0] <= t0 <= time_range[1] if is_valid: if hasattr(block, 'adc'): adc = block.adc t = adc.delay + [(x * adc.dwell) for x in range(0, int(adc.num_samples))] sp11.plot((t0 + t), np.zeros(len(t)), 'rx') if hasattr(block, 'rf'): rf = block.rf tc, ic = calc_rf_center(rf) t = rf.t + rf.delay tc = tc + rf.delay sp12.plot(t_factor * (t0 + t), abs(rf.signal)) sp13.plot(t_factor * (t0 + t), np.angle( rf.signal * np.exp(1j * rf.phase_offset) * np.exp(1j * 2 * math.pi * rf.t * rf.freq_offset)), t_factor * (t0 + tc), np.angle(rf.signal[ic] * np.exp(1j * rf.phase_offset) * np.exp( 1j * 2 * math.pi * rf.t[ic] * rf.freq_offset)), 'xb') grad_channels = ['gx', 'gy', 'gz'] for x in range(0, len(grad_channels)): if hasattr(block, grad_channels[x]): grad = getattr(block, grad_channels[x]) if grad.type == 'grad': # In place unpacking of grad.t with the starred expression t = grad.delay + [0, *(grad.t + (grad.t[1] - grad.t[0]) / 2), grad.t[-1] + grad.t[1] - grad.t[0]] waveform = np.array([grad.first, grad.last]) waveform = 1e-3 * np.insert(waveform, 1, grad.waveform) else: t = np.cumsum([0, grad.delay, grad.rise_time, grad.flat_time, grad.fall_time]) waveform = 1e-3 * grad.amplitude * np.array([0, 0, 1, 1, 0]) fig2_sp_list[x].plot(t_factor * (t0 + t), waveform) t0 += calc_duration(block) grad_plot_labels = ['x', 'y', 'z'] sp11.set_ylabel('ADC') sp12.set_ylabel('RF mag (Hz)') sp13.set_ylabel('RF phase (rad)') [fig2_sp_list[x].set_ylabel(f'G{grad_plot_labels[x]} (kHz/m)') for x in range(3)] # Setting display limits disp_range = t_factor * np.array([time_range[0], min(t0, time_range[1])]) sp11.set_xlim(disp_range) sp12.set_xlim(disp_range) sp13.set_xlim(disp_range) [x.set_xlim(disp_range) for x in fig2_sp_list] fig1.tight_layout() fig2.tight_layout() if save: fig1.savefig('seq_plot1.jpg') fig2.savefig('seq_plot2.jpg') plt.show() diff --git a/pypulseq/Sequence/write_seq.py b/pypulseq/Sequence/write_seq.py index 89f747f..e02b59d 100755 --- a/pypulseq/Sequence/write_seq.py +++ b/pypulseq/Sequence/write_seq.py @@ -1,135 +1,154 @@ import numpy as np def write(self, file_name): """ Writes the calling `Sequence` object as a `.seq` file with filename `file_name`. Parameters ---------- file_name : str File name of `.seq` file to be written to disk. """ # `>.0f` is used when only decimals have to be displayed. # `>g` is used when insignificant zeros have to be truncated. file_name += '.seq' if file_name[-4:] != '.seq' not in file_name else '' output_file = open(file_name, 'w') output_file.write('# Pulseq sequence file\n') output_file.write('# Created by PyPulseq\n\n') output_file.write('[VERSION]\n') output_file.write(f'major {self.version_major}\n') output_file.write(f'minor {self.version_minor}\n') output_file.write(f'revision {self.version_revision}\n') output_file.write('\n') if len(self.definitions) != 0: output_file.write('[DEFINITIONS]\n') keys = list(self.definitions.keys()) values = list(self.definitions.values()) for i in range(len(keys)): output_file.write(f'{keys[i]} ') if isinstance(values[i], str): output_file.write(f'{values[i]} ') else: for t in range(len(values[i])): output_file.write(f'{values[i][t]:.6g} ') output_file.write('\n') output_file.write('\n') output_file.write('# Format of blocks:\n') - output_file.write('# # D RF GX GY GZ ADC\n') + output_file.write('# # D RF GX GY GZ ADC EXT\n') # EXT not part of pypulseq - inserted for trigger support by mveldmann output_file.write('[BLOCKS]\n') id_format_width = '{:' + str(len(str(len(self.block_events)))) + 'd} ' - id_format_str = id_format_width + '{:2d} {:2d} {:3d} {:3d} {:3d} {:2d}\n' + id_format_str = id_format_width + '{:2d} {:2d} {:3d} {:3d} {:3d} {:2d} {:2d}\n' # EXT not part of pypulseq - inserted for trigger support by mveldmann for i in range(len(self.block_events)): s = id_format_str.format(*np.insert(self.block_events[i + 1], 0, (i + 1))) output_file.write(s) output_file.write('\n') if len(self.rf_library.keys) != 0: output_file.write('# Format of RF events:\n') output_file.write('# id amplitude mag_id phase_id delay freq phase\n') output_file.write('# .. Hz .... .... us Hz rad\n') output_file.write('[RF]\n') rf_lib_keys = self.rf_library.keys # id_format_str = '{:>1.0f} {:>12g} {:>1.0f} {:>1.0f} {:>g} {:>g}\n' id_format_str = '{:.0f} {:12g} {:.0f} {:.0f} {:g} {:g} {:g}\n' for k in rf_lib_keys.keys(): lib_data1 = self.rf_library.data[k][0:3] lib_data2 = self.rf_library.data[k][4:7] delay = round(self.rf_library.data[k][3] * 1e6) s = id_format_str.format(k, *lib_data1, delay, *lib_data2) output_file.write(s) output_file.write('\n') grad_lib_values = np.array(list(self.grad_library.type.values())) arb_grad_mask = np.where(grad_lib_values == 'g')[0] trap_grad_mask = np.where(grad_lib_values == 't')[0] if len(arb_grad_mask) != 0: output_file.write('# Format of arbitrary gradients:\n') output_file.write('# id amplitude shape_id delay\n') output_file.write('# .. Hz/m .... us\n') output_file.write('[GRADIENTS]\n') id_format_str = '{:.0f} {:12g} {:.0f} {:.0f}\n' keys = np.array(list(self.grad_library.keys.keys())) for k in keys[arb_grad_mask]: s = id_format_str.format(k, *self.grad_library.data[k][:2], round(self.grad_library.data[k][2] * 1e6)) output_file.write(s) output_file.write('\n') if len(trap_grad_mask) != 0: output_file.write('# Format of trapezoid gradients:\n') output_file.write('# id amplitude rise flat fall delay\n') output_file.write('# .. Hz/m us us us us\n') output_file.write('[TRAP]\n') keys = np.array(list(self.grad_library.keys.keys())) id_format_str = '{:2g} {:12g} {:3g} {:4g} {:3g} {:3g}\n' for k in keys[trap_grad_mask]: data = np.copy(self.grad_library.data[k]) # Make a copy to leave the original untouched data[1:] = np.round(1e6 * data[1:]) s = id_format_str.format(k, *data) output_file.write(s) output_file.write('\n') if len(self.adc_library.keys) != 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('[ADC]\n') keys = self.adc_library.keys id_format_str = '{:.0f} {:.0f} {:.0f} {:.0f} {:.0f} {:.6g}\n' for k in keys.values(): data = np.multiply(self.adc_library.data[k][0:5], [1, 1e9, 1e6, 1, 1]) s = id_format_str.format(k, *data) output_file.write(s) output_file.write('\n') if len(self.delay_library.keys) != 0: output_file.write('# Format of delays:\n') output_file.write('# id delay (us)\n') output_file.write('[DELAYS]\n') keys = self.delay_library.keys id_format_str = '{:.0f} {:.0f}\n' for k in keys.values(): s = id_format_str.format(k, *np.round(1e6 * self.delay_library.data[k])) output_file.write(s) output_file.write('\n') + # inserted for trigger support by mveldmann + if len(self.extension_library.keys) != 0: + output_file.write('# Format of extensions:\n') + output_file.write('# id type ref next\n') + output_file.write('[EXTENSIONS]\n') + keys = self.extension_library.keys + id_format_str = '{:.0f} {:.0f} {:.0f} {:.0f}\n' + for k in keys.values(): + s = id_format_str.format(k, 1, k, 0) # always use extension type 1 - TRIGGERS + output_file.write(s) + output_file.write('\n') + output_file.write('extension TRIGGERS 1\n') + ext_format_str = '{:.0f} {:.0f} {:.0f} {:.0f} {:.0f}\n' + for k in keys.values(): + s = ext_format_str.format(k, 1, 3, *np.round(1e6 * self.extension_library.data[k])) # use always digital output trigger (trigger type:1, channel: 3) + output_file.write(s) + output_file.write('\n') + + if len(self.shape_library.keys) != 0: output_file.write('# Sequence Shapes\n') output_file.write('[SHAPES]\n\n') keys = self.shape_library.keys for k in keys.values(): shape_data = self.shape_library.data[k] s = 'shape_id {:.0f}\n' s = s.format(k) output_file.write(s) s = 'num_samples {:.0f}\n' s = s.format(shape_data[0]) output_file.write(s) s = '{:.9g}\n' * len(shape_data[1:]) s = s.format(*shape_data[1:]) output_file.write(s) output_file.write('\n') diff --git a/pypulseq/make_trig.py b/pypulseq/make_trig.py new file mode 100644 index 0000000..4ebddd6 --- /dev/null +++ b/pypulseq/make_trig.py @@ -0,0 +1,33 @@ +# inserted for trigger support by mveldmann + +from types import SimpleNamespace +from pypulseq.opts import Opts + +def make_trig(delay: float = 0, duration: float = 1e-5 , system: Opts = Opts()) -> SimpleNamespace: + """ + Creates a trigger event. + + Parameters + ---------- + delay : float + Delay in seconds + duration: float + Duration in seconds + + Returns + ------- + trigger : SimpleNamespace + Trigger event. + """ + + trig = SimpleNamespace() + if delay < 0: + raise ValueError('Delay {:.2f} ms is invalid'.format(delay * 1e3)) + trig.type = 'trigger' + trig.delay = delay + if duration < system.grad_raster_time: + print('Duration too short, set to gradient raster time (minimum duration)') + trig.duration = system.grad_raster_time + else: + trig.duration = duration + return trig