Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bug fixes in Sequence and make_trapezoid #128

Merged
merged 2 commits into from
Aug 18, 2023
Merged

Bug fixes in Sequence and make_trapezoid #128

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ee9210c
Bug fixes in Sequence and make_trapezoid
FrankZijlstra Aug 17, 2023
4aa12a1
More bug fixes
FrankZijlstra Aug 17, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 22 additions & 22 deletions pypulseq/Sequence/sequence.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand All @@ -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)):
Expand All @@ -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
Expand Down Expand Up @@ -1391,13 +1388,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)
Expand Down Expand Up @@ -1502,14 +1502,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]+eps < wave_data_local[0, 0]
):
wave_data[j][
:, wave_cnt[j] + np.arange(length)
Expand All @@ -1521,11 +1521,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,:wave_cnt[j]])
if np.any(rftdiff < eps):
raise Warning(
"Time vector elements are not monotonically increasing."
)



Expand Down
4 changes: 2 additions & 2 deletions pypulseq/make_trapezoid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 "
Expand All @@ -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.
Expand Down