From 50bcfec8319f340999ec418a0b00395f26f8be6f Mon Sep 17 00:00:00 2001 From: Mojtaba Date: Mon, 9 Mar 2026 09:59:09 +0100 Subject: [PATCH 1/5] Add rotate3D implementation --- src/pypulseq/rotate3D.py | 232 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100644 src/pypulseq/rotate3D.py diff --git a/src/pypulseq/rotate3D.py b/src/pypulseq/rotate3D.py new file mode 100644 index 00000000..b65f5c64 --- /dev/null +++ b/src/pypulseq/rotate3D.py @@ -0,0 +1,232 @@ +from types import SimpleNamespace +from typing import List, Union, Optional, Sequence, Any +import copy + +import numpy as np + +from pypulseq.add_gradients import add_gradients +from pypulseq.opts import Opts +from pypulseq.scale_grad import scale_grad + + +def __get_grad_abs_mag(grad: SimpleNamespace) -> float: + """Magnitude used for thresholding output components (MATLAB getGradAbsMag).""" + if grad.type == "trap": + return float(abs(grad.amplitude)) + return float(np.max(np.abs(grad.waveform))) + + +def align_gradient_to_raster(g: SimpleNamespace, raster: float) -> SimpleNamespace: + """ + Align gradient timing parameters to the gradient raster and verify rounding error. + This is used to avoid add_gradients() complaining about non-raster time points. + """ + tol = 1e-9 # seconds + + if getattr(g, "type", None) == "grad": + original_tt = np.array(g.tt, dtype=np.float64) + rounded_tt = np.round(original_tt / raster) * raster + if np.any(np.abs(original_tt - rounded_tt) > tol): + raise ValueError( + f"'grad' tt values not aligned to raster (>{tol}s): {original_tt} -> {rounded_tt}" + ) + g.tt = rounded_tt + g.waveform = np.array(g.waveform, dtype=np.float64) + + elif getattr(g, "type", None) == "trap": + original_delay = float(getattr(g, "delay", 0.0)) + original_rise = float(getattr(g, "rise_time", 0.0)) + original_flat = float(getattr(g, "flat_time", 0.0)) + original_fall = float(getattr(g, "fall_time", 0.0)) + + rounded_delay = float(np.round(original_delay / raster) * raster) + rounded_rise = float(np.round(original_rise / raster) * raster) + rounded_flat = float(np.round(original_flat / raster) * raster) + rounded_fall = float(np.round(original_fall / raster) * raster) + + if any( + abs(x - y) > tol + for x, y in [ + (original_delay, rounded_delay), + (original_rise, rounded_rise), + (original_flat, rounded_flat), + (original_fall, rounded_fall), + ] + ): + raise ValueError( + f"'trap' timing values not aligned to raster (>{tol}s): " + f"delay={original_delay}→{rounded_delay}, " + f"rise={original_rise}→{rounded_rise}, " + f"flat={original_flat}→{rounded_flat}, " + f"fall={original_fall}→{rounded_fall}" + ) + + g.delay = rounded_delay + g.rise_time = rounded_rise + g.flat_time = rounded_flat + g.fall_time = rounded_fall + + # always align delay if present + original_delay = float(getattr(g, "delay", 0.0)) + rounded_delay = float(np.round(original_delay / raster) * raster) + if abs(original_delay - rounded_delay) > tol: + raise ValueError( + f"'delay' value not aligned to raster (>{tol}s): {original_delay} → {rounded_delay}" + ) + g.delay = rounded_delay + + return g + + +def _quat_wxyz_to_rotmat(q_wxyz: np.ndarray) -> np.ndarray: + """ + Convert quaternion [w,x,y,z] (scalar first) to a 3x3 rotation matrix. + MATLAB: mr.aux.quat.toRotMat(rotation) + """ + q = np.asarray(q_wxyz, dtype=float).reshape(-1) + if q.size != 4: + raise ValueError("Quaternion must have length 4 [w, x, y, z].") + + w, x, y, z = q + n = np.sqrt(w * w + x * x + y * y + z * z) + if n == 0: + raise ValueError("Quaternion norm must be non-zero.") + # normalize (MATLAB expects unit quaternion; we tolerate small numeric drift) + w, x, y, z = w / n, x / n, y / n, z / n + + return np.array( + [ + [1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w)], + [2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w)], + [2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)], + ], + dtype=float, + ) + + +def _parse_rotation_to_matrix( + *, + rotation: Optional[Any] = None, + rotation_matrix: Optional[np.ndarray] = None, +) -> np.ndarray: + """ + MATLAB parity: + - rotation can be 3x3 matrix OR quaternion length-4 [w,x,y,z] + Backward compatibility: + - rotation_matrix keyword (your current API) still works. + """ + if rotation is None and rotation_matrix is None: + raise ValueError("You must provide either 'rotation' (matrix or quaternion) or 'rotation_matrix' (3x3).") + + if rotation_matrix is not None: + rot = np.asarray(rotation_matrix, dtype=float) + else: + rot = np.asarray(rotation, dtype=float) + + if rot.shape == (3, 3): + return rot + + # quaternion case + if rot.size == 4: + return _quat_wxyz_to_rotmat(rot) + + raise ValueError("rotation must be either a 3x3 matrix or a quaternion [w, x, y, z] (scalar first).") + + +def rotate3D( + *args: SimpleNamespace, + rotation: Optional[Any] = None, + rotation_matrix: Optional[np.ndarray] = None, + system: Union[Opts, None] = None, +) -> List[SimpleNamespace]: + """ + Rotate gradient events (trap/grad) by a 3x3 rotation matrix (or quaternion like MATLAB). + Non-gradient events (rf/adc/delay/etc) are passed through unchanged. + + Accepts mixed inputs like: + ms_rotate3D(gx, gy, rf, adc, rotation_matrix=R, system=sys) + ms_rotate3D(gx, gy, rf, adc, rotation=R, system=sys) + ms_rotate3D(gx, gy, rf, adc, rotation=[w,x,y,z], system=sys) # quaternion + + Also accepts list inputs: + ms_rotate3D([gx, gy, gz], rf, adc, rotation_matrix=R, system=sys) + """ + # ---- flatten list/tuple inputs inside args ---- + flat_args: List[SimpleNamespace] = [] + for a in args: + if a is None: + continue + if isinstance(a, (list, tuple)): + for x in a: + if x is not None: + flat_args.append(x) + else: + flat_args.append(a) + args = tuple(flat_args) + + if system is None: + system = Opts.default + + rotMat = _parse_rotation_to_matrix(rotation=rotation, rotation_matrix=rotation_matrix) + + axes = ["x", "y", "z"] + events_to_rotate: dict[str, SimpleNamespace] = {} + bypass: List[SimpleNamespace] = [] + + # ---- split into gradients vs bypass ---- + for event in args: + if not hasattr(event, "type") or event.type not in ("grad", "trap"): + bypass.append(event) + continue + + if not hasattr(event, "channel") or event.channel not in axes: + raise ValueError(f"Invalid or missing gradient channel: {getattr(event, 'channel', None)}") + + if event.channel in events_to_rotate: + raise ValueError(f"More than one gradient for channel: {event.channel}") + + # deepcopy to avoid mutating caller objects (matching your current code) + events_to_rotate[event.channel] = copy.deepcopy(event) + + # nothing to rotate -> just return bypass + if len(events_to_rotate) == 0: + return list(bypass) + + # ---- thresholding like MATLAB ---- + max_mag = 0.0 + for ax in axes: + if ax in events_to_rotate: + max_mag = max(max_mag, __get_grad_abs_mag(events_to_rotate[ax])) + fthresh = 1e-6 + thresh = fthresh * max_mag + + rotated_grads: List[SimpleNamespace] = [] + for j in range(3): + grad_out_curr = None + for i in range(3): + ax_in = axes[i] + if ax_in not in events_to_rotate: + continue + + w = float(rotMat[j, i]) + if abs(w) < fthresh: + continue + + scaled = scale_grad(grad=events_to_rotate[ax_in], scale=w) + scaled.channel = axes[j] + + if grad_out_curr is None: + grad_out_curr = scaled + else: + # raster-align before summing (your current approach) + grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) + scaled = align_gradient_to_raster(scaled, system.grad_raster_time) + + # add_gradients expects a list + grad_out_curr = add_gradients([grad_out_curr, scaled], system=system) + + # only output non-zero amplitude gradients + if grad_out_curr is not None and __get_grad_abs_mag(grad_out_curr) >= thresh: + rotated_grads.append(grad_out_curr) + + return [*bypass, *rotated_grads] \ No newline at end of file From 63ed6ddb8b46be9c77150bdb7ab3e5c2346be880 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 9 Mar 2026 09:49:02 +0000 Subject: [PATCH 2/5] [pre-commit] auto fixes from pre-commit hooks --- src/pypulseq/rotate3D.py | 54 +++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 29 deletions(-) diff --git a/src/pypulseq/rotate3D.py b/src/pypulseq/rotate3D.py index b65f5c64..6de225ae 100644 --- a/src/pypulseq/rotate3D.py +++ b/src/pypulseq/rotate3D.py @@ -1,6 +1,6 @@ -from types import SimpleNamespace -from typing import List, Union, Optional, Sequence, Any import copy +from types import SimpleNamespace +from typing import Any, List, Optional, Union import numpy as np @@ -11,7 +11,7 @@ def __get_grad_abs_mag(grad: SimpleNamespace) -> float: """Magnitude used for thresholding output components (MATLAB getGradAbsMag).""" - if grad.type == "trap": + if grad.type == 'trap': return float(abs(grad.amplitude)) return float(np.max(np.abs(grad.waveform))) @@ -23,21 +23,19 @@ def align_gradient_to_raster(g: SimpleNamespace, raster: float) -> SimpleNamespa """ tol = 1e-9 # seconds - if getattr(g, "type", None) == "grad": + if getattr(g, 'type', None) == 'grad': original_tt = np.array(g.tt, dtype=np.float64) rounded_tt = np.round(original_tt / raster) * raster if np.any(np.abs(original_tt - rounded_tt) > tol): - raise ValueError( - f"'grad' tt values not aligned to raster (>{tol}s): {original_tt} -> {rounded_tt}" - ) + raise ValueError(f"'grad' tt values not aligned to raster (>{tol}s): {original_tt} -> {rounded_tt}") g.tt = rounded_tt g.waveform = np.array(g.waveform, dtype=np.float64) - elif getattr(g, "type", None) == "trap": - original_delay = float(getattr(g, "delay", 0.0)) - original_rise = float(getattr(g, "rise_time", 0.0)) - original_flat = float(getattr(g, "flat_time", 0.0)) - original_fall = float(getattr(g, "fall_time", 0.0)) + elif getattr(g, 'type', None) == 'trap': + original_delay = float(getattr(g, 'delay', 0.0)) + original_rise = float(getattr(g, 'rise_time', 0.0)) + original_flat = float(getattr(g, 'flat_time', 0.0)) + original_fall = float(getattr(g, 'fall_time', 0.0)) rounded_delay = float(np.round(original_delay / raster) * raster) rounded_rise = float(np.round(original_rise / raster) * raster) @@ -55,10 +53,10 @@ def align_gradient_to_raster(g: SimpleNamespace, raster: float) -> SimpleNamespa ): raise ValueError( f"'trap' timing values not aligned to raster (>{tol}s): " - f"delay={original_delay}→{rounded_delay}, " - f"rise={original_rise}→{rounded_rise}, " - f"flat={original_flat}→{rounded_flat}, " - f"fall={original_fall}→{rounded_fall}" + f'delay={original_delay}→{rounded_delay}, ' + f'rise={original_rise}→{rounded_rise}, ' + f'flat={original_flat}→{rounded_flat}, ' + f'fall={original_fall}→{rounded_fall}' ) g.delay = rounded_delay @@ -67,12 +65,10 @@ def align_gradient_to_raster(g: SimpleNamespace, raster: float) -> SimpleNamespa g.fall_time = rounded_fall # always align delay if present - original_delay = float(getattr(g, "delay", 0.0)) + original_delay = float(getattr(g, 'delay', 0.0)) rounded_delay = float(np.round(original_delay / raster) * raster) if abs(original_delay - rounded_delay) > tol: - raise ValueError( - f"'delay' value not aligned to raster (>{tol}s): {original_delay} → {rounded_delay}" - ) + raise ValueError(f"'delay' value not aligned to raster (>{tol}s): {original_delay} → {rounded_delay}") g.delay = rounded_delay return g @@ -85,12 +81,12 @@ def _quat_wxyz_to_rotmat(q_wxyz: np.ndarray) -> np.ndarray: """ q = np.asarray(q_wxyz, dtype=float).reshape(-1) if q.size != 4: - raise ValueError("Quaternion must have length 4 [w, x, y, z].") + raise ValueError('Quaternion must have length 4 [w, x, y, z].') w, x, y, z = q n = np.sqrt(w * w + x * x + y * y + z * z) if n == 0: - raise ValueError("Quaternion norm must be non-zero.") + raise ValueError('Quaternion norm must be non-zero.') # normalize (MATLAB expects unit quaternion; we tolerate small numeric drift) w, x, y, z = w / n, x / n, y / n, z / n @@ -130,7 +126,7 @@ def _parse_rotation_to_matrix( if rot.size == 4: return _quat_wxyz_to_rotmat(rot) - raise ValueError("rotation must be either a 3x3 matrix or a quaternion [w, x, y, z] (scalar first).") + raise ValueError('rotation must be either a 3x3 matrix or a quaternion [w, x, y, z] (scalar first).') def rotate3D( @@ -169,21 +165,21 @@ def rotate3D( rotMat = _parse_rotation_to_matrix(rotation=rotation, rotation_matrix=rotation_matrix) - axes = ["x", "y", "z"] + axes = ['x', 'y', 'z'] events_to_rotate: dict[str, SimpleNamespace] = {} bypass: List[SimpleNamespace] = [] # ---- split into gradients vs bypass ---- for event in args: - if not hasattr(event, "type") or event.type not in ("grad", "trap"): + if not hasattr(event, 'type') or event.type not in ('grad', 'trap'): bypass.append(event) continue - if not hasattr(event, "channel") or event.channel not in axes: - raise ValueError(f"Invalid or missing gradient channel: {getattr(event, 'channel', None)}") + if not hasattr(event, 'channel') or event.channel not in axes: + raise ValueError(f'Invalid or missing gradient channel: {getattr(event, "channel", None)}') if event.channel in events_to_rotate: - raise ValueError(f"More than one gradient for channel: {event.channel}") + raise ValueError(f'More than one gradient for channel: {event.channel}') # deepcopy to avoid mutating caller objects (matching your current code) events_to_rotate[event.channel] = copy.deepcopy(event) @@ -229,4 +225,4 @@ def rotate3D( if grad_out_curr is not None and __get_grad_abs_mag(grad_out_curr) >= thresh: rotated_grads.append(grad_out_curr) - return [*bypass, *rotated_grads] \ No newline at end of file + return [*bypass, *rotated_grads] From 295b6895eed0c4dfa146da1f5c6dd5df5a366d1c Mon Sep 17 00:00:00 2001 From: Mojtaba Date: Fri, 27 Mar 2026 13:10:48 +0100 Subject: [PATCH 3/5] address review comments for rotate3D --- src/pypulseq/rotate3D.py | 133 +++++++++++++++++++++++---------------- 1 file changed, 80 insertions(+), 53 deletions(-) diff --git a/src/pypulseq/rotate3D.py b/src/pypulseq/rotate3D.py index 6de225ae..83bb69a5 100644 --- a/src/pypulseq/rotate3D.py +++ b/src/pypulseq/rotate3D.py @@ -1,17 +1,17 @@ -import copy from types import SimpleNamespace -from typing import Any, List, Optional, Union - +from typing import List, Union, Any +import copy import numpy as np - from pypulseq.add_gradients import add_gradients from pypulseq.opts import Opts from pypulseq.scale_grad import scale_grad +from scipy.spatial.transform import Rotation + -def __get_grad_abs_mag(grad: SimpleNamespace) -> float: - """Magnitude used for thresholding output components (MATLAB getGradAbsMag).""" - if grad.type == 'trap': +def _get_grad_abs_mag(grad: SimpleNamespace) -> float: + """Magnitude used for thresholding output components.""" + if grad.type == "trap": return float(abs(grad.amplitude)) return float(np.max(np.abs(grad.waveform))) @@ -74,42 +74,29 @@ def align_gradient_to_raster(g: SimpleNamespace, raster: float) -> SimpleNamespa return g + def _quat_wxyz_to_rotmat(q_wxyz: np.ndarray) -> np.ndarray: """ Convert quaternion [w,x,y,z] (scalar first) to a 3x3 rotation matrix. - MATLAB: mr.aux.quat.toRotMat(rotation) """ q = np.asarray(q_wxyz, dtype=float).reshape(-1) if q.size != 4: raise ValueError('Quaternion must have length 4 [w, x, y, z].') - w, x, y, z = q - n = np.sqrt(w * w + x * x + y * y + z * z) - if n == 0: - raise ValueError('Quaternion norm must be non-zero.') - # normalize (MATLAB expects unit quaternion; we tolerate small numeric drift) - w, x, y, z = w / n, x / n, y / n, z / n + if np.linalg.norm(q) == 0: + raise ValueError("Quaternion norm must be non-zero.") - return np.array( - [ - [1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w)], - [2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w)], - [2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)], - ], - dtype=float, - ) + return Rotation.from_quat(q, scalar_first=True).as_matrix() + def _parse_rotation_to_matrix( *, - rotation: Optional[Any] = None, - rotation_matrix: Optional[np.ndarray] = None, + rotation: Union[Any, None] = None, + rotation_matrix: Union[np.ndarray, None] = None, ) -> np.ndarray: """ - MATLAB parity: - - rotation can be 3x3 matrix OR quaternion length-4 [w,x,y,z] - Backward compatibility: - - rotation_matrix keyword (your current API) still works. + rotation can be 3x3 matrix OR quaternion length-4 [w,x,y,z] """ if rotation is None and rotation_matrix is None: raise ValueError("You must provide either 'rotation' (matrix or quaternion) or 'rotation_matrix' (3x3).") @@ -131,23 +118,68 @@ def _parse_rotation_to_matrix( def rotate3D( *args: SimpleNamespace, - rotation: Optional[Any] = None, - rotation_matrix: Optional[np.ndarray] = None, + rotation: Union[Any, None] = None, + rotation_matrix: Union[np.ndarray, None] = None, system: Union[Opts, None] = None, ) -> List[SimpleNamespace]: """ - Rotate gradient events (trap/grad) by a 3x3 rotation matrix (or quaternion like MATLAB). - Non-gradient events (rf/adc/delay/etc) are passed through unchanged. - - Accepts mixed inputs like: - ms_rotate3D(gx, gy, rf, adc, rotation_matrix=R, system=sys) - ms_rotate3D(gx, gy, rf, adc, rotation=R, system=sys) - ms_rotate3D(gx, gy, rf, adc, rotation=[w,x,y,z], system=sys) # quaternion - - Also accepts list inputs: - ms_rotate3D([gx, gy, gz], rf, adc, rotation_matrix=R, system=sys) + Rotate gradient events by a 3*3 rotation matrix or quaternion. + + Non-gradient events (e.g. RF, ADC, delay) are passed through unchanged. + + Parameters + args : SimpleNamespace or list of SimpleNamespace + Input events. Can include gradient events ('grad', 'trap') and + non-gradient events (rf, adc, delay, etc.). Lists or tuples of events + are also accepted and will be flattened. + + rotation : array_like or None, Union + Rotation specified either as a 3*3 matrix or a quaternion + [w, x, y, z] (scalar-first convention). Default is None. + + rotation_matrix : np.ndarray or None, Union + Explicit 3*3 rotation matrix. Provided for backward compatibility. + If both 'rotation' and 'rotation_matrix' are given, 'rotation_matrix' + takes precedence. Default is None. + + system : Opts or None, Union + PyPulseq system limits. If None, 'Opts.default' is used. + + Returns + out : list of SimpleNamespace + List of events including rotated gradient components and unchanged + non-gradient events. + + Notes + Only one gradient per axis (x, y, z) is allowed in the input. + + Examples + >>> import numpy as np + >>> import pypulseq as pp + >>> from rotate3D import rotate3D + >>> + >>> system = pp.Opts() + >>> gx = pp.make_trapezoid(channel='x', system=system, amplitude=100, duration=2e-3) + >>> gy = pp.make_trapezoid(channel='y', system=system, amplitude=50, duration=2e-3) + >>> gz = pp.make_trapezoid(channel="z", system=system, amplitude=20, duration=2e-3) + >>> rf = pp.make_block_pulse(flip_angle=np.pi/2, duration=1e-3, system=system) + >>> adc = pp.make_adc(num_samples=64, dwell=10e-6, delay=0, system=system) + >>> + >>> angle = np.pi / 4 + >>> Rz = np.array([ + ... [np.cos(angle), -np.sin(angle), 0], + ... [np.sin(angle), np.cos(angle), 0], + ... [0, 0, 1], + ... ]) + >>> + >>> q = np.array([np.cos(angle / 2), 0.0, 0.0, np.sin(angle / 2)], dtype=float) + >>> + >>> out1 = rotate3D(gx, gy, gz, rf, adc, rotation_matrix=Rz, system=system) + >>> out2 = rotate3D(gx, gy, gz, rf, adc, rotation=Rz, system=system) + >>> out3 = rotate3D(gx, gy, gz, rf, adc, rotation=q, system=system) """ - # ---- flatten list/tuple inputs inside args ---- + + # flatten list/tuple inputs inside args flat_args: List[SimpleNamespace] = [] for a in args: if a is None: @@ -169,7 +201,7 @@ def rotate3D( events_to_rotate: dict[str, SimpleNamespace] = {} bypass: List[SimpleNamespace] = [] - # ---- split into gradients vs bypass ---- + # split into gradients vs bypass for event in args: if not hasattr(event, 'type') or event.type not in ('grad', 'trap'): bypass.append(event) @@ -181,18 +213,16 @@ def rotate3D( if event.channel in events_to_rotate: raise ValueError(f'More than one gradient for channel: {event.channel}') - # deepcopy to avoid mutating caller objects (matching your current code) events_to_rotate[event.channel] = copy.deepcopy(event) - # nothing to rotate -> just return bypass if len(events_to_rotate) == 0: return list(bypass) - # ---- thresholding like MATLAB ---- + # thresholding max_mag = 0.0 for ax in axes: if ax in events_to_rotate: - max_mag = max(max_mag, __get_grad_abs_mag(events_to_rotate[ax])) + max_mag = max(max_mag, _get_grad_abs_mag(events_to_rotate[ax])) fthresh = 1e-6 thresh = fthresh * max_mag @@ -214,15 +244,12 @@ def rotate3D( if grad_out_curr is None: grad_out_curr = scaled else: - # raster-align before summing (your current approach) - grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) - scaled = align_gradient_to_raster(scaled, system.grad_raster_time) - - # add_gradients expects a list + #grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) + #scaled = align_gradient_to_raster(scaled, system.grad_raster_time) grad_out_curr = add_gradients([grad_out_curr, scaled], system=system) + grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) - # only output non-zero amplitude gradients - if grad_out_curr is not None and __get_grad_abs_mag(grad_out_curr) >= thresh: + if grad_out_curr is not None and _get_grad_abs_mag(grad_out_curr) >= thresh: rotated_grads.append(grad_out_curr) return [*bypass, *rotated_grads] From 5461d32f289905ed2e0c9d8b499f29f25c610bef Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 27 Mar 2026 12:29:28 +0000 Subject: [PATCH 4/5] [pre-commit] auto fixes from pre-commit hooks --- src/pypulseq/rotate3D.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/pypulseq/rotate3D.py b/src/pypulseq/rotate3D.py index 83bb69a5..5f33cefc 100644 --- a/src/pypulseq/rotate3D.py +++ b/src/pypulseq/rotate3D.py @@ -1,17 +1,18 @@ -from types import SimpleNamespace -from typing import List, Union, Any import copy +from types import SimpleNamespace +from typing import Any, List, Union + import numpy as np +from scipy.spatial.transform import Rotation + from pypulseq.add_gradients import add_gradients from pypulseq.opts import Opts from pypulseq.scale_grad import scale_grad -from scipy.spatial.transform import Rotation - def _get_grad_abs_mag(grad: SimpleNamespace) -> float: """Magnitude used for thresholding output components.""" - if grad.type == "trap": + if grad.type == 'trap': return float(abs(grad.amplitude)) return float(np.max(np.abs(grad.waveform))) @@ -74,7 +75,6 @@ def align_gradient_to_raster(g: SimpleNamespace, raster: float) -> SimpleNamespa return g - def _quat_wxyz_to_rotmat(q_wxyz: np.ndarray) -> np.ndarray: """ Convert quaternion [w,x,y,z] (scalar first) to a 3x3 rotation matrix. @@ -84,10 +84,9 @@ def _quat_wxyz_to_rotmat(q_wxyz: np.ndarray) -> np.ndarray: raise ValueError('Quaternion must have length 4 [w, x, y, z].') if np.linalg.norm(q) == 0: - raise ValueError("Quaternion norm must be non-zero.") + raise ValueError('Quaternion norm must be non-zero.') return Rotation.from_quat(q, scalar_first=True).as_matrix() - def _parse_rotation_to_matrix( @@ -173,8 +172,8 @@ def rotate3D( ... ]) >>> >>> q = np.array([np.cos(angle / 2), 0.0, 0.0, np.sin(angle / 2)], dtype=float) - >>> - >>> out1 = rotate3D(gx, gy, gz, rf, adc, rotation_matrix=Rz, system=system) + >>> + >>> out1 = rotate3D(gx, gy, gz, rf, adc, rotation_matrix=Rz, system=system) >>> out2 = rotate3D(gx, gy, gz, rf, adc, rotation=Rz, system=system) >>> out3 = rotate3D(gx, gy, gz, rf, adc, rotation=q, system=system) """ @@ -244,8 +243,8 @@ def rotate3D( if grad_out_curr is None: grad_out_curr = scaled else: - #grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) - #scaled = align_gradient_to_raster(scaled, system.grad_raster_time) + # grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) + # scaled = align_gradient_to_raster(scaled, system.grad_raster_time) grad_out_curr = add_gradients([grad_out_curr, scaled], system=system) grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) From f39394554334e6cebc075db9bb45723cc0914ae9 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 14 Oct 2026 16:21:21 +0200 Subject: [PATCH 5/5] renaming, major refactoring, adding tests, ... --- src/pypulseq/__init__.py | 1 + src/pypulseq/rotate3D.py | 254 -------------------------------------- src/pypulseq/rotate_3d.py | 148 ++++++++++++++++++++++ tests/test_rotate_3d.py | 237 +++++++++++++++++++++++++++++++++++ 4 files changed, 386 insertions(+), 254 deletions(-) delete mode 100644 src/pypulseq/rotate3D.py create mode 100644 src/pypulseq/rotate_3d.py create mode 100644 tests/test_rotate_3d.py diff --git a/src/pypulseq/__init__.py b/src/pypulseq/__init__.py index 2aa0171e..5932976f 100644 --- a/src/pypulseq/__init__.py +++ b/src/pypulseq/__init__.py @@ -59,6 +59,7 @@ def round_half_up(n, decimals=0): from pypulseq.opts import Opts from pypulseq.points_to_waveform import points_to_waveform from pypulseq.rotate import rotate +from pypulseq.rotate_3d import rotate_3d from pypulseq.scale_grad import scale_grad from pypulseq.split_gradient import split_gradient from pypulseq.split_gradient_at import split_gradient_at diff --git a/src/pypulseq/rotate3D.py b/src/pypulseq/rotate3D.py deleted file mode 100644 index 5f33cefc..00000000 --- a/src/pypulseq/rotate3D.py +++ /dev/null @@ -1,254 +0,0 @@ -import copy -from types import SimpleNamespace -from typing import Any, List, Union - -import numpy as np -from scipy.spatial.transform import Rotation - -from pypulseq.add_gradients import add_gradients -from pypulseq.opts import Opts -from pypulseq.scale_grad import scale_grad - - -def _get_grad_abs_mag(grad: SimpleNamespace) -> float: - """Magnitude used for thresholding output components.""" - if grad.type == 'trap': - return float(abs(grad.amplitude)) - return float(np.max(np.abs(grad.waveform))) - - -def align_gradient_to_raster(g: SimpleNamespace, raster: float) -> SimpleNamespace: - """ - Align gradient timing parameters to the gradient raster and verify rounding error. - This is used to avoid add_gradients() complaining about non-raster time points. - """ - tol = 1e-9 # seconds - - if getattr(g, 'type', None) == 'grad': - original_tt = np.array(g.tt, dtype=np.float64) - rounded_tt = np.round(original_tt / raster) * raster - if np.any(np.abs(original_tt - rounded_tt) > tol): - raise ValueError(f"'grad' tt values not aligned to raster (>{tol}s): {original_tt} -> {rounded_tt}") - g.tt = rounded_tt - g.waveform = np.array(g.waveform, dtype=np.float64) - - elif getattr(g, 'type', None) == 'trap': - original_delay = float(getattr(g, 'delay', 0.0)) - original_rise = float(getattr(g, 'rise_time', 0.0)) - original_flat = float(getattr(g, 'flat_time', 0.0)) - original_fall = float(getattr(g, 'fall_time', 0.0)) - - rounded_delay = float(np.round(original_delay / raster) * raster) - rounded_rise = float(np.round(original_rise / raster) * raster) - rounded_flat = float(np.round(original_flat / raster) * raster) - rounded_fall = float(np.round(original_fall / raster) * raster) - - if any( - abs(x - y) > tol - for x, y in [ - (original_delay, rounded_delay), - (original_rise, rounded_rise), - (original_flat, rounded_flat), - (original_fall, rounded_fall), - ] - ): - raise ValueError( - f"'trap' timing values not aligned to raster (>{tol}s): " - f'delay={original_delay}→{rounded_delay}, ' - f'rise={original_rise}→{rounded_rise}, ' - f'flat={original_flat}→{rounded_flat}, ' - f'fall={original_fall}→{rounded_fall}' - ) - - g.delay = rounded_delay - g.rise_time = rounded_rise - g.flat_time = rounded_flat - g.fall_time = rounded_fall - - # always align delay if present - original_delay = float(getattr(g, 'delay', 0.0)) - rounded_delay = float(np.round(original_delay / raster) * raster) - if abs(original_delay - rounded_delay) > tol: - raise ValueError(f"'delay' value not aligned to raster (>{tol}s): {original_delay} → {rounded_delay}") - g.delay = rounded_delay - - return g - - -def _quat_wxyz_to_rotmat(q_wxyz: np.ndarray) -> np.ndarray: - """ - Convert quaternion [w,x,y,z] (scalar first) to a 3x3 rotation matrix. - """ - q = np.asarray(q_wxyz, dtype=float).reshape(-1) - if q.size != 4: - raise ValueError('Quaternion must have length 4 [w, x, y, z].') - - if np.linalg.norm(q) == 0: - raise ValueError('Quaternion norm must be non-zero.') - - return Rotation.from_quat(q, scalar_first=True).as_matrix() - - -def _parse_rotation_to_matrix( - *, - rotation: Union[Any, None] = None, - rotation_matrix: Union[np.ndarray, None] = None, -) -> np.ndarray: - """ - rotation can be 3x3 matrix OR quaternion length-4 [w,x,y,z] - """ - if rotation is None and rotation_matrix is None: - raise ValueError("You must provide either 'rotation' (matrix or quaternion) or 'rotation_matrix' (3x3).") - - if rotation_matrix is not None: - rot = np.asarray(rotation_matrix, dtype=float) - else: - rot = np.asarray(rotation, dtype=float) - - if rot.shape == (3, 3): - return rot - - # quaternion case - if rot.size == 4: - return _quat_wxyz_to_rotmat(rot) - - raise ValueError('rotation must be either a 3x3 matrix or a quaternion [w, x, y, z] (scalar first).') - - -def rotate3D( - *args: SimpleNamespace, - rotation: Union[Any, None] = None, - rotation_matrix: Union[np.ndarray, None] = None, - system: Union[Opts, None] = None, -) -> List[SimpleNamespace]: - """ - Rotate gradient events by a 3*3 rotation matrix or quaternion. - - Non-gradient events (e.g. RF, ADC, delay) are passed through unchanged. - - Parameters - args : SimpleNamespace or list of SimpleNamespace - Input events. Can include gradient events ('grad', 'trap') and - non-gradient events (rf, adc, delay, etc.). Lists or tuples of events - are also accepted and will be flattened. - - rotation : array_like or None, Union - Rotation specified either as a 3*3 matrix or a quaternion - [w, x, y, z] (scalar-first convention). Default is None. - - rotation_matrix : np.ndarray or None, Union - Explicit 3*3 rotation matrix. Provided for backward compatibility. - If both 'rotation' and 'rotation_matrix' are given, 'rotation_matrix' - takes precedence. Default is None. - - system : Opts or None, Union - PyPulseq system limits. If None, 'Opts.default' is used. - - Returns - out : list of SimpleNamespace - List of events including rotated gradient components and unchanged - non-gradient events. - - Notes - Only one gradient per axis (x, y, z) is allowed in the input. - - Examples - >>> import numpy as np - >>> import pypulseq as pp - >>> from rotate3D import rotate3D - >>> - >>> system = pp.Opts() - >>> gx = pp.make_trapezoid(channel='x', system=system, amplitude=100, duration=2e-3) - >>> gy = pp.make_trapezoid(channel='y', system=system, amplitude=50, duration=2e-3) - >>> gz = pp.make_trapezoid(channel="z", system=system, amplitude=20, duration=2e-3) - >>> rf = pp.make_block_pulse(flip_angle=np.pi/2, duration=1e-3, system=system) - >>> adc = pp.make_adc(num_samples=64, dwell=10e-6, delay=0, system=system) - >>> - >>> angle = np.pi / 4 - >>> Rz = np.array([ - ... [np.cos(angle), -np.sin(angle), 0], - ... [np.sin(angle), np.cos(angle), 0], - ... [0, 0, 1], - ... ]) - >>> - >>> q = np.array([np.cos(angle / 2), 0.0, 0.0, np.sin(angle / 2)], dtype=float) - >>> - >>> out1 = rotate3D(gx, gy, gz, rf, adc, rotation_matrix=Rz, system=system) - >>> out2 = rotate3D(gx, gy, gz, rf, adc, rotation=Rz, system=system) - >>> out3 = rotate3D(gx, gy, gz, rf, adc, rotation=q, system=system) - """ - - # flatten list/tuple inputs inside args - flat_args: List[SimpleNamespace] = [] - for a in args: - if a is None: - continue - if isinstance(a, (list, tuple)): - for x in a: - if x is not None: - flat_args.append(x) - else: - flat_args.append(a) - args = tuple(flat_args) - - if system is None: - system = Opts.default - - rotMat = _parse_rotation_to_matrix(rotation=rotation, rotation_matrix=rotation_matrix) - - axes = ['x', 'y', 'z'] - events_to_rotate: dict[str, SimpleNamespace] = {} - bypass: List[SimpleNamespace] = [] - - # split into gradients vs bypass - for event in args: - if not hasattr(event, 'type') or event.type not in ('grad', 'trap'): - bypass.append(event) - continue - - if not hasattr(event, 'channel') or event.channel not in axes: - raise ValueError(f'Invalid or missing gradient channel: {getattr(event, "channel", None)}') - - if event.channel in events_to_rotate: - raise ValueError(f'More than one gradient for channel: {event.channel}') - - events_to_rotate[event.channel] = copy.deepcopy(event) - - if len(events_to_rotate) == 0: - return list(bypass) - - # thresholding - max_mag = 0.0 - for ax in axes: - if ax in events_to_rotate: - max_mag = max(max_mag, _get_grad_abs_mag(events_to_rotate[ax])) - fthresh = 1e-6 - thresh = fthresh * max_mag - - rotated_grads: List[SimpleNamespace] = [] - for j in range(3): - grad_out_curr = None - for i in range(3): - ax_in = axes[i] - if ax_in not in events_to_rotate: - continue - - w = float(rotMat[j, i]) - if abs(w) < fthresh: - continue - - scaled = scale_grad(grad=events_to_rotate[ax_in], scale=w) - scaled.channel = axes[j] - - if grad_out_curr is None: - grad_out_curr = scaled - else: - # grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) - # scaled = align_gradient_to_raster(scaled, system.grad_raster_time) - grad_out_curr = add_gradients([grad_out_curr, scaled], system=system) - grad_out_curr = align_gradient_to_raster(grad_out_curr, system.grad_raster_time) - - if grad_out_curr is not None and _get_grad_abs_mag(grad_out_curr) >= thresh: - rotated_grads.append(grad_out_curr) - - return [*bypass, *rotated_grads] diff --git a/src/pypulseq/rotate_3d.py b/src/pypulseq/rotate_3d.py new file mode 100644 index 00000000..fd5b45c9 --- /dev/null +++ b/src/pypulseq/rotate_3d.py @@ -0,0 +1,148 @@ +from types import SimpleNamespace +from typing import List, Union + +import numpy as np +from scipy.spatial.transform import Rotation + +from pypulseq.add_gradients import add_gradients +from pypulseq.block_to_events import block_to_events +from pypulseq.make_delay import make_delay +from pypulseq.opts import Opts +from pypulseq.scale_grad import scale_grad +from pypulseq.utils.tracing import trace, trace_enabled + + +def _get_grad_abs_mag(grad: SimpleNamespace) -> float: + if grad.type == 'trap': + return abs(grad.amplitude) + return np.max(np.abs(grad.waveform)) + + +def _rotation_to_matrix(rotation: Union[np.ndarray, list, tuple, float, Rotation]) -> np.ndarray: + """Convert the supported rotation representations to a 3x3 rotation matrix.""" + if isinstance(rotation, Rotation): + return rotation.as_matrix() + + rotation = np.asarray(rotation, dtype=float) + if rotation.shape == (3, 3): + return rotation + + values = rotation.ravel() + if values.size == 4: + # Unit quaternion with the scalar component first [w, x, y, z]; scipy expects [x, y, z, w] + return Rotation.from_quat(values[[1, 2, 3, 0]]).as_matrix() + + if values.size in (1, 2): + phi = values[0] + if not -np.pi <= phi < 2 * np.pi: + raise ValueError(f'Rotation angle phi ({phi:.2f}) is invalid. Must be within [-pi, 2*pi) radians.') + rot = Rotation.from_euler('z', phi) + if values.size == 2: + theta = values[1] + if not -np.pi <= theta <= np.pi: + raise ValueError(f'Rotation angle theta ({theta:.2f}) is invalid. Must be within [-pi, pi] radians.') + rot = rot * Rotation.from_euler('y', theta) + return rot.as_matrix() + + raise ValueError( + 'The parameter `rotation` must be a 3x3 rotation matrix, a unit quaternion [w, x, y, z], ' + 'a pair of angles [phi, theta], a single angle phi or a scipy Rotation object.' + ) + + +def rotate_3d( + *args: SimpleNamespace, + rotation: Union[np.ndarray, list, tuple, float, Rotation], + system: Union[Opts, None] = None, +) -> List[SimpleNamespace]: + """ + Rotates the corresponding gradient(s) by the provided 3D rotation. Non-gradient events are not affected. + + The rotation can be given as a 3x3 rotation matrix, a unit quaternion with the scalar component first + [w, x, y, z], a pair of polar angles [phi, theta] in radians (rotation about y by theta followed by a rotation + about z by phi), a single angle phi in radians (rotation about z) or a `scipy.spatial.transform.Rotation` object. + For a single rotation about the z-axis `pypulseq.rotate()` can be faster. + + Gradients that end up on the same axis are added together with `pypulseq.add_gradients()`. Rotated gradient + components with negligible amplitude (relative to the largest input gradient) are dropped. Events that are not + modified are passed through unchanged, including their `id` attribute if present; the `id` attribute is removed + from the rotated gradients. + + See also `pypulseq.rotate.rotate()` and `pypulseq.Sequence.sequence.add_block()`. + + Parameters + ---------- + args : SimpleNamespace + Events to be rotated. Can also be a single block (as returned by `Sequence.get_block()`), in which case + the block duration is preserved as a delay event, or a single list of events, e.g. the output of a + previous `rotate_3d()` call. + rotation : numpy.ndarray, list, tuple, float or scipy.spatial.transform.Rotation + Rotation as a 3x3 matrix, a unit quaternion [w, x, y, z], polar angles [phi, theta] in radians, a single + angle phi in radians or a scipy `Rotation` object. + system : Opts, default=Opts() + System limits, used when adding gradient components on the same axis. + + Returns + ------- + rotated_grads : list[SimpleNamespace] + Non-gradient events followed by the rotated gradient(s). Can be passed directly to `Sequence.add_block()`. + + Raises + ------ + ValueError + If `rotation` has an unsupported shape or the angles are out of range. + If more than one gradient is provided for the same axis. + If a gradient has an invalid channel. + """ + if system is None: + system = Opts.default + + rot_mat = _rotation_to_matrix(rotation) + + events = block_to_events(*args) + if len(events) == 1 and isinstance(events[0], (list, tuple)): + events = events[0] # Allow nested calls, e.g. rotate_3d(rotate_3d(...), rotation=...) + + axes = ['x', 'y', 'z'] + grads_in = {} + bypass = [] + for event in events: + if event is None: + continue + if isinstance(event, (float, int)): + bypass.append(make_delay(event)) # block_duration of a block passed as input + continue + if getattr(event, 'type', None) not in ('grad', 'trap'): + bypass.append(event) + elif event.channel not in axes: + raise ValueError(f'Invalid gradient channel. Expected one of {axes}, got: {event.channel}') + elif event.channel in grads_in: + raise ValueError(f'More than one gradient on the same axis {event.channel} provided.') + else: + grads_in[event.channel] = event + + # Measure of the relevant amplitude, used to drop negligible rotated components + fthresh = 1e-6 + thresh = fthresh * max((_get_grad_abs_mag(g) for g in grads_in.values()), default=0) + + rotated_grads = [] + for j, out_axis in enumerate(axes): + components = [] + for i, in_axis in enumerate(axes): + if in_axis not in grads_in or abs(rot_mat[j, i]) < fthresh: + continue + g = scale_grad(grad=grads_in[in_axis], scale=rot_mat[j, i]) + g.channel = out_axis + components.append(g) + + if len(components) == 0: + continue + g = components[0] if len(components) == 1 else add_gradients(grads=components, system=system) + if _get_grad_abs_mag(g) >= thresh: + rotated_grads.append(g) + + if trace_enabled(): + for grad in rotated_grads: + grad.trace = trace() + + return [*bypass, *rotated_grads] diff --git a/tests/test_rotate_3d.py b/tests/test_rotate_3d.py new file mode 100644 index 00000000..4b2efdf2 --- /dev/null +++ b/tests/test_rotate_3d.py @@ -0,0 +1,237 @@ +from copy import copy +from types import SimpleNamespace + +import numpy as np +import pypulseq as pp +import pytest +from pypulseq.rotate_3d import rotate_3d +from scipy.spatial.transform import Rotation + +gx_trap = pp.make_trapezoid('x', area=1000, duration=1e-3) +gy_trap = pp.make_trapezoid('y', area=-500, duration=1e-3) +gz_trap = pp.make_trapezoid('z', area=200, duration=1e-3) +gx_ext = pp.make_extended_trapezoid('x', amplitudes=[0, 1e5, 5e4, 0], times=[0, 1e-4, 4e-4, 6e-4]) +gy_arb = pp.make_extended_trapezoid( + 'y', amplitudes=[0, 8e4, 8e4, 0], times=[0, 1e-4, 3e-4, 5e-4], convert_to_arbitrary=True +) +rf = pp.make_block_pulse(flip_angle=np.pi / 2, duration=1e-3) +adc = pp.make_adc(num_samples=100, duration=1e-3) +delay = pp.make_delay(2e-3) + + +def rot_mat(axis, angle): # standard right-handed rotation matrix + return Rotation.from_euler(axis, angle).as_matrix() + + +def by_channel(events): # dict channel -> gradient, non-gradients ignored + grads = {} + for event in events: + if getattr(event, 'type', None) not in ('grad', 'trap'): + continue + assert event.channel not in grads, f'duplicate channel {event.channel}' + grads[event.channel] = event + return grads + + +def assert_grads_equal(a, b, rel=1e-6, abs_=1e-6): + assert a.type == b.type and a.channel == b.channel + assert a.delay == pytest.approx(b.delay, abs=1e-12) + assert a.area == pytest.approx(b.area, rel=rel, abs=abs_) + if a.type == 'trap': + for f in ('amplitude', 'rise_time', 'flat_time', 'fall_time', 'flat_area'): + assert getattr(a, f) == pytest.approx(getattr(b, f), rel=rel, abs=abs_), f + else: + np.testing.assert_allclose(a.tt, b.tt, rtol=rel, atol=1e-12) + np.testing.assert_allclose(a.waveform, b.waveform, rtol=rel, atol=abs_) + for f in ('first', 'last', 'shape_dur'): + assert getattr(a, f) == pytest.approx(getattr(b, f), rel=rel, abs=abs_), f + + +def test_identity(): + amplitudes = [g.amplitude for g in (gx_trap, gy_trap, gz_trap)] + out = rotate_3d(gx_trap, gy_trap, gz_trap, rotation=np.eye(3)) + assert len(out) == 3 + assert [g.channel for g in out] == ['x', 'y', 'z'] + for rotated, original in zip(out, (gx_trap, gy_trap, gz_trap)): + assert_grads_equal(rotated, original) + # Inputs are not modified + assert [g.amplitude for g in (gx_trap, gy_trap, gz_trap)] == amplitudes + + +def test_rotation_input_forms(): + angle = 0.3 + c, s = np.cos(angle / 2), np.sin(angle / 2) + reference = by_channel(rotate_3d(gx_trap, gy_trap, rotation=rot_mat('z', angle))) + forms = [ + rot_mat('z', angle), + [c, 0, 0, s], + angle, + [angle, 0.0], + Rotation.from_euler('z', angle), + ] + for r in forms: + out = by_channel(rotate_3d(gx_trap, gy_trap, rotation=r)) + assert set(out) == set(reference) + for channel in reference: + assert_grads_equal(out[channel], reference[channel]) + + # [phi, theta] means rotation about y by theta followed by rotation about z by phi + phi, theta = 0.4, -0.7 + a = by_channel(rotate_3d(gx_trap, gy_trap, gz_trap, rotation=[phi, theta])) + b = by_channel(rotate_3d(gx_trap, gy_trap, gz_trap, rotation=rot_mat('z', phi) @ rot_mat('y', theta))) + assert set(a) == set(b) + for channel in a: + assert_grads_equal(a[channel], b[channel]) + + +@pytest.mark.parametrize('grads', [[g] for g in (gx_trap, gy_trap, gz_trap, gx_ext, gy_arb)] + [[gx_trap, gy_trap]]) +@pytest.mark.parametrize('angle', [0.0, 0.1, np.pi / 4, np.pi / 2, np.pi, -np.pi / 3, 2 * np.pi]) +@pytest.mark.parametrize('axis', ['x', 'y', 'z']) +def test_vs_rotate(axis, angle, grads): + # pypulseq's rotate() about the y-axis currently uses the opposite handedness to a right-handed + # rotation matrix (and to MATLAB's mr.rotate); compensate here until rotate() is fixed + sign = -1 if axis == 'y' else 1 + a = by_channel(pp.rotate(*[copy(g) for g in grads], axis=axis, angle=angle)) + b = by_channel(rotate_3d(*grads, rotation=rot_mat(axis, sign * angle))) + assert set(a) == set(b) + for channel in a: + assert_grads_equal(a[channel], b[channel], rel=1e-6, abs_=1e-3) + + +def test_oblique_rotation_and_inverse(): + R = Rotation.from_rotvec(0.9 * np.array([0.5, 0.3, 0.6]) / np.linalg.norm([0.5, 0.3, 0.6])).as_matrix() + out = rotate_3d(gx_trap, rotation=R) + assert len(out) == 3 + assert all(g.type == 'trap' for g in out) + for j, g in enumerate(out): + assert g.amplitude == pytest.approx(R[j, 0] * gx_trap.amplitude, rel=1e-9) + amplitudes = np.array([g.amplitude for g in out]) + assert np.linalg.norm(amplitudes) == pytest.approx(abs(gx_trap.amplitude), rel=1e-9) + + back = rotate_3d(*out, rotation=R.T) + assert len(back) == 1 + assert back[0].channel == 'x' + assert_grads_equal(back[0], gx_trap, abs_=1e-6 * abs(gx_trap.amplitude)) + + +def test_passthrough_and_id(): + gx = copy(gx_trap) + gx.id = 7 + rf_with_id = copy(rf) + rf_with_id.id = 3 + out = rotate_3d(rf_with_id, gx, adc, delay, rotation=rot_mat('z', np.pi / 2)) + assert out[0] is rf_with_id + assert out[1] is adc + assert out[2] is delay + assert out[3].type == 'trap' and out[3].channel == 'y' + assert not hasattr(out[3], 'id') + assert rf_with_id.id == 3 + assert gx.id == 7 + assert gx.amplitude == pytest.approx(gx_trap.amplitude) + + # None inputs are skipped + out = rotate_3d(None, gx_trap, None, rotation=np.eye(3)) + assert len(out) == 1 + assert_grads_equal(out[0], gx_trap) + + +def test_block_and_list_input(): + seq = pp.Sequence() + seq.add_block(rf, gx_trap, gz_trap, adc) + block = seq.get_block(1) + R = rot_mat('x', 0.5) + a = rotate_3d(block, rotation=R) + b = rotate_3d(rf, gx_trap, gz_trap, adc, rotation=R) + # The block's block_duration is preserved as one extra delay event + delays = [e for e in a if getattr(e, 'type', None) == 'delay'] + assert len(delays) == 1 + assert delays[0].delay == pytest.approx(block.block_duration) + a_events = [e for e in a if getattr(e, 'type', None) != 'delay'] + assert len(a_events) == len(b) + assert sorted(getattr(e, 'type', None) for e in a_events) == sorted(getattr(e, 'type', None) for e in b) + ca, cb = by_channel(a_events), by_channel(b) + assert set(ca) == set(cb) + for channel in ca: + assert_grads_equal(ca[channel], cb[channel]) + + # The rotated block can be passed back to add_block with the block duration preserved + seq2 = pp.Sequence() + seq2.add_block(*a) + assert seq2.block_durations[1] == pytest.approx(seq.block_durations[1]) + + # Nested call: rotate_3d(rotate_3d(...)) equals a single rotation by the combined matrix + R1, R2 = rot_mat('z', 0.4), rot_mat('x', -0.6) + nested = by_channel(rotate_3d(rotate_3d(gx_trap, gy_trap, rotation=R1), rotation=R2)) + combined = by_channel(rotate_3d(gx_trap, gy_trap, rotation=R2 @ R1)) + assert set(nested) == set(combined) + for channel in nested: + assert_grads_equal(nested[channel], combined[channel], abs_=1e-3) + + +def test_negligible_components_dropped(): + out = rotate_3d(gx_trap, rotation=rot_mat('z', np.pi / 2)) + assert len(out) == 1 + assert out[0].channel == 'y' + assert out[0].amplitude == pytest.approx(gx_trap.amplitude, rel=1e-6) + + out = rotate_3d(gz_trap, rotation=rot_mat('z', 1.0)) + assert len(out) == 1 + assert_grads_equal(out[0], gz_trap) + + +def test_mixed_gradient_types_are_added(): + out = by_channel(rotate_3d(gx_trap, gy_arb, rotation=rot_mat('z', np.pi / 4))) + assert set(out) == {'x', 'y'} + assert out['x'].type == 'grad' and out['y'].type == 'grad' + c, s = np.cos(np.pi / 4), np.sin(np.pi / 4) + assert out['x'].area == pytest.approx(c * gx_trap.area - s * gy_arb.area, rel=1e-6) + assert out['y'].area == pytest.approx(s * gx_trap.area + c * gy_arb.area, rel=1e-6) + + +def test_errors(): + with pytest.raises(ValueError, match='More than one gradient'): + rotate_3d(gx_trap, gx_trap, rotation=np.eye(3)) + with pytest.raises(ValueError, match='Invalid gradient channel'): + rotate_3d(SimpleNamespace(type='trap', channel='q'), rotation=np.eye(3)) + with pytest.raises(ValueError, match='must be a 3x3'): + rotate_3d(gx_trap, rotation=[1, 2, 3]) + with pytest.raises(ValueError, match='phi'): + rotate_3d(gx_trap, rotation=7.0) + with pytest.raises(ValueError, match='theta'): + rotate_3d(gx_trap, rotation=[0.0, 4.0]) + with pytest.raises(TypeError): + rotate_3d(gx_trap) + + +def test_system_limits_enforced(): + sys_lim = pp.Opts(max_grad=30, grad_unit='mT/m', max_slew=100, slew_unit='T/m/s') + g1 = pp.make_trapezoid('x', amplitude=0.95 * sys_lim.max_grad, flat_time=1e-3, system=sys_lim) + g2 = pp.make_trapezoid('y', amplitude=0.95 * sys_lim.max_grad, flat_time=1e-3, system=sys_lim) + with pytest.raises(ValueError): + rotate_3d(g1, g2, rotation=rot_mat('z', np.pi / 4), system=sys_lim) + + sys_relaxed = pp.Opts(max_grad=60, grad_unit='mT/m', max_slew=200, slew_unit='T/m/s') + out = rotate_3d(g1, g2, rotation=rot_mat('z', np.pi / 4), system=sys_relaxed) + # Equal-amplitude x and y gradients rotated by 45 deg cancel on x and combine to a single y trap + assert len(out) == 1 + assert out[0].type == 'trap' and out[0].channel == 'y' + assert out[0].amplitude == pytest.approx(np.sqrt(2) * 0.95 * sys_lim.max_grad, rel=1e-6) + + +def test_sequence_roundtrip(tmp_path): + seq = pp.Sequence() + seq.add_block(*rotate_3d(rf, gx_trap, gz_trap, rotation=Rotation.from_rotvec([0.3, -0.2, 0.5]).as_matrix())) + seq.add_block(*rotate_3d(gx_ext, gy_arb, rotation=[0.7, 0.2])) + ok, err = seq.check_timing() + assert ok, err + + seq.write(tmp_path / 'r.seq') + seq2 = pp.Sequence() + seq2.read(tmp_path / 'r.seq') + + for block_id in (1, 2): + a = by_channel(list(seq.get_block(block_id).__dict__.values())) + b = by_channel(list(seq2.get_block(block_id).__dict__.values())) + assert set(a) == set(b) + for channel in a: + assert a[channel].area == pytest.approx(b[channel].area, rel=1e-4)