From ee9210cf6140e1d677e7e45489153ca044d6a1e3 Mon Sep 17 00:00:00 2001 From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com> Date: Thu, 17 Aug 2023 11:57:31 +0200 Subject: [PATCH 1/2] Bug fixes in Sequence and make_trapezoid Sequence.waveforms_and_times: - Fixed missing implementation for recompressing trapezoids from arbitrary gradients. Now proceeds without recompressing and added a TODO to implement the MATLAB code. - Fixed issue where rounding errors could create non-monotonically increasing times. Also moved the non-monotonicity check outside the loop. make_trapezoid: - Fixed bug that would give an 'requested area is too large' error when explicitly specifiying a triangular gradient with duration, rise_time, and area. - Fixed check for area == 0, while an unspecified area would be None --- pypulseq/Sequence/sequence.py | 33 ++++++++++++++++++--------------- pypulseq/make_trapezoid.py | 4 ++-- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py index 28aaebd..c35a3cd 100755 --- a/pypulseq/Sequence/sequence.py +++ b/pypulseq/Sequence/sequence.py @@ -1391,13 +1391,16 @@ def waveforms_and_times( but first we have to restore samples on the edges of the gradient raster intervals for that we need the first sample. """ - max_abs = np.max(np.abs(grad.waveform)) - odd_step1 = np.array([grad.first, *(2 * grad.waveform)]) - odd_step2 = odd_step1 * ( - np.mod(np.arange(1, len(odd_step1) + 1), 2) * 2 - 1 - ) - waveform_odd_rest = np.cumsum(odd_step2) * ( - np.mod(np.arange(1, len(odd_step2)), 2) * 2 - 1 + + # TODO: Implement restoreAdditionalShapeSamples + # https://github.com/pulseq/pulseq/blob/master/matlab/%2Bmr/restoreAdditionalShapeSamples.m + + out_len[j] += len(grad.tt) + shape_pieces[j, block_counter] = np.array( + [ + curr_dur + grad.delay + grad.tt, + grad.waveform, + ] ) else: # Extended trapezoid out_len[j] += len(grad.tt) @@ -1502,14 +1505,14 @@ def waveforms_and_times( wave_data[j] = np.zeros((2, int(out_len[j])), dtype=np.complex128) wave_cnt = np.zeros(shape_channels, dtype=int) - for block_counter in range(num_blocks): - for j in range(shape_channels): + for j in range(shape_channels): + for block_counter in range(num_blocks): if shape_pieces[j, block_counter] is not None: wave_data_local = shape_pieces[j, block_counter] length = wave_data_local.shape[1] if ( wave_cnt[j] == 0 - or wave_data[j][0, wave_cnt[j] - 1] != wave_data_local[0, 0] + or wave_data[j][0, wave_cnt[j] - 1] < wave_data_local[0, 0] ): wave_data[j][ :, wave_cnt[j] + np.arange(length) @@ -1521,11 +1524,11 @@ def waveforms_and_times( ] = wave_data_local[:, 1:] wave_cnt[j] += length - 1 - rftdiff = np.diff(wave_data[j][0, : wave_cnt[j]]) - if np.any(rftdiff < eps): - raise Warning( - "Time vector elements are not monotonically increasing." - ) + rftdiff = np.diff(wave_data[j][0]) + if np.any(rftdiff < eps): + raise Warning( + "Time vector elements are not monotonically increasing." + ) diff --git a/pypulseq/make_trapezoid.py b/pypulseq/make_trapezoid.py index b670a96..0112f06 100644 --- a/pypulseq/make_trapezoid.py +++ b/pypulseq/make_trapezoid.py @@ -130,7 +130,7 @@ def make_trapezoid( fall_time = rise_time amplitude2 = area / (duration - 0.5 * rise_time - 0.5 * fall_time) possible = ( - duration > (rise_time + fall_time) and np.abs(amplitude2) < max_grad + duration >= (rise_time + fall_time) and np.abs(amplitude2) <= max_grad ) assert possible, ( f"Requested area is too large for this gradient. Probably amplitude is violated " @@ -153,7 +153,7 @@ def make_trapezoid( # Adjust amplitude (after rounding) to match area amplitude2 = area / (rise_time / 2 + fall_time / 2 + flat_time) else: - if area == 0: + if area == None: raise ValueError("Must supply area or duration.") else: # Find the shortest possible duration. First check if the area can be realized as a triangle. From 4aa12a1cfa13847b75b4f21b295ff6a5dd6782fa Mon Sep 17 00:00:00 2001 From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com> Date: Thu, 17 Aug 2023 17:11:14 +0200 Subject: [PATCH 2/2] More bug fixes - Actually fixed the non-monotonic time problem - Two small fixes in sequence plot() that cut off the start of the sequence when time_range is specified, and coloured labels wrongly. --- pypulseq/Sequence/sequence.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/pypulseq/Sequence/sequence.py b/pypulseq/Sequence/sequence.py index c35a3cd..c48cbe7 100755 --- a/pypulseq/Sequence/sequence.py +++ b/pypulseq/Sequence/sequence.py @@ -956,7 +956,8 @@ def plot( 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))) - label_colors_to_plot = np.roll(label_colors_to_plot, -1, axis=0).tolist() + cycler = mpl.cycler(color=label_colors_to_plot) + sp11.set_prop_cycle(cycler) # Block timings block_edges = np.cumsum([0, *self.block_durations]) @@ -970,7 +971,8 @@ def plot( for block_counter in range(len(self.block_events)): block = self.get_block(block_counter + 1) - is_valid = time_range[0] <= t0 <= time_range[1] + 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)): @@ -997,11 +999,6 @@ def plot( ) if label_defined and len(label_idx_to_plot) != 0: - cycler = mpl.cycler(color=label_colors_to_plot) - sp11.set_prop_cycle(cycler) - label_colors_to_plot = np.roll( - label_colors_to_plot, -1, axis=0 - ).tolist() 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 @@ -1512,7 +1509,7 @@ def waveforms_and_times( length = wave_data_local.shape[1] if ( wave_cnt[j] == 0 - or wave_data[j][0, wave_cnt[j] - 1] < wave_data_local[0, 0] + or wave_data[j][0, wave_cnt[j] - 1]+eps < wave_data_local[0, 0] ): wave_data[j][ :, wave_cnt[j] + np.arange(length) @@ -1524,7 +1521,7 @@ def waveforms_and_times( ] = wave_data_local[:, 1:] wave_cnt[j] += length - 1 - rftdiff = np.diff(wave_data[j][0]) + rftdiff = np.diff(wave_data[j][0,:wave_cnt[j]]) if np.any(rftdiff < eps): raise Warning( "Time vector elements are not monotonically increasing."