diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 4b302d6..1052ba8 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -11,6 +11,7 @@ on:
permissions:
contents: write
+ id-token: write
jobs:
build:
@@ -71,10 +72,59 @@ jobs:
```
Requires MATLAB R2022a+. See [DEPLOY.md](matlab/DEPLOY.md) for details.
+ **Python:** Install from PyPI:
+ ```bash
+ pip install seqeyes-python
+ ```
+
## Changes
See [commits since last release](https://github.com/${{ github.repository }}/compare/${{ github.event.before }}...${{ github.ref_name }}).
files: |
seqeyes-web.vsix
seqeyes-*.mltbx
+ python/dist/*
generate_release_notes: true
+
+ python-build:
+ name: Build Python package
+ runs-on: ubuntu-latest
+ steps:
+ - name: Checkout
+ uses: actions/checkout@v4
+
+ - name: Setup Python
+ uses: actions/setup-python@v5
+ with:
+ python-version: '3.12'
+ cache: pip
+ cache-dependency-path: python/pyproject.toml
+
+ - name: Install build tools
+ run: |
+ python -m pip install --upgrade pip
+ pip install build
+
+ - name: Build wheel and sdist
+ run: cd python && python -m build
+
+ - name: Upload Python artifacts
+ uses: actions/upload-artifact@v4
+ with:
+ name: seqeyes-python-dist
+ path: python/dist/*
+
+ publish-pypi:
+ name: Publish to PyPI
+ needs: python-build
+ runs-on: ubuntu-latest
+ if: startsWith(github.ref, 'refs/tags/v')
+ steps:
+ - name: Download Python artifacts
+ uses: actions/download-artifact@v4
+ with:
+ name: seqeyes-python-dist
+ path: dist/
+
+ - name: Publish to PyPI
+ uses: pypa/gh-action-pypi-publish@release/v1
diff --git a/.github/workflows/python-tests.yml b/.github/workflows/python-tests.yml
new file mode 100644
index 0000000..2918c0c
--- /dev/null
+++ b/.github/workflows/python-tests.yml
@@ -0,0 +1,45 @@
+name: Python Package Tests
+
+on:
+ pull_request:
+ branches: [main, master]
+ push:
+ branches: [main, master, 'feature/**', 'fix/**']
+ workflow_dispatch:
+
+jobs:
+ test:
+ name: Python ${{ matrix.python-version }}
+ runs-on: ubuntu-latest
+ timeout-minutes: 15
+ strategy:
+ fail-fast: false
+ matrix:
+ python-version: ['3.10', '3.11', '3.12']
+
+ steps:
+ - name: Checkout
+ uses: actions/checkout@v5
+
+ - name: Setup Python ${{ matrix.python-version }}
+ uses: actions/setup-python@v5
+ with:
+ python-version: ${{ matrix.python-version }}
+ cache: pip
+ cache-dependency-path: python/pyproject.toml
+
+ - name: Install package with test dependencies
+ run: |
+ python -m pip install --upgrade pip
+ pip install -e "./python[test,pypulseq]"
+
+ - name: Run tests
+ run: python -m pytest python/tests/ -v
+
+ - name: Upload test results
+ uses: actions/upload-artifact@v5
+ if: ${{ !cancelled() }}
+ with:
+ name: python-test-results-py${{ matrix.python-version }}
+ path: python/.pytest_results/
+ retention-days: 14
diff --git a/.gitignore b/.gitignore
index be14e4b..62e2b35 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,8 +8,11 @@ dist/
*.mltbx
*.asc
.vscode-test/
+.vscode/
test-results/
playwright-report/
performance-results/
matlab/pulseq-bundle.js
.DS_Store
+__pycache__/
+.pytest_cache/
\ No newline at end of file
diff --git a/.vscode/launch.json b/.vscode/launch.json
deleted file mode 100644
index 1d975cc..0000000
--- a/.vscode/launch.json
+++ /dev/null
@@ -1,30 +0,0 @@
-{
- "version": "0.2.0",
- "configurations": [
- {
- "name": "Run Extension",
- "type": "extensionHost",
- "request": "launch",
- "args": [
- "--extensionDevelopmentPath=${workspaceFolder}"
- ],
- "outFiles": [
- "${workspaceFolder}/out/**/*.js"
- ],
- "preLaunchTask": "${defaultBuildTask}"
- },
- {
- "name": "Extension Tests",
- "type": "extensionHost",
- "request": "launch",
- "args": [
- "--extensionDevelopmentPath=${workspaceFolder}",
- "--extensionTestsPath=${workspaceFolder}/out/test/suite/index"
- ],
- "outFiles": [
- "${workspaceFolder}/out/test/**/*.js"
- ],
- "preLaunchTask": "${defaultBuildTask}"
- }
- ]
-}
diff --git a/.vscode/tasks.json b/.vscode/tasks.json
deleted file mode 100644
index 9a4fd8e..0000000
--- a/.vscode/tasks.json
+++ /dev/null
@@ -1,24 +0,0 @@
-{
- "version": "2.0.0",
- "tasks": [
- {
- "type": "npm",
- "script": "watch",
- "problemMatcher": "$tsc-watch",
- "isBackground": true,
- "presentation": {
- "reveal": "never"
- },
- "group": {
- "kind": "build",
- "isDefault": true
- }
- },
- {
- "type": "npm",
- "script": "compile",
- "problemMatcher": "$tsc",
- "group": "build"
- }
- ]
-}
diff --git a/README.md b/README.md
index aa1b9e3..5fc8fe3 100644
--- a/README.md
+++ b/README.md
@@ -73,6 +73,41 @@ open('spiral_inout.seq') % or doubleβclick in Current Folder
All the same features as the browser & VS Code versions β 7 channels, kβspace viewer, themes, tooltips β rendered inside a native MATLAB figure. Requires R2022a+.
+## π Python Package
+
+Interactive Pulseq sequence viewer for Jupyter notebooks and Python scripts β a dropβin replacement for `pypulseq.Sequence.plot()`. Renders directly in notebook cell output or opens in your default browser.
+
+### Install
+
+```bash
+pip install seqeyes-python
+```
+
+For pypulseq integration:
+
+```bash
+pip install seqeyes-python[pypulseq]
+```
+
+### Usage
+
+```python
+import seqeyes
+
+# Enable SeqEyes (once per session) β seq.plot() is now interactive
+seqeyes.set(theme="dark", time_disp="ms")
+
+# Build your sequence with pypulseq as usual
+seq.plot() # interactive viewer in Jupyter
+seq.plot(show_blocks=True) # perβcall overrides
+seq.plot(time_range=(0, 0.05)) # zoom to first 50 ms
+
+# Restore matplotlib at any time
+seqeyes.reset()
+```
+
+All the same features as the other versions β interactive waveforms, kβspace viewer, themes, tooltips β rendered directly in Jupyter or your browser. Requires Python β₯ 3.9.
+
## Features
- **Custom editor for `.seq` files** β opens automatically on doubleβclick
diff --git a/package.json b/package.json
index 86d09b3..6e4d2fc 100644
--- a/package.json
+++ b/package.json
@@ -2,7 +2,7 @@
"name": "seqeyes-web",
"displayName": "SeqEyes",
"description": "Visualize Pulseq MRI sequences inside VS Code β inspect sequence diagrams, k-space trajectorieswith interactive user interface.",
- "version": "0.1.17",
+ "version": "0.2.0",
"publisher": "SeqEyesDeveloper",
"license": "MIT",
"icon": "images/logo.png",
diff --git a/python/README.md b/python/README.md
new file mode 100644
index 0000000..e6406ef
--- /dev/null
+++ b/python/README.md
@@ -0,0 +1,94 @@
+# SeqEyes β Interactive Pulseq MRI Sequence Viewer for Python
+
+**SeqEyes** is a lightweight Python package that provides interactive
+visualization of Pulseq (.seq) MRI sequences in Jupyter notebooks.
+It works as a dropβin replacement for `pypulseq.Sequence.plot()`,
+rendering an interactive viewer directly in Jupyter notebook cell
+output β just like Plotly.
+
+## Features
+
+- ποΈ **Interactive waveform viewer** β zoom, pan, perβchannel amplitude zoom
+- π **Tooltip** with block details (RF amplitude, gradient strength, ADC params)
+- πΊοΈ **3D kβspace trajectory viewer** β rotate, zoom, depthβsorted rendering
+- π¨ **8 colour themes** β system, light, dark, dracula, nord, and more
+- π **Unit conversion** β time (s / ms / Β΅s), gradient (Hz/m / mT/m / G/cm)
+- π **Minimap** with TR/TE overlay and viewport indicator
+- πΎ **Export to standalone HTML** β shareable, no Python needed
+- π **Dropβin pypulseq integration** β `seq.plot()` just works
+
+## Installation
+
+```bash
+pip install seqeyes-python
+```
+
+For pypulseq integration:
+```bash
+pip install seqeyes-python[pypulseq]
+```
+
+## Quick Start
+
+```python
+import seqeyes
+
+# Enable SeqEyes (once per session) β seq.plot() is now interactive
+seqeyes.set(theme="dark", time_disp="ms")
+
+# Build your sequence with pypulseq as usual
+seq.plot() # interactive viewer in Jupyter
+seq.plot(show_blocks=True) # perβcall overrides
+seq.plot(time_range=(0, 0.05)) # zoom to first 50 ms
+
+# Restore matplotlib at any time
+seqeyes.reset()
+```
+
+In a plain `.py` script (no Jupyter), `seq.plot()` opens the viewer
+in a desktop popβup window (requires `pywebview`) or falls back to
+your default browser.
+
+### Using without pypulseq
+
+```python
+from seqeyes import SeqEyesViewer
+
+with open('my_sequence.seq') as f:
+ viewer = SeqEyesViewer(f.read(), theme="dark")
+
+viewer # renders inline in Jupyter
+```
+
+## API Reference
+
+| Function | Description |
+|---|---|
+| `seqeyes.set(**kwargs)` | Enable SeqEyes and set global defaults (`theme`, `show_blocks`, `time_disp`, `grad_disp`, `time_range`) |
+| `seqeyes.reset()` | Restore matplotlib `seq.plot()` and clear all defaults |
+| `SeqEyesViewer(seq_text, ...)` | Lowβlevel viewer for raw `.seq` content (no pypulseq needed) |
+
+## Viewer Controls
+
+| Action | How |
+|---|---|
+| Zoom | Scroll wheel |
+| Pan | Click + drag |
+| Amplitude zoom (per channel) | Ctrl + scroll wheel |
+| Tooltip | Hover over waveform |
+| Toggle channels | Click legend labels |
+| KβSpace viewer | Click "KβSpace" button |
+| Rotate kβspace | Click + drag in panel |
+| Minimap navigation | Click on minimap strip |
+| Open another file | π Open button |
+
+## Requirements
+
+- Python β₯ 3.9
+- numpy β₯ 1.21
+- pypulseq β₯ 1.4 (optional, for `seq.plot()` integration)
+- pywebview β₯ 5 (optional, for native desktop popβup windows)
+
+## License
+
+MIT
diff --git a/python/examples/demo.ipynb b/python/examples/demo.ipynb
new file mode 100644
index 0000000..370aa9d
--- /dev/null
+++ b/python/examples/demo.ipynb
@@ -0,0 +1,394 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "e697bb42",
+ "metadata": {},
+ "source": [
+ "# SeqEyes β Interactive Pulseq MRI Sequence Viewer\n",
+ "\n",
+ "Dropβin replacement for `pypulseq.Sequence.plot()`.\n",
+ "\n",
+ "```python\n",
+ "import seqeyes\n",
+ "seqeyes.set(theme=\"dark\") # enable + set defaults\n",
+ "\n",
+ "seq.plot() # interactive viewer inline\n",
+ "seq.plot(show_blocks=True) # perβcall overrides\n",
+ "seq.plot(time_range=(0, 0.05)) # zoom to first 50 ms\n",
+ "\n",
+ "seqeyes.reset() # back to matplotlib\n",
+ "```\n",
+ "\n",
+ "**Plain `.py` script?** `seq.plot()` opens your browser automatically."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "699af581",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "d:\\MIniconda\\envs\\NUM\\Lib\\site-packages\\sigpy\\config.py:27: UserWarning: Importing cupy.cuda.cudnn failed. For more details, see the error stack below:\n",
+ "DLL load failed while importing cudnn: The specified module could not be found.\n",
+ " warnings.warn(\n"
+ ]
+ },
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "SeqEyes 0.1.0 β seq.plot() is now interactive!\n",
+ "GRE: 320 blocks | EPI: 390 blocks | Radial: 405 blocks\n"
+ ]
+ }
+ ],
+ "source": [
+ "import seqeyes\n",
+ "import numpy as np\n",
+ "import pypulseq as pp\n",
+ "\n",
+ "# ββ Enable SeqEyes (once per session) ββ\n",
+ "seqeyes.set(time_disp=\"ms\", grad_disp=\"kHz/m\")\n",
+ "print(f\"SeqEyes {seqeyes.__version__} β seq.plot() is now interactive!\")\n",
+ "\n",
+ "# βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ\n",
+ "# 1. GRE (from pypulseq/examples/scripts/write_gre.py)\n",
+ "# βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ\n",
+ "fov_gre = 256e-3; n_x_gre = 64; n_y_gre = 64\n",
+ "flip_gre = 10; slice_thickness = 3e-3; tr_gre = 12e-3; te_gre = 5e-3\n",
+ "\n",
+ "system_gre = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=150, slew_unit='T/m/s',\n",
+ " rf_ringdown_time=20e-6, rf_dead_time=100e-6, adc_dead_time=10e-6)\n",
+ "seq_gre = pp.Sequence(system_gre)\n",
+ "\n",
+ "rf, gz, _ = pp.make_sinc_pulse(flip_angle=np.deg2rad(flip_gre), duration=3e-3,\n",
+ " slice_thickness=slice_thickness, apodization=0.42,\n",
+ " time_bw_product=4, system=system_gre, return_gz=True,\n",
+ " delay=system_gre.rf_dead_time, use='excitation')\n",
+ "\n",
+ "delta_kx = 1 / fov_gre; delta_ky = 1 / fov_gre\n",
+ "gx = pp.make_trapezoid(channel='x', flat_area=n_x_gre * delta_kx, flat_time=3.2e-3, system=system_gre)\n",
+ "adc = pp.make_adc(num_samples=n_x_gre, duration=gx.flat_time, delay=gx.rise_time, system=system_gre)\n",
+ "gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2, duration=1e-3, system=system_gre)\n",
+ "gz_reph = pp.make_trapezoid(channel='z', area=-gz.area / 2, duration=1e-3, system=system_gre)\n",
+ "gx_spoil = pp.make_trapezoid(channel='x', area=2 * n_x_gre * delta_kx, system=system_gre)\n",
+ "gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system_gre)\n",
+ "phase_areas = (np.arange(n_y_gre) - n_y_gre / 2) * delta_ky\n",
+ "\n",
+ "te_delay = te_gre - (pp.calc_duration(gz, rf) - pp.calc_rf_center(rf)[0] - rf.delay) - pp.calc_duration(gx_pre) - pp.calc_duration(gx) / 2 - pp.eps\n",
+ "te_delay = np.ceil(te_delay / seq_gre.grad_raster_time) * seq_gre.grad_raster_time\n",
+ "tr_delay = tr_gre - pp.calc_duration(gz, rf) - pp.calc_duration(gx_pre) - pp.calc_duration(gx) - te_delay\n",
+ "tr_delay = np.ceil(tr_delay / seq_gre.grad_raster_time) * seq_gre.grad_raster_time\n",
+ "\n",
+ "rf_phase = 0; rf_inc = 0\n",
+ "for i_phase in range(n_y_gre):\n",
+ " rf.phase_offset = rf_phase / 180 * np.pi\n",
+ " adc.phase_offset = rf_phase / 180 * np.pi\n",
+ " rf_inc = divmod(rf_inc + 117, 360.0)[1]\n",
+ " rf_phase = divmod(rf_phase + rf_inc, 360.0)[1]\n",
+ " seq_gre.add_block(rf, gz)\n",
+ " gy_pre = pp.make_trapezoid(channel='y', area=phase_areas[i_phase], duration=pp.calc_duration(gx_pre), system=system_gre)\n",
+ " seq_gre.add_block(gx_pre, gy_pre, gz_reph)\n",
+ " seq_gre.add_block(pp.make_delay(te_delay))\n",
+ " seq_gre.add_block(gx, adc)\n",
+ " gy_pre.amplitude = -gy_pre.amplitude\n",
+ " seq_gre.add_block(pp.make_delay(tr_delay), gx_spoil, gy_pre, gz_spoil)\n",
+ "\n",
+ "# βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ\n",
+ "# 2. EPI (from pypulseq/examples/scripts/write_epi.py)\n",
+ "# βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ\n",
+ "fov_epi = 220e-3; n_x_epi = 64; n_y_epi = 64; n_slices = 3\n",
+ "\n",
+ "system_epi = pp.Opts(max_grad=32, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s',\n",
+ " rf_ringdown_time=30e-6, rf_dead_time=100e-6)\n",
+ "seq_epi = pp.Sequence(system_epi)\n",
+ "\n",
+ "rf_epi, gz_epi, _ = pp.make_sinc_pulse(flip_angle=np.pi / 2, system=system_epi, duration=3e-3,\n",
+ " slice_thickness=slice_thickness, apodization=0.5,\n",
+ " time_bw_product=4, return_gz=True,\n",
+ " delay=system_epi.rf_dead_time, use='excitation')\n",
+ "\n",
+ "dkx = 1 / fov_epi; dky = 1 / fov_epi; k_width = n_x_epi * dkx\n",
+ "adc_dwell = 4e-6; adc_duration = n_x_epi * adc_dwell\n",
+ "gx_flat = np.ceil(adc_duration * 1e5) * 1e-5\n",
+ "gx_epi = pp.make_trapezoid(channel='x', system=system_epi, amplitude=k_width / adc_duration, flat_time=gx_flat)\n",
+ "adc_epi = pp.make_adc(num_samples=n_x_epi, duration=adc_duration,\n",
+ " delay=gx_epi.rise_time + gx_flat / 2 - (adc_duration - adc_dwell) / 2)\n",
+ "pre_time = 8e-4\n",
+ "gx_pre_epi = pp.make_trapezoid(channel='x', system=system_epi, area=-gx_epi.area / 2, duration=pre_time)\n",
+ "gz_reph_epi = pp.make_trapezoid(channel='z', system=system_epi, area=-gz_epi.area / 2, duration=pre_time)\n",
+ "gy_pre_epi = pp.make_trapezoid(channel='y', system=system_epi, area=-n_y_epi / 2 * dky, duration=pre_time)\n",
+ "gy_blip_dur = np.ceil(2 * np.sqrt(dky / system_epi.max_slew) / 10e-6) * 10e-6\n",
+ "gy_epi = pp.make_trapezoid(channel='y', system=system_epi, area=dky, duration=gy_blip_dur)\n",
+ "\n",
+ "for i_slice in range(n_slices):\n",
+ " rf_epi.freq_offset = gz_epi.amplitude * slice_thickness * (i_slice - (n_slices - 1) / 2)\n",
+ " seq_epi.add_block(rf_epi, gz_epi)\n",
+ " seq_epi.add_block(gx_pre_epi, gy_pre_epi, gz_reph_epi)\n",
+ " for _ in range(n_y_epi):\n",
+ " seq_epi.add_block(gx_epi, adc_epi)\n",
+ " seq_epi.add_block(gy_epi)\n",
+ " gx_epi.amplitude = -gx_epi.amplitude\n",
+ "\n",
+ "# βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ\n",
+ "# 3. Radial GRE (from pypulseq/examples/scripts/write_radial_gre.py)\n",
+ "# βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ\n",
+ "fov_rad = 260e-3; n_x_rad = 64; flip_rad = 10; n_spokes = 60; n_dummy = 20\n",
+ "tr_rad = 20e-3; te_rad = 8e-3\n",
+ "\n",
+ "system_rad = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=120, slew_unit='T/m/s',\n",
+ " rf_ringdown_time=20e-6, rf_dead_time=100e-6, adc_dead_time=10e-6)\n",
+ "seq_rad = pp.Sequence(system_rad)\n",
+ "\n",
+ "rf_rad, gz_rad, _ = pp.make_sinc_pulse(apodization=0.5, duration=4e-3,\n",
+ " flip_angle=np.deg2rad(flip_rad),\n",
+ " slice_thickness=slice_thickness, system=system_rad,\n",
+ " time_bw_product=4, return_gz=True,\n",
+ " delay=system_rad.rf_dead_time, use='excitation')\n",
+ "\n",
+ "dkx_rad = 1 / fov_rad\n",
+ "gx_rad = pp.make_trapezoid(channel='x', flat_area=n_x_rad * dkx_rad, flat_time=6.4e-3 / 5, system=system_rad)\n",
+ "adc_rad = pp.make_adc(num_samples=n_x_rad, duration=gx_rad.flat_time, delay=gx_rad.rise_time, system=system_rad)\n",
+ "gx_pre_rad = pp.make_trapezoid(channel='x', area=-gx_rad.area / 2 - dkx_rad / 2, duration=2e-3, system=system_rad)\n",
+ "gz_reph_rad = pp.make_trapezoid(channel='z', area=-gz_rad.area / 2, duration=2e-3, system=system_rad)\n",
+ "gx_spoil_rad = pp.make_trapezoid(channel='x', area=0.5 * n_x_rad * dkx_rad, system=system_rad)\n",
+ "gz_spoil_rad = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system_rad)\n",
+ "\n",
+ "te_delay_rad = te_rad - pp.calc_duration(gx_pre_rad) - gz_rad.fall_time - gz_rad.flat_time / 2 - pp.calc_duration(gx_rad) / 2\n",
+ "te_delay_rad = np.ceil(te_delay_rad / seq_rad.grad_raster_time) * seq_rad.grad_raster_time\n",
+ "tr_delay_rad = tr_rad - pp.calc_duration(gx_pre_rad) - pp.calc_duration(gz_rad) - pp.calc_duration(gx_rad) - te_delay_rad\n",
+ "tr_delay_rad = np.ceil(tr_delay_rad / seq_rad.grad_raster_time) * seq_rad.grad_raster_time\n",
+ "\n",
+ "spoke_angle = np.pi / n_spokes\n",
+ "rf_phase = 0; rf_inc = 0\n",
+ "for i_spoke in range(-n_dummy, n_spokes + 1):\n",
+ " rf_rad.phase_offset = rf_phase / 180 * np.pi\n",
+ " adc_rad.phase_offset = rf_phase / 180 * np.pi\n",
+ " rf_inc = divmod(rf_inc + 117, 360.0)[1]\n",
+ " rf_phase = divmod(rf_inc + rf_phase, 360.0)[1]\n",
+ " seq_rad.add_block(rf_rad, gz_rad)\n",
+ " phi = spoke_angle * (i_spoke - 1)\n",
+ " seq_rad.add_block(*pp.rotate(gx_pre_rad, gz_reph_rad, angle=phi, axis='z'))\n",
+ " seq_rad.add_block(pp.make_delay(te_delay_rad))\n",
+ " if i_spoke > 0:\n",
+ " seq_rad.add_block(*pp.rotate(gx_rad, adc_rad, angle=phi, axis='z'))\n",
+ " else:\n",
+ " seq_rad.add_block(*pp.rotate(gx_rad, angle=phi, axis='z'))\n",
+ " seq_rad.add_block(*pp.rotate(gx_spoil_rad, gz_spoil_rad, pp.make_delay(tr_delay_rad), angle=phi, axis='z'))\n",
+ "\n",
+ "print(f\"GRE: {len(seq_gre.block_events)} blocks | \"\n",
+ " f\"EPI: {len(seq_epi.block_events)} blocks | \"\n",
+ " f\"Radial: {len(seq_rad.block_events)} blocks\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "304160f9",
+ "metadata": {},
+ "source": [
+ "## 1. Gradient Echo β `write_gre.py`\n",
+ "\n",
+ "Classic Cartesian GRE with sinc excitation, phase encoding, and\n",
+ "gradient spoiling. 64Γ64 matrix, TE=5 ms, TR=12 ms."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "9c2645c0",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "d:\\MIniconda\\envs\\NUM\\Lib\\site-packages\\IPython\\core\\display.py:447: UserWarning: Consider using IPython.display.IFrame instead\n",
+ " warnings.warn(\"Consider using IPython.display.IFrame instead\")\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "seq_gre.plot(time_range=(0, tr_gre), time_disp=\"ms\", grad_disp=\"kHz/m\", theme=\"dark\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "31cffaa6",
+ "metadata": {},
+ "source": [
+ "## 2. Echo Planar Imaging β `write_epi.py`\n",
+ "\n",
+ "Singleβshot EPI with blipped phase encoding, alternating readout\n",
+ "polarity. 64Γ64, 3 slices, 220 mm FOV."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "4fa512a7",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "seq_epi.plot(time_disp=\"ms\", grad_disp=\"kHz/m\", theme=\"dark\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fdb2b647",
+ "metadata": {},
+ "source": [
+ "## 3. Radial GRE β `write_radial_gre.py`\n",
+ "\n",
+ "NonβCartesian radial trajectory with 60 spokes, rotated readout\n",
+ "gradients, and RF spoiling. Open the **KβSpace** panel!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "2750f1cf",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "seq_rad.plot(time_disp=\"ms\", grad_disp=\"kHz/m\", theme=\"dark\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0dc32d72",
+ "metadata": {},
+ "source": [
+ "## Customize the view\n",
+ "\n",
+ "`seq.plot()` accepts the same arguments as pypulseq's original plot,\n",
+ "plus a `theme` keyword:\n",
+ "\n",
+ "| Argument | Default | Description |\n",
+ "|---|---|---|\n",
+ "| `show_blocks` | `False` | Show blockβboundary lines |\n",
+ "| `time_range` | `(0, inf)` | Zoom to `(start_sec, end_sec)` |\n",
+ "| `time_disp` | `\"s\"` | Time unit: `\"s\"`, `\"ms\"`, `\"us\"` |\n",
+ "| `grad_disp` | `\"kHz/m\"` | Gradient unit: `\"Hz/m\"`, `\"kHz/m\"`, `\"mT/m\"`, `\"G/cm\"` |\n",
+ "| `theme` | `\"system\"` | `\"light\"`, `\"dark\"`, `\"dracula\"`, `\"nord\"`, ... |"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "a8ef3961",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# Zoom to first TR with block boundaries + custom units\n",
+ "seqeyes.set(show_blocks=True, time_disp=\"ms\", grad_disp=\"kHz/m\")\n",
+ "seq_gre.plot(time_range=(0, tr_gre), theme=\"dark\")\n",
+ "seqeyes.reset() # restore defaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ec80c99c",
+ "metadata": {},
+ "source": [
+ "## Viewer Controls\n",
+ "\n",
+ "| Action | How |\n",
+ "|---|---|\n",
+ "| **Zoom** | Scroll wheel |\n",
+ "| **Pan** | Click + drag |\n",
+ "| **Amplitude zoom (per channel)** | Ctrl + scroll wheel |\n",
+ "| **Tooltip** | Hover over waveform |\n",
+ "| **Toggle channels** | Click legend labels |\n",
+ "| **Block boundaries** | Check \"Blocks\" in toolbar |\n",
+ "| **KβSpace 3D viewer** | Click \"KβSpace\" button |\n",
+ "| **Rotate kβspace** | Click + drag in panel |\n",
+ "| **Change time unit** | Dropdown (s / ms / Β΅s) |\n",
+ "| **Change gradient unit** | Dropdown (Hz/m / mT/m / G/cm) |\n",
+ "| **Change theme** | Dropdown in toolbar |\n",
+ "| **Open another .seq** | π Open button |\n",
+ "\n",
+ "## Plain `.py` script?\n",
+ "\n",
+ "In a standalone Python script (no Jupyter), `seq.plot()` opens the\n",
+ "viewer in your default web browser. See `examples/demo.py`."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "NUM",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.11.13"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/python/examples/demo.py b/python/examples/demo.py
new file mode 100644
index 0000000..2ecbc40
--- /dev/null
+++ b/python/examples/demo.py
@@ -0,0 +1,153 @@
+#!/usr/bin/env python3
+"""
+demo.py β SeqEyes with 3 official pypulseq example sequences.
+Run: python demo.py
+Each ``seq.plot()`` opens an interactive viewer in your default browser.
+"""
+
+import seqeyes
+import numpy as np
+import pypulseq as pp
+
+# ββ Enable SeqEyes (once per session) βββββββββββββββββββββββββββββββββ
+seqeyes.set(time_disp="ms", grad_disp="kHz/m")
+print(f"SeqEyes {seqeyes.__version__} β seq.plot() is now interactive!")
+
+
+# ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+# 1. GRE (from pypulseq/examples/scripts/write_gre.py)
+# ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+fov_gre = 256e-3; n_x_gre = 64; n_y_gre = 64
+flip_gre = 10; slice_thickness = 3e-3; tr_gre = 12e-3; te_gre = 5e-3
+
+system_gre = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=150, slew_unit='T/m/s',
+ rf_ringdown_time=20e-6, rf_dead_time=100e-6, adc_dead_time=10e-6)
+seq_gre = pp.Sequence(system_gre)
+
+rf, gz, _ = pp.make_sinc_pulse(flip_angle=np.deg2rad(flip_gre), duration=3e-3,
+ slice_thickness=slice_thickness, apodization=0.42,
+ time_bw_product=4, system=system_gre, return_gz=True,
+ delay=system_gre.rf_dead_time, use='excitation')
+
+delta_kx = 1 / fov_gre; delta_ky = 1 / fov_gre
+gx = pp.make_trapezoid(channel='x', flat_area=n_x_gre * delta_kx, flat_time=3.2e-3, system=system_gre)
+adc = pp.make_adc(num_samples=n_x_gre, duration=gx.flat_time, delay=gx.rise_time, system=system_gre)
+gx_pre = pp.make_trapezoid(channel='x', area=-gx.area / 2, duration=1e-3, system=system_gre)
+gz_reph = pp.make_trapezoid(channel='z', area=-gz.area / 2, duration=1e-3, system=system_gre)
+gx_spoil = pp.make_trapezoid(channel='x', area=2 * n_x_gre * delta_kx, system=system_gre)
+gz_spoil = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system_gre)
+phase_areas = (np.arange(n_y_gre) - n_y_gre / 2) * delta_ky
+
+te_delay = te_gre - (pp.calc_duration(gz, rf) - pp.calc_rf_center(rf)[0] - rf.delay) - pp.calc_duration(gx_pre) - pp.calc_duration(gx) / 2 - pp.eps
+te_delay = np.ceil(te_delay / seq_gre.grad_raster_time) * seq_gre.grad_raster_time
+tr_delay = tr_gre - pp.calc_duration(gz, rf) - pp.calc_duration(gx_pre) - pp.calc_duration(gx) - te_delay
+tr_delay = np.ceil(tr_delay / seq_gre.grad_raster_time) * seq_gre.grad_raster_time
+
+rf_phase = 0; rf_inc = 0
+for i_phase in range(n_y_gre):
+ rf.phase_offset = rf_phase / 180 * np.pi
+ adc.phase_offset = rf_phase / 180 * np.pi
+ rf_inc = divmod(rf_inc + 117, 360.0)[1]
+ rf_phase = divmod(rf_phase + rf_inc, 360.0)[1]
+ seq_gre.add_block(rf, gz)
+ gy_pre = pp.make_trapezoid(channel='y', area=phase_areas[i_phase], duration=pp.calc_duration(gx_pre), system=system_gre)
+ seq_gre.add_block(gx_pre, gy_pre, gz_reph)
+ seq_gre.add_block(pp.make_delay(te_delay))
+ seq_gre.add_block(gx, adc)
+ gy_pre.amplitude = -gy_pre.amplitude
+ seq_gre.add_block(pp.make_delay(tr_delay), gx_spoil, gy_pre, gz_spoil)
+
+print("1/3 GRE sequence β opening viewer ...")
+seq_gre.plot(time_range=(0, 0.05)) # zoom to first 50 ms
+
+
+# ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+# 2. EPI (from pypulseq/examples/scripts/write_epi.py)
+# ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+fov_epi = 220e-3; n_x_epi = 64; n_y_epi = 64; n_slices = 3
+
+system_epi = pp.Opts(max_grad=32, grad_unit='mT/m', max_slew=130, slew_unit='T/m/s',
+ rf_ringdown_time=30e-6, rf_dead_time=100e-6)
+seq_epi = pp.Sequence(system_epi)
+
+rf_epi, gz_epi, _ = pp.make_sinc_pulse(flip_angle=np.pi / 2, system=system_epi, duration=3e-3,
+ slice_thickness=slice_thickness, apodization=0.5,
+ time_bw_product=4, return_gz=True,
+ delay=system_epi.rf_dead_time, use='excitation')
+
+dkx = 1 / fov_epi; dky = 1 / fov_epi; k_width = n_x_epi * dkx
+adc_dwell = 4e-6; adc_duration = n_x_epi * adc_dwell
+gx_flat = np.ceil(adc_duration * 1e5) * 1e-5
+gx_epi = pp.make_trapezoid(channel='x', system=system_epi, amplitude=k_width / adc_duration, flat_time=gx_flat)
+adc_epi = pp.make_adc(num_samples=n_x_epi, duration=adc_duration,
+ delay=gx_epi.rise_time + gx_flat / 2 - (adc_duration - adc_dwell) / 2)
+pre_time = 8e-4
+gx_pre_epi = pp.make_trapezoid(channel='x', system=system_epi, area=-gx_epi.area / 2, duration=pre_time)
+gz_reph_epi = pp.make_trapezoid(channel='z', system=system_epi, area=-gz_epi.area / 2, duration=pre_time)
+gy_pre_epi = pp.make_trapezoid(channel='y', system=system_epi, area=-n_y_epi / 2 * dky, duration=pre_time)
+gy_blip_dur = np.ceil(2 * np.sqrt(dky / system_epi.max_slew) / 10e-6) * 10e-6
+gy_epi = pp.make_trapezoid(channel='y', system=system_epi, area=dky, duration=gy_blip_dur)
+
+for i_slice in range(n_slices):
+ rf_epi.freq_offset = gz_epi.amplitude * slice_thickness * (i_slice - (n_slices - 1) / 2)
+ seq_epi.add_block(rf_epi, gz_epi)
+ seq_epi.add_block(gx_pre_epi, gy_pre_epi, gz_reph_epi)
+ for _ in range(n_y_epi):
+ seq_epi.add_block(gx_epi, adc_epi)
+ seq_epi.add_block(gy_epi)
+ gx_epi.amplitude = -gx_epi.amplitude
+
+print("2/3 EPI sequence β opening viewer ...")
+seq_epi.plot(show_blocks=True) # show block boundaries
+
+
+# ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+# 3. Radial GRE (from pypulseq/examples/scripts/write_radial_gre.py)
+# ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+fov_rad = 260e-3; n_x_rad = 64; flip_rad = 10; n_spokes = 60; n_dummy = 20
+tr_rad = 20e-3; te_rad = 8e-3
+
+system_rad = pp.Opts(max_grad=28, grad_unit='mT/m', max_slew=120, slew_unit='T/m/s',
+ rf_ringdown_time=20e-6, rf_dead_time=100e-6, adc_dead_time=10e-6)
+seq_rad = pp.Sequence(system_rad)
+
+rf_rad, gz_rad, _ = pp.make_sinc_pulse(apodization=0.5, duration=4e-3,
+ flip_angle=np.deg2rad(flip_rad),
+ slice_thickness=slice_thickness, system=system_rad,
+ time_bw_product=4, return_gz=True,
+ delay=system_rad.rf_dead_time, use='excitation')
+
+dkx_rad = 1 / fov_rad
+gx_rad = pp.make_trapezoid(channel='x', flat_area=n_x_rad * dkx_rad, flat_time=6.4e-3 / 5, system=system_rad)
+adc_rad = pp.make_adc(num_samples=n_x_rad, duration=gx_rad.flat_time, delay=gx_rad.rise_time, system=system_rad)
+gx_pre_rad = pp.make_trapezoid(channel='x', area=-gx_rad.area / 2 - dkx_rad / 2, duration=2e-3, system=system_rad)
+gz_reph_rad = pp.make_trapezoid(channel='z', area=-gz_rad.area / 2, duration=2e-3, system=system_rad)
+gx_spoil_rad = pp.make_trapezoid(channel='x', area=0.5 * n_x_rad * dkx_rad, system=system_rad)
+gz_spoil_rad = pp.make_trapezoid(channel='z', area=4 / slice_thickness, system=system_rad)
+
+te_delay_rad = te_rad - pp.calc_duration(gx_pre_rad) - gz_rad.fall_time - gz_rad.flat_time / 2 - pp.calc_duration(gx_rad) / 2
+te_delay_rad = np.ceil(te_delay_rad / seq_rad.grad_raster_time) * seq_rad.grad_raster_time
+tr_delay_rad = tr_rad - pp.calc_duration(gx_pre_rad) - pp.calc_duration(gz_rad) - pp.calc_duration(gx_rad) - te_delay_rad
+tr_delay_rad = np.ceil(tr_delay_rad / seq_rad.grad_raster_time) * seq_rad.grad_raster_time
+
+spoke_angle = np.pi / n_spokes
+rf_phase = 0; rf_inc = 0
+for i_spoke in range(-n_dummy, n_spokes + 1):
+ rf_rad.phase_offset = rf_phase / 180 * np.pi
+ adc_rad.phase_offset = rf_phase / 180 * np.pi
+ rf_inc = divmod(rf_inc + 117, 360.0)[1]
+ rf_phase = divmod(rf_inc + rf_phase, 360.0)[1]
+ seq_rad.add_block(rf_rad, gz_rad)
+ phi = spoke_angle * (i_spoke - 1)
+ seq_rad.add_block(*pp.rotate(gx_pre_rad, gz_reph_rad, angle=phi, axis='z'))
+ seq_rad.add_block(pp.make_delay(te_delay_rad))
+ if i_spoke > 0:
+ seq_rad.add_block(*pp.rotate(gx_rad, adc_rad, angle=phi, axis='z'))
+ else:
+ seq_rad.add_block(*pp.rotate(gx_rad, angle=phi, axis='z'))
+ seq_rad.add_block(*pp.rotate(gx_spoil_rad, gz_spoil_rad, pp.make_delay(tr_delay_rad), angle=phi, axis='z'))
+
+print("3/3 Radial GRE sequence β opening viewer ...")
+seq_rad.plot(theme="dark") # overrides 'system' default from set()
+
+print("\nDone! 3 viewer tabs opened.")
diff --git a/python/pyproject.toml b/python/pyproject.toml
new file mode 100644
index 0000000..4de487a
--- /dev/null
+++ b/python/pyproject.toml
@@ -0,0 +1,47 @@
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = "seqeyes-python"
+version = "0.2.0"
+description = "Interactive Pulseq MRI sequence viewer for Jupyter β lightweight, beautiful, replaces seq.plot() in pypulseq."
+readme = "README.md"
+license = {text = "MIT"}
+authors = [
+ {name = "SeqEyes Developers"},
+]
+keywords = ["pulseq", "mri", "sequence", "visualization", "k-space", "gradient", "seqeyes"]
+classifiers = [
+ "Development Status :: 3 - Alpha",
+ "Intended Audience :: Science/Research",
+ "License :: OSI Approved :: MIT License",
+ "Programming Language :: Python :: 3",
+ "Programming Language :: Python :: 3.9",
+ "Programming Language :: Python :: 3.10",
+ "Programming Language :: Python :: 3.11",
+ "Programming Language :: Python :: 3.12",
+ "Topic :: Scientific/Engineering :: Medical Science Apps.",
+ "Topic :: Scientific/Engineering :: Visualization",
+ "Framework :: Jupyter",
+]
+requires-python = ">=3.9"
+dependencies = [
+ "numpy>=1.21",
+]
+
+[project.optional-dependencies]
+pypulseq = ["pypulseq>=1.4"]
+desktop = ["pywebview>=5"]
+test = ["pytest>=7"]
+
+[project.urls]
+Homepage = "https://github.com/bughht/seqeyes_plugin"
+Repository = "https://github.com/bughht/seqeyes_plugin"
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/seqeyes"]
+artifacts = ["src/seqeyes/resources/*"]
+
+[tool.hatch.build.targets.sdist]
+artifacts = ["src/seqeyes/resources/*"]
diff --git a/python/src/seqeyes/__init__.py b/python/src/seqeyes/__init__.py
new file mode 100644
index 0000000..1cde2fe
--- /dev/null
+++ b/python/src/seqeyes/__init__.py
@@ -0,0 +1,32 @@
+"""
+SeqEyes β Interactive Pulseq MRI Sequence Viewer
+=================================================
+
+Dropβin replacement for ``pypulseq.Sequence.plot()``.
+
+Quick start
+-----------
+>>> import seqeyes
+>>> seqeyes.set(theme="dark") # enable SeqEyes + set defaults
+
+>>> seq.plot() # interactive viewer in Jupyter
+>>> seq.plot(show_blocks=True) # perβcall overrides
+>>> seq.plot(time_range=(0, 0.05)) # zoom to first 50 ms
+
+>>> seqeyes.reset() # back to matplotlib
+
+In a plain ``.py`` script (no Jupyter):
+>>> seq.plot() # opens the viewer in your default web browser
+
+API
+---
+- :func:`set` β enable SeqEyes + configure defaults
+- :func:`reset` β restore matplotlib
+- :class:`SeqEyesViewer` β lowβlevel viewer
+"""
+
+from seqeyes._plot import set, reset
+from seqeyes._renderer import SeqEyesViewer
+
+__all__ = ["set", "reset", "SeqEyesViewer"]
+__version__ = "0.2.0"
diff --git a/python/src/seqeyes/_plot.py b/python/src/seqeyes/_plot.py
new file mode 100644
index 0000000..434aef8
--- /dev/null
+++ b/python/src/seqeyes/_plot.py
@@ -0,0 +1,221 @@
+"""
+_plot.py β SeqEyes plotting API.
+
+Call ``seqeyes.set()`` to switch ``seq.plot()`` to the interactive
+SeqEyes viewer. Call ``seqeyes.reset()`` to restore matplotlib.
+
+In Jupyter the viewer renders inline. In a ``.py`` script it opens
+a native desktop popβup window (requires ``pywebview``: ``pip install pywebview``).
+"""
+
+from __future__ import annotations
+
+import os
+import tempfile
+from pathlib import Path
+from typing import Any, Optional, Tuple, Union
+
+
+# ββ Module state ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+_defaults: dict[str, Any] = {}
+_original_plot = None
+_patched = False
+
+
+# ββ Public API ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+def set(
+ *,
+ theme: Optional[str] = None,
+ show_blocks: Optional[bool] = None,
+ time_disp: Optional[str] = None,
+ grad_disp: Optional[str] = None,
+ time_range: Optional[Tuple[float, float]] = None,
+) -> None:
+ """Enable SeqEyes for all ``seq.plot()`` calls and set global defaults.
+
+ Call this once at the top of your notebook / script β with or
+ without arguments. Any keywords become defaults applied to every
+ subsequent ``seq.plot()`` (overridable perβcall).
+
+ >>> seqeyes.set() # enable, system defaults
+ >>> seqeyes.set(theme="dark", time_disp="ms") # enable + defaults
+ >>> seq.plot() # uses dark + ms
+ >>> seq.plot(theme="light") # overrides theme only
+ """
+ _store_defaults(theme, show_blocks, time_disp, grad_disp, time_range)
+ _ensure_patched()
+
+
+def reset() -> None:
+ """Restore matplotlib ``seq.plot()`` and clear all SeqEyes defaults.
+
+ Call ``seqeyes.set()`` again to reβenable at any time.
+ """
+ global _patched
+ _defaults.clear()
+ _patched = False
+
+ try:
+ from pypulseq import Sequence as _Seq
+ except ImportError:
+ return
+
+ if _original_plot is not None:
+ _Seq.plot = _original_plot # type: ignore[attr-defined]
+
+
+# ββ Internal helpers ββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+def _store_defaults(
+ theme: Optional[str],
+ show_blocks: Optional[bool],
+ time_disp: Optional[str],
+ grad_disp: Optional[str],
+ time_range: Optional[Tuple[float, float]],
+) -> None:
+ for name, val in [
+ ("theme", theme),
+ ("show_blocks", show_blocks),
+ ("time_disp", time_disp),
+ ("grad_disp", grad_disp),
+ ("time_range", time_range),
+ ]:
+ if val is not None:
+ _defaults[name] = val
+
+
+def _ensure_patched() -> None:
+ global _patched, _original_plot
+
+ if _patched:
+ return
+ _patched = True
+
+ try:
+ from pypulseq import Sequence as _Seq
+ except ImportError:
+ raise ImportError("pypulseq is not installed. Install with: pip install pypulseq")
+
+ if _original_plot is None:
+ _original_plot = getattr(_Seq, "plot", None)
+
+ # ββ Replacement .plot() ββββββββββββββββββββββββββββββββββββββββββ
+ def _seqeyes_plot(
+ self: object,
+ label: str = "",
+ show_blocks: bool = False,
+ time_range: Any = (0, float("inf")),
+ time_disp: str = "s",
+ grad_disp: str = "kHz/m",
+ **kwargs: Any,
+ ) -> None:
+ # Merge globals β perβcall args take precedence
+ sb = _defaults.get("show_blocks", show_blocks)
+ td = _defaults.get("time_disp", time_disp)
+ gd = _defaults.get("grad_disp", grad_disp)
+ tr = _defaults.get("time_range", time_range)
+ th = str(kwargs.pop("theme", _defaults.get("theme", "system")))
+
+ seq_text = _seq_to_text(self)
+ if not seq_text:
+ raise RuntimeError("Could not write .seq text from the pypulseq Sequence.")
+
+ from seqeyes._renderer import SeqEyesViewer
+
+ viewer = SeqEyesViewer(
+ seq_text, label=label, show_blocks=sb,
+ time_range=tr, time_disp=td, grad_disp=gd, theme=th,
+ )
+
+ if _in_jupyter():
+ from IPython.display import display
+ display(viewer)
+ else:
+ _open_in_browser(viewer)
+
+ # Preserve original as _plot_matplotlib
+ if _original_plot is not None and not hasattr(_Seq, "_plot_matplotlib"):
+ _Seq._plot_matplotlib = _original_plot # type: ignore[attr-defined]
+
+ _Seq.plot = _seqeyes_plot # type: ignore[attr-defined]
+
+ # _repr_html_ so bare ``seq`` autoβrenders in Jupyter
+ def _seq_repr_html_(self: object) -> str:
+ from seqeyes._renderer import SeqEyesViewer
+ return SeqEyesViewer(_seq_to_text(self))._repr_html_()
+
+ _Seq._repr_html_ = _seq_repr_html_ # type: ignore[attr-defined]
+
+
+def _in_jupyter() -> bool:
+ try:
+ from IPython import get_ipython
+ s = get_ipython()
+ return s is not None and ("ZMQ" in s.__class__.__name__ or "Shell" in s.__class__.__name__)
+ except ImportError:
+ return False
+
+
+def _open_in_browser(viewer: "SeqEyesViewer") -> None: # noqa: F821
+ """Open the viewer in a native desktop popβup window (pywebview).
+
+ Spawns a subprocess so ``webview.start()`` runs on its own main
+ thread without blocking the calling script. Falls back to the
+ system browser if pywebview is not installed.
+ """
+ html = viewer.to_html()
+
+ # Write HTML to a temp file β file:// URLs are more reliable with
+ # pywebview than inline html= on some backends.
+ p = os.path.join(tempfile.gettempdir(), "seqeyes_viewer.html")
+ Path(p).write_text(html, encoding="utf-8")
+ url = f"file:///{p.replace(os.sep, '/')}"
+
+ try:
+ import webview # type: ignore[import-untyped]
+ except ImportError:
+ _fallback_browser(url)
+ return
+
+ # Run webview in a subprocess so it gets its own main thread.
+ # We shell out to a tiny inline script β this avoids multiprocessing
+ # "spawn" issues (reβimporting the parent script on Windows).
+ import subprocess
+ import sys
+
+ script = (
+ "import webview;"
+ f"webview.create_window(title='SeqEyes β Pulseq MRI Sequence Viewer',url={url!r},width=1280,height=800,resizable=True,easy_drag=False);"
+ "webview.start()"
+ )
+ try:
+ subprocess.Popen(
+ [sys.executable, "-c", script],
+ stdout=subprocess.DEVNULL,
+ stderr=subprocess.DEVNULL,
+ )
+ except Exception:
+ _fallback_browser(url)
+
+
+def _fallback_browser(url: str) -> None:
+ """Last resort: open the viewer in the system web browser."""
+ import webbrowser
+ webbrowser.open(url)
+
+
+def _seq_to_text(seq: object) -> str:
+ if hasattr(seq, "write"):
+ with tempfile.NamedTemporaryFile(mode="r", suffix=".seq", delete=False, encoding="utf-8") as f:
+ tp = f.name
+ try:
+ seq.write(tp) # type: ignore[union-attr]
+ return Path(tp).read_text(encoding="utf-8")
+ finally:
+ try:
+ os.unlink(tp)
+ except OSError:
+ pass
+ return ""
diff --git a/python/src/seqeyes/_renderer.py b/python/src/seqeyes/_renderer.py
new file mode 100644
index 0000000..7400960
--- /dev/null
+++ b/python/src/seqeyes/_renderer.py
@@ -0,0 +1,236 @@
+"""
+_renderer.py β Jupyter HTML renderer for the SeqEyes interactive viewer.
+
+Generates a selfβcontained HTML string that Jupyter displays inline via
+``_repr_html_()``, exactly the same mechanism Plotly uses. The viewer is
+an iframe embedding the SeqEyes web UI with sequence data injected.
+"""
+
+from __future__ import annotations
+
+import base64
+import json
+import os
+from pathlib import Path
+from typing import Optional
+
+
+# ββ Paths βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
+_RESOURCES_DIR = Path(__file__).resolve().parent / "resources"
+_VIEWER_TEMPLATE_PATH = _RESOURCES_DIR / "viewer.html"
+
+# pulseq-bundle.js lives in the repo root; during dev we read it directly.
+# When installed as a package, a copy should be placed in resources/.
+_BUNDLE_CANDIDATES = [
+ _RESOURCES_DIR / "pulseq-bundle.js",
+ Path(__file__).resolve().parent.parent.parent.parent / "web" / "pulseq-bundle.js",
+ Path.cwd() / "web" / "pulseq-bundle.js",
+]
+
+
+def _find_bundle() -> str:
+ """Locate and read the pulseq-bundle.js file."""
+ for p in _BUNDLE_CANDIDATES:
+ if p.is_file():
+ return p.read_text(encoding="utf-8")
+ raise FileNotFoundError(
+ "pulseq-bundle.js not found. "
+ "Build it with: npm run build:web (from the repo root), "
+ "or copy it to python/src/seqeyes/resources/"
+ )
+
+
+def _read_viewer_template() -> str:
+ """Read the standalone viewer HTML template."""
+ if _VIEWER_TEMPLATE_PATH.is_file():
+ return _VIEWER_TEMPLATE_PATH.read_text(encoding="utf-8")
+ alt = Path.cwd() / "python" / "src" / "seqeyes" / "resources" / "viewer.html"
+ if alt.is_file():
+ return alt.read_text(encoding="utf-8")
+ raise FileNotFoundError(
+ f"Viewer template not found at {_VIEWER_TEMPLATE_PATH}. "
+ "Ensure the package is installed correctly."
+ )
+
+
+def _build_html(
+ seq_text: str,
+ *,
+ theme: Optional[str] = None,
+ inject_bundle: bool = True,
+ label: str = "",
+ show_blocks: bool = False,
+ time_range: tuple = (0, float("inf")),
+ time_disp: str = "s",
+ grad_disp: str = "kHz/m",
+) -> str:
+ """Assemble the complete viewer HTML.
+
+ Parameters
+ ----------
+ seq_text : str
+ Raw .seq file content.
+ theme : str or None
+ CSS theme class to apply to ````.
+ inject_bundle : bool
+ If True, inline the pulseq-bundle.js parser into the HTML.
+ label : str
+ Display label (injected as JS variable, not yet rendered).
+ show_blocks : bool
+ Initial blockβboundary visibility.
+ time_range : tuple[float, float]
+ Initial viewport range in seconds.
+ time_disp : str
+ Initial time display unit (``"s"``, ``"ms"``, ``"us"``).
+ grad_disp : str
+ Initial gradient display unit (``"Hz/m"``, ``"kHz/m"``,
+ ``"mT/m"``, ``"G/cm"``).
+ """
+ template = _read_viewer_template()
+ seq_b64 = base64.b64encode(seq_text.encode("utf-8")).decode("ascii")
+
+ # 1. Inject the pulseq-bundle.js
+ bundle_placeholder = ""
+ if inject_bundle and bundle_placeholder in template:
+ bundle_js = _find_bundle()
+ template = template.replace(
+ bundle_placeholder,
+ f"",
+ )
+ elif bundle_placeholder in template:
+ template = template.replace(bundle_placeholder, "")
+
+ # 2. Build options injection block (sequence data + display opts)
+ opts = [
+ f"window.SEQEYES_RAW_B64 = {json.dumps(seq_b64)};",
+ f"window.SEQEYES_LABEL = {json.dumps(label)};",
+ f"window.SEQEYES_SHOW_BLOCKS = {json.dumps(show_blocks)};",
+ f"window.SEQEYES_TIME_RANGE = {json.dumps(list(time_range))};",
+ f"window.SEQEYES_TIME_DISP = {json.dumps(time_disp)};",
+ f"window.SEQEYES_GRAD_DISP = {json.dumps(grad_disp)};",
+ ]
+ opts_block = "\n".join(opts)
+
+ data_placeholder = "/* SEQEYES_DATA_PLACEHOLDER */"
+ if data_placeholder in template:
+ template = template.replace(data_placeholder, opts_block)
+ else:
+ template = template.replace(
+ "",
+ f"\n
",
+ f'',
+ )
+
+ return template
+
+
+class SeqEyesViewer:
+ """Interactive Pulseq sequence viewer for Jupyter notebook output.
+
+ Normally you don't create this directly β call ``seqeyes.set()``
+ once, then use ``seq.plot(...)`` on any ``pypulseq.Sequence``.
+
+ Parameters
+ ----------
+ seq_text : str
+ Raw .seq file content as a string.
+ label : str
+ Display label (not yet rendered on the viewer).
+ show_blocks : bool
+ Whether to show blockβboundary lines initially.
+ time_range : tuple[float, float]
+ Initial time viewport ``(start_sec, end_sec)``. ``(0, inf)``
+ shows the whole sequence.
+ time_disp : str
+ Time axis unit β ``"s"``, ``"ms"``, or ``"us"``.
+ grad_disp : str
+ Gradient axis unit β ``"Hz/m"``, ``"kHz/m"``, ``"mT/m"``, ``"G/cm"``.
+ theme : str or None
+ Colour theme.
+ width : str
+ CSS width of the iframe (e.g. ``"100%"``).
+ height : str
+ CSS height of the iframe (e.g. ``"550px"``).
+ """
+
+ def __init__(
+ self,
+ seq_text: str,
+ *,
+ label: str = "",
+ show_blocks: bool = False,
+ time_range: tuple = (0, float("inf")),
+ time_disp: str = "s",
+ grad_disp: str = "kHz/m",
+ theme: Optional[str] = None,
+ width: str = "100%",
+ height: str = "550px",
+ ) -> None:
+ self._seq_text = seq_text
+ self._label = label
+ self._show_blocks = show_blocks
+ self._time_range = time_range
+ self._time_disp = time_disp
+ self._grad_disp = grad_disp
+ self._theme = theme
+ self._width = width
+ self._height = height
+
+ # ββ Jupyter integration βββββββββββββββββββββββββββββββββββββββββββ
+
+ def _repr_html_(self) -> str:
+ """Return an HTML iframe that Jupyter renders inline."""
+ html = _build_html(
+ self._seq_text,
+ theme=self._theme,
+ label=self._label,
+ show_blocks=self._show_blocks,
+ time_range=self._time_range,
+ time_disp=self._time_disp,
+ grad_disp=self._grad_disp,
+ )
+ b64 = base64.b64encode(html.encode("utf-8")).decode("ascii")
+ return (
+ f'"
+ )
+
+ def _ipython_display_(self) -> None:
+ """IPython display hook."""
+ from IPython.display import HTML, display
+
+ display(HTML(self._repr_html_()))
+
+ # ββ Convenience ββββββββββββββββββββββββββββββββββββββββββββββββββββ
+
+ def show(self) -> "SeqEyesViewer":
+ """Display the viewer (useful in IPython when autoβdisplay is off)."""
+ from IPython.display import display
+
+ display(self)
+ return self
+
+ # ββ Direct HTML access (for debugging / custom embedding) ββββββββββ
+
+ def to_html(self, *, inject_bundle: bool = True) -> str:
+ """Return the full standalone HTML string (for saving to a file)."""
+ return _build_html(
+ self._seq_text,
+ theme=self._theme,
+ inject_bundle=inject_bundle,
+ label=self._label,
+ show_blocks=self._show_blocks,
+ time_range=self._time_range,
+ time_disp=self._time_disp,
+ grad_disp=self._grad_disp,
+ )
diff --git a/python/src/seqeyes/resources/pulseq-bundle.js b/python/src/seqeyes/resources/pulseq-bundle.js
new file mode 100644
index 0000000..b94ba49
--- /dev/null
+++ b/python/src/seqeyes/resources/pulseq-bundle.js
@@ -0,0 +1,2296 @@
+"use strict";
+var Pulseq = (() => {
+ var __defProp = Object.defineProperty;
+ var __getOwnPropDesc = Object.getOwnPropertyDescriptor;
+ var __getOwnPropNames = Object.getOwnPropertyNames;
+ var __hasOwnProp = Object.prototype.hasOwnProperty;
+ var __export = (target, all) => {
+ for (var name in all)
+ __defProp(target, name, { get: all[name], enumerable: true });
+ };
+ var __copyProps = (to, from, except, desc) => {
+ if (from && typeof from === "object" || typeof from === "function") {
+ for (let key of __getOwnPropNames(from))
+ if (!__hasOwnProp.call(to, key) && key !== except)
+ __defProp(to, key, { get: () => from[key], enumerable: !(desc = __getOwnPropDesc(from, key)) || desc.enumerable });
+ }
+ return to;
+ };
+ var __toCommonJS = (mod) => __copyProps(__defProp({}, "__esModule", { value: true }), mod);
+
+ // web/pulseq-browser.ts
+ var pulseq_browser_exports = {};
+ __export(pulseq_browser_exports, {
+ PACKAGE_VERSION: () => PACKAGE_VERSION,
+ calculateKspace: () => calculateKspace,
+ calculateM1: () => calculateM1,
+ calculatePns: () => calculatePns,
+ decodeAllBlocks: () => decodeAllBlocks,
+ detectSequenceTiming: () => detectSequenceTiming,
+ exportKspaceArtifacts: () => exportKspaceArtifacts,
+ formatTrajectoryText: () => formatTrajectoryText,
+ getTotalDuration: () => getTotalDuration,
+ parsePnsHardwareAsc: () => parsePnsHardwareAsc,
+ parseSequenceText: () => parseSequenceText,
+ safePnsModel: () => safePnsModel
+ });
+
+ // package.json
+ var version = "0.1.17";
+
+ // src/pulseq/decompressor.ts
+ function decompressShape(compressed, numSamples) {
+ const packedLen = compressed.length;
+ if (!Number.isInteger(numSamples) || numSamples <= 0) {
+ throw new Error(`Invalid shape sample count: ${numSamples}`);
+ }
+ if (packedLen === numSamples) {
+ return new Float64Array(compressed);
+ }
+ const result = new Float64Array(numSamples);
+ let iPacked = 0;
+ let iUnpacked = 0;
+ while (iPacked < packedLen && iUnpacked < numSamples) {
+ if (iPacked + 1 >= packedLen) {
+ result[iUnpacked] = compressed[iPacked];
+ iPacked++;
+ iUnpacked++;
+ break;
+ }
+ if (compressed[iPacked] !== compressed[iPacked + 1]) {
+ result[iUnpacked] = compressed[iPacked];
+ iPacked++;
+ iUnpacked++;
+ } else {
+ if (iPacked + 2 >= packedLen) {
+ throw new Error("Malformed compressed shape: repeat marker is missing its count");
+ }
+ const value = compressed[iPacked];
+ const rawRepeat = compressed[iPacked + 2];
+ const repeatCount = Math.round(rawRepeat) + 2;
+ if (Math.abs(rawRepeat + 2 - repeatCount) > 1e-6 || repeatCount < 2) {
+ throw new Error(`Malformed compressed shape: invalid repeat count ${rawRepeat}`);
+ }
+ if (iUnpacked + repeatCount > numSamples) {
+ throw new Error("Malformed compressed shape: repeat block exceeds expected sample count");
+ }
+ iPacked += 3;
+ const end = iUnpacked + repeatCount;
+ while (iUnpacked < end) {
+ result[iUnpacked] = value;
+ iUnpacked++;
+ }
+ }
+ }
+ if (iUnpacked !== numSamples) {
+ throw new Error(`Malformed compressed shape: expected ${numSamples} samples, decoded ${iUnpacked}`);
+ }
+ let cumSum = 0;
+ for (let i = 0; i < numSamples; i++) {
+ cumSum += result[i];
+ result[i] = cumSum;
+ }
+ return result;
+ }
+
+ // src/pulseq/types.ts
+ var VER_PRE_14 = 1004e3;
+ var VER_V15 = 1005e3;
+ var VER_V15001 = 1005001;
+ function makeVersionCombined(major, minor, revision) {
+ return major * 1e6 + minor * 1e3 + revision;
+ }
+
+ // src/pulseq/reader.ts
+ function parseSequenceText(text) {
+ const lines = text.split(/\r?\n/);
+ const seq = createEmptySequence();
+ const seenSections = /* @__PURE__ */ new Set();
+ let sectionName = null;
+ let sectionLines = [];
+ for (const line of lines) {
+ const m = line.match(/^\[(\w+)\]$/);
+ if (m) {
+ if (sectionName) dispatchSection(seq, sectionName, sectionLines);
+ sectionName = m[1];
+ seenSections.add(sectionName);
+ sectionLines = [];
+ } else {
+ sectionLines.push(line);
+ }
+ }
+ if (sectionName) dispatchSection(seq, sectionName, sectionLines);
+ seq.versionCombined = makeVersionCombined(
+ seq.version.major,
+ seq.version.minor,
+ seq.version.revision
+ );
+ extractRasterTimes(seq);
+ validateSequence(seq, seenSections);
+ return seq;
+ }
+ function dispatchSection(seq, name, lines) {
+ const valid = lines.filter((l) => {
+ const t = l.trim();
+ return t && !t.startsWith("#");
+ });
+ switch (name) {
+ case "VERSION":
+ parseVersion(seq, valid);
+ break;
+ case "DEFINITIONS":
+ parseDefinitions(seq, valid);
+ break;
+ case "BLOCKS":
+ parseBlocks(seq, valid);
+ break;
+ case "RF":
+ parseRF(seq, valid);
+ break;
+ case "GRADIENTS":
+ parseArbitraryGrads(seq, valid);
+ break;
+ case "TRAP":
+ parseTrapGrads(seq, valid);
+ break;
+ case "ADC":
+ parseADC(seq, valid);
+ break;
+ case "EXTENSIONS":
+ parseExtensions(seq, valid);
+ break;
+ case "SHAPES":
+ parseShapes(seq, lines);
+ break;
+ }
+ }
+ function createEmptySequence() {
+ return {
+ version: { major: 1, minor: 0, revision: 0 },
+ versionCombined: 0,
+ definitions: /* @__PURE__ */ new Map(),
+ definitionsRaw: /* @__PURE__ */ new Map(),
+ blocks: [],
+ rfs: /* @__PURE__ */ new Map(),
+ arbitraryGrads: /* @__PURE__ */ new Map(),
+ trapGrads: /* @__PURE__ */ new Map(),
+ adcs: /* @__PURE__ */ new Map(),
+ extensions: /* @__PURE__ */ new Map(),
+ extensionNames: /* @__PURE__ */ new Map(),
+ extensionTypes: /* @__PURE__ */ new Map(),
+ triggers: [],
+ ncos: [],
+ rotations: [],
+ labelSets: [],
+ labelIncs: [],
+ softDelays: [],
+ rfShims: [],
+ shapes: /* @__PURE__ */ new Map(),
+ rasterTimes: { blockDurationRaster: 1e-5, gradientRaster: 1e-5, rfRaster: 1e-6, adcRaster: 1e-7 }
+ };
+ }
+ function ver(seq) {
+ if (seq.versionCombined > 0) return seq.versionCombined;
+ return makeVersionCombined(seq.version.major, seq.version.minor, seq.version.revision);
+ }
+ function parseError(message) {
+ throw new Error(`Pulseq parse error: ${message}`);
+ }
+ function requireFieldCount(section, line, count, allowed) {
+ const allowedCounts = Array.isArray(allowed) ? allowed : [allowed];
+ if (!allowedCounts.includes(count)) {
+ parseError(`${section} row has ${count} fields, expected ${allowedCounts.join(" or ")}: ${line}`);
+ }
+ }
+ function toNumber(value, section, line) {
+ const n = Number(value);
+ if (!Number.isFinite(n)) parseError(`${section} row contains a non-numeric field '${value}': ${line}`);
+ return n;
+ }
+ function toInt(value, section, line) {
+ const n = toNumber(value, section, line);
+ if (!Number.isInteger(n)) parseError(`${section} row contains a non-integer field '${value}': ${line}`);
+ return n;
+ }
+ function splitFields(line) {
+ return line.trim().split(/\s+/);
+ }
+ function parseVersion(seq, lines) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ requireFieldCount("VERSION", line, p.length, 2);
+ const [k, v] = p;
+ const n = toInt(v, "VERSION", line);
+ if (k === "major") seq.version.major = n;
+ else if (k === "minor") seq.version.minor = n;
+ else if (k === "revision") seq.version.revision = n;
+ }
+ seq.versionCombined = makeVersionCombined(
+ seq.version.major,
+ seq.version.minor,
+ seq.version.revision
+ );
+ }
+ function parseDefinitions(seq, lines) {
+ for (const line of lines) {
+ const idx = line.search(/\s/);
+ if (idx < 0) {
+ seq.definitions.set(line.trim(), []);
+ continue;
+ }
+ const key = line.substring(0, idx);
+ const vals = line.substring(idx + 1).trim().split(/\s+/).map(Number).filter((n) => !isNaN(n));
+ seq.definitions.set(key, vals);
+ seq.definitionsRaw.set(key, line.substring(idx + 1).trim());
+ }
+ }
+ function parseBlocks(seq, lines) {
+ const vc = ver(seq);
+ for (const line of lines) {
+ const p = splitFields(line);
+ requireFieldCount("BLOCKS", line, p.length, [7, 8]);
+ const num = toInt(p[0], "BLOCKS", line);
+ const extId = p.length === 8 ? toInt(p[7], "BLOCKS", line) : 0;
+ if (vc < VER_PRE_14) {
+ seq.blocks.push({
+ num,
+ dur: toNumber(p[1], "BLOCKS", line),
+ rfId: toInt(p[2], "BLOCKS", line),
+ gxId: toInt(p[3], "BLOCKS", line),
+ gyId: toInt(p[4], "BLOCKS", line),
+ gzId: toInt(p[5], "BLOCKS", line),
+ adcId: toInt(p[6], "BLOCKS", line),
+ extId
+ });
+ } else {
+ seq.blocks.push({
+ num,
+ dur: toNumber(p[1], "BLOCKS", line),
+ rfId: toInt(p[2], "BLOCKS", line),
+ gxId: toInt(p[3], "BLOCKS", line),
+ gyId: toInt(p[4], "BLOCKS", line),
+ gzId: toInt(p[5], "BLOCKS", line),
+ adcId: toInt(p[6], "BLOCKS", line),
+ extId
+ });
+ }
+ }
+ }
+ function parseRF(seq, lines) {
+ const vc = ver(seq);
+ for (const line of lines) {
+ const parts = splitFields(line);
+ const id = toInt(parts[0], "RF", line);
+ const amp = toNumber(parts[1], "RF", line);
+ const magId = toInt(parts[2], "RF", line);
+ const phId = toInt(parts[3], "RF", line);
+ if (vc >= VER_V15) {
+ requireFieldCount("RF", line, parts.length, 12);
+ const use = parts[11].toLowerCase();
+ if (!/^[erisu]$/.test(use)) parseError(`RF row has invalid use flag '${parts[11]}': ${line}`);
+ seq.rfs.set(id, {
+ id,
+ amplitude: amp,
+ magShapeId: magId,
+ phaseShapeId: phId,
+ timeShapeId: toInt(parts[4], "RF", line),
+ center: toNumber(parts[5], "RF", line),
+ delay: toNumber(parts[6], "RF", line),
+ freqPPM: toNumber(parts[7], "RF", line),
+ phasePPM: toNumber(parts[8], "RF", line),
+ freqOffset: toNumber(parts[9], "RF", line),
+ phaseOffset: toNumber(parts[10], "RF", line),
+ phaseModShapeId: 0,
+ use
+ });
+ } else if (vc >= VER_PRE_14) {
+ requireFieldCount("RF", line, parts.length, 8);
+ seq.rfs.set(id, {
+ id,
+ amplitude: amp,
+ magShapeId: magId,
+ phaseShapeId: phId,
+ timeShapeId: toInt(parts[4], "RF", line),
+ center: -1,
+ // not in v1.4.x
+ delay: toNumber(parts[5], "RF", line),
+ freqPPM: 0,
+ phasePPM: 0,
+ freqOffset: toNumber(parts[6], "RF", line),
+ phaseOffset: toNumber(parts[7], "RF", line),
+ phaseModShapeId: 0,
+ use: "u"
+ });
+ } else {
+ requireFieldCount("RF", line, parts.length, 7);
+ seq.rfs.set(id, {
+ id,
+ amplitude: amp,
+ magShapeId: magId,
+ phaseShapeId: phId,
+ timeShapeId: 0,
+ center: -1,
+ delay: toNumber(parts[4], "RF", line),
+ freqPPM: 0,
+ phasePPM: 0,
+ freqOffset: toNumber(parts[5], "RF", line),
+ phaseOffset: toNumber(parts[6], "RF", line),
+ phaseModShapeId: 0,
+ use: "u"
+ });
+ }
+ }
+ }
+ function parseArbitraryGrads(seq, lines) {
+ const vc = ver(seq);
+ for (const line of lines) {
+ const p = splitFields(line);
+ const id = toInt(p[0], "GRADIENTS", line);
+ if (vc >= VER_V15) {
+ requireFieldCount("GRADIENTS", line, p.length, 7);
+ seq.arbitraryGrads.set(id, {
+ id,
+ amplitude: toNumber(p[1], "GRADIENTS", line),
+ first: toNumber(p[2], "GRADIENTS", line),
+ last: toNumber(p[3], "GRADIENTS", line),
+ shapeId: toInt(p[4], "GRADIENTS", line),
+ timeId: toInt(p[5], "GRADIENTS", line),
+ delay: toNumber(p[6], "GRADIENTS", line)
+ });
+ } else if (vc >= VER_PRE_14) {
+ requireFieldCount("GRADIENTS", line, p.length, 5);
+ seq.arbitraryGrads.set(id, {
+ id,
+ amplitude: toNumber(p[1], "GRADIENTS", line),
+ first: NaN,
+ last: NaN,
+ shapeId: toInt(p[2], "GRADIENTS", line),
+ timeId: toInt(p[3], "GRADIENTS", line),
+ delay: toNumber(p[4], "GRADIENTS", line)
+ });
+ } else {
+ requireFieldCount("GRADIENTS", line, p.length, 4);
+ seq.arbitraryGrads.set(id, {
+ id,
+ amplitude: toNumber(p[1], "GRADIENTS", line),
+ first: NaN,
+ last: NaN,
+ shapeId: toInt(p[2], "GRADIENTS", line),
+ timeId: 0,
+ delay: toNumber(p[3], "GRADIENTS", line)
+ });
+ }
+ }
+ }
+ function parseTrapGrads(seq, lines) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ requireFieldCount("TRAP", line, p.length, 6);
+ const id = toInt(p[0], "TRAP", line);
+ seq.trapGrads.set(id, {
+ id,
+ amplitude: toNumber(p[1], "TRAP", line),
+ rise: toNumber(p[2], "TRAP", line),
+ flat: toNumber(p[3], "TRAP", line),
+ fall: toNumber(p[4], "TRAP", line),
+ delay: toNumber(p[5], "TRAP", line)
+ });
+ }
+ }
+ function parseADC(seq, lines) {
+ const vc = ver(seq);
+ for (const line of lines) {
+ const p = splitFields(line);
+ const id = toInt(p[0], "ADC", line);
+ if (vc >= VER_V15) {
+ requireFieldCount("ADC", line, p.length, 9);
+ seq.adcs.set(id, {
+ id,
+ numSamples: toInt(p[1], "ADC", line),
+ dwell: toNumber(p[2], "ADC", line),
+ delay: toNumber(p[3], "ADC", line),
+ freqPPM: toNumber(p[4], "ADC", line),
+ phasePPM: toNumber(p[5], "ADC", line),
+ freqOffset: toNumber(p[6], "ADC", line),
+ phaseOffset: toNumber(p[7], "ADC", line),
+ deadTime: 0,
+ discardPre: 0,
+ discardPost: 0,
+ phaseModShapeId: toInt(p[8], "ADC", line)
+ });
+ } else {
+ requireFieldCount("ADC", line, p.length, 6);
+ seq.adcs.set(id, {
+ id,
+ numSamples: toInt(p[1], "ADC", line),
+ dwell: toNumber(p[2], "ADC", line),
+ delay: toNumber(p[3], "ADC", line),
+ freqPPM: 0,
+ phasePPM: 0,
+ freqOffset: toNumber(p[4], "ADC", line),
+ phaseOffset: toNumber(p[5], "ADC", line),
+ deadTime: 0,
+ discardPre: 0,
+ discardPost: 0,
+ phaseModShapeId: 0
+ });
+ }
+ }
+ }
+ function parseExtensions(seq, valid) {
+ const vc = ver(seq);
+ _unknownLabelCounter = 0;
+ _unknownLabels.clear();
+ let i = 0;
+ while (i < valid.length) {
+ const line = valid[i].trim();
+ if (line.startsWith("extension ")) break;
+ const p = splitFields(line);
+ requireFieldCount("EXTENSIONS", line, p.length, 4);
+ const id = toInt(p[0], "EXTENSIONS", line);
+ seq.extensions.set(id, {
+ id,
+ type: toInt(p[1], "EXTENSIONS", line),
+ ref: toInt(p[2], "EXTENSIONS", line),
+ nextId: toInt(p[3], "EXTENSIONS", line)
+ });
+ i++;
+ }
+ while (i < valid.length) {
+ const line = valid[i].trim();
+ const extM = line.match(/^extension\s+(\w+)\s+(\d+)/i);
+ if (!extM) {
+ i++;
+ continue;
+ }
+ const extName = extM[1].toUpperCase();
+ const extId = +extM[2];
+ seq.extensionNames.set(extId, extName);
+ seq.extensionTypes.set(extId, extensionNameToType(extName));
+ i++;
+ const dataLines = [];
+ while (i < valid.length && !valid[i].trim().startsWith("extension ")) {
+ dataLines.push(valid[i].trim());
+ i++;
+ }
+ switch (extName) {
+ case "TRIGGERS":
+ parseTriggerSpecs(seq, dataLines);
+ break;
+ case "NCO":
+ parseNCOSpecs(seq, dataLines);
+ break;
+ case "ROTATIONS":
+ parseRotationSpecs(seq, dataLines, vc);
+ break;
+ case "LABELSET":
+ parseLabelSpecs(seq, dataLines, true);
+ break;
+ case "LABELINC":
+ parseLabelSpecs(seq, dataLines, false);
+ break;
+ case "DELAYS":
+ parseSoftDelaySpecs(seq, dataLines);
+ break;
+ case "RF_SHIMS":
+ parseRFShimSpecs(seq, dataLines);
+ break;
+ default:
+ break;
+ }
+ }
+ }
+ function extensionNameToType(name) {
+ switch (name.toUpperCase()) {
+ case "TRIGGERS":
+ return 1 /* EXT_TRIGGER */;
+ case "ROTATIONS":
+ return 2 /* EXT_ROTATION */;
+ case "LABELSET":
+ return 3 /* EXT_LABELSET */;
+ case "LABELINC":
+ return 4 /* EXT_LABELINC */;
+ case "DELAYS":
+ return 5 /* EXT_DELAY */;
+ case "RF_SHIMS":
+ return 6 /* EXT_RF_SHIM */;
+ case "NCO":
+ return 100 /* EXT_NCO */;
+ default:
+ return 999 /* EXT_UNKNOWN */;
+ }
+ }
+ function parseTriggerSpecs(seq, lines) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ requireFieldCount("TRIGGERS", line, p.length, 5);
+ seq.triggers.push({
+ id: toInt(p[0], "TRIGGERS", line),
+ triggerType: toInt(p[1], "TRIGGERS", line),
+ channel: toInt(p[2], "TRIGGERS", line),
+ delay: toNumber(p[3], "TRIGGERS", line),
+ duration: toNumber(p[4], "TRIGGERS", line)
+ });
+ }
+ }
+ function parseNCOSpecs(seq, lines) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ requireFieldCount("NCO", line, p.length, 6);
+ seq.ncos.push({
+ id: toInt(p[0], "NCO", line),
+ channel: toInt(p[1], "NCO", line),
+ frequency: toNumber(p[2], "NCO", line),
+ phase: toNumber(p[3], "NCO", line),
+ delay: toNumber(p[4], "NCO", line),
+ duration: toNumber(p[5], "NCO", line)
+ });
+ }
+ }
+ function parseRotationSpecs(seq, lines, vc) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ if (vc >= VER_V15) {
+ requireFieldCount("ROTATIONS", line, p.length, 5);
+ const [q0, q1, q2, q3] = [
+ toNumber(p[1], "ROTATIONS", line),
+ toNumber(p[2], "ROTATIONS", line),
+ toNumber(p[3], "ROTATIONS", line),
+ toNumber(p[4], "ROTATIONS", line)
+ ];
+ const norm = Math.sqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3);
+ if (Math.abs(norm - 1) > 1e-3 || norm === 0) {
+ parseError(`ROTATIONS row has a non-normalized quaternion: ${line}`);
+ }
+ seq.rotations.push({
+ id: toInt(p[0], "ROTATIONS", line),
+ values: [q0 / norm, q1 / norm, q2 / norm, q3 / norm]
+ });
+ } else {
+ requireFieldCount("ROTATIONS", line, p.length, 10);
+ seq.rotations.push({
+ id: toInt(p[0], "ROTATIONS", line),
+ values: p.slice(1, 10).map((v) => toNumber(v, "ROTATIONS", line))
+ });
+ }
+ }
+ }
+ var KNOWN_LABELS = {
+ "SLC": { labelId: 0, flagId: 0 },
+ "SEG": { labelId: 1, flagId: 0 },
+ "REP": { labelId: 2, flagId: 0 },
+ "AVG": { labelId: 3, flagId: 0 },
+ "ECO": { labelId: 4, flagId: 0 },
+ "PHS": { labelId: 5, flagId: 0 },
+ "SET": { labelId: 6, flagId: 0 },
+ "ACQ": { labelId: 7, flagId: 0 },
+ "LIN": { labelId: 8, flagId: 0 },
+ "PAR": { labelId: 9, flagId: 0 },
+ "ONCE": { labelId: 10, flagId: 0 },
+ "NAV": { labelId: 0, flagId: 1 },
+ "REV": { labelId: 0, flagId: 2 },
+ "SMS": { labelId: 0, flagId: 4 },
+ "REF": { labelId: 0, flagId: 8 },
+ "IMA": { labelId: 0, flagId: 16 },
+ "OFF": { labelId: 0, flagId: 32 },
+ "NOISE": { labelId: 0, flagId: 64 },
+ "PMC": { labelId: 0, flagId: 128 },
+ "NOPOS": { labelId: 0, flagId: 256 },
+ "NOROT": { labelId: 0, flagId: 512 },
+ "NOSCL": { labelId: 0, flagId: 1024 }
+ };
+ var _unknownLabelCounter = 0;
+ var _unknownLabels = /* @__PURE__ */ new Map();
+ function decodeLabel(name) {
+ const known = KNOWN_LABELS[name];
+ if (known) return known;
+ let id = _unknownLabels.get(name);
+ if (id === void 0) {
+ id = 1e3 + _unknownLabelCounter++;
+ _unknownLabels.set(name, id);
+ }
+ return { labelId: id, flagId: 0 };
+ }
+ function parseLabelSpecs(seq, lines, isSet) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ requireFieldCount(isSet ? "LABELSET" : "LABELINC", line, p.length, 3);
+ const { labelId, flagId } = decodeLabel(p[2]);
+ const spec = {
+ id: toInt(p[0], isSet ? "LABELSET" : "LABELINC", line),
+ value: toNumber(p[1], isSet ? "LABELSET" : "LABELINC", line),
+ labelId,
+ flagId
+ };
+ if (isSet) seq.labelSets.push(spec);
+ else seq.labelIncs.push(spec);
+ }
+ }
+ function parseSoftDelaySpecs(seq, lines) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ if (p.length < 4) parseError(`DELAYS row has ${p.length} fields, expected at least 4: ${line}`);
+ const hintMatch = line.match(/^\s*\S+\s+\S+\s+\S+\s+\S+\s*(.*)$/);
+ seq.softDelays.push({
+ id: toInt(p[0], "DELAYS", line),
+ numId: toInt(p[1], "DELAYS", line),
+ offset: toNumber(p[2], "DELAYS", line),
+ factor: toNumber(p[3], "DELAYS", line),
+ hint: hintMatch ? hintMatch[1].trim() : ""
+ });
+ }
+ }
+ function parseRFShimSpecs(seq, lines) {
+ for (const line of lines) {
+ const p = splitFields(line);
+ if (p.length < 2) parseError(`RF_SHIMS row has ${p.length} fields, expected at least 2: ${line}`);
+ const nChan = toInt(p[1], "RF_SHIMS", line);
+ requireFieldCount("RF_SHIMS", line, p.length, 2 + nChan * 2);
+ const amps = [];
+ const phases = [];
+ for (let c = 0; c < nChan; c++) {
+ amps.push(toNumber(p[2 + c * 2], "RF_SHIMS", line));
+ phases.push(toNumber(p[2 + c * 2 + 1], "RF_SHIMS", line));
+ }
+ seq.rfShims.push({ id: toInt(p[0], "RF_SHIMS", line), nChannels: nChan, amplitudes: amps, phases });
+ }
+ }
+ function parseShapes(seq, lines) {
+ let i = 0;
+ while (i < lines.length) {
+ const t = lines[i].trim();
+ if (!t || t.startsWith("#") || t.startsWith("[")) {
+ i++;
+ continue;
+ }
+ const m = t.match(/^shape_id\s+(\d+)/);
+ if (!m) {
+ i++;
+ continue;
+ }
+ const shapeId = +m[1];
+ i++;
+ let numSamples = 0;
+ while (i < lines.length) {
+ const l = lines[i].trim();
+ if (!l || l.startsWith("#")) {
+ i++;
+ continue;
+ }
+ const nm = l.match(/^num_samples\s+(\d+)/);
+ if (nm) {
+ numSamples = +nm[1];
+ i++;
+ break;
+ }
+ if (l.match(/^shape_id\s+\d+/) || l.startsWith("[")) break;
+ i++;
+ }
+ if (numSamples <= 0) continue;
+ const vals = [];
+ while (i < lines.length && vals.length < numSamples) {
+ const l = lines[i].trim();
+ if (l.match(/^shape_id\s+\d+/) || l.startsWith("[")) break;
+ if (!l || l.startsWith("#")) {
+ i++;
+ continue;
+ }
+ for (const n of l.split(/\s+/).map(Number).filter((x) => !isNaN(x))) {
+ if (vals.length < numSamples) vals.push(n);
+ }
+ i++;
+ }
+ if (vals.length === 0) continue;
+ storeShape(seq, shapeId, numSamples, vals);
+ }
+ }
+ function storeShape(seq, id, num, raw) {
+ const decompressed = raw.length === num ? new Float64Array(raw) : decompressShape(raw, num);
+ seq.shapes.set(id, { numSamples: num, samples: decompressed });
+ }
+ function extractRasterTimes(seq) {
+ const set = (key, field) => {
+ const v = seq.definitions.get(key);
+ if (v?.length) seq.rasterTimes[field] = v[0];
+ };
+ set("BlockDurationRaster", "blockDurationRaster");
+ set("GradientRasterTime", "gradientRaster");
+ set("RadiofrequencyRasterTime", "rfRaster");
+ set("AdcRasterTime", "adcRaster");
+ }
+ function validateSequence(seq, seenSections) {
+ if (!seenSections.has("VERSION")) parseError("Required [VERSION] section is missing");
+ if (seq.version.major !== 1 || seq.version.minor > 5) {
+ parseError(`Unsupported Pulseq version ${seq.version.major}.${seq.version.minor}.${seq.version.revision}`);
+ }
+ const vc = ver(seq);
+ if (vc >= VER_PRE_14) {
+ requireNumericDefinition(seq, "AdcRasterTime");
+ requireNumericDefinition(seq, "GradientRasterTime");
+ requireNumericDefinition(seq, "RadiofrequencyRasterTime");
+ requireNumericDefinition(seq, "BlockDurationRaster");
+ }
+ if (vc >= VER_V15001) {
+ const required = seq.definitionsRaw.get("RequiredExtensions")?.split(/\s+/).filter(Boolean) ?? [];
+ for (const name of required) {
+ if (extensionNameToType(name) === 999 /* EXT_UNKNOWN */) {
+ parseError(`Unknown required extension '${name}'`);
+ }
+ }
+ }
+ if (!seenSections.has("BLOCKS")) parseError("Required [BLOCKS] section is missing");
+ for (const block of seq.blocks) {
+ if (block.rfId > 0 && !seq.rfs.has(block.rfId)) {
+ parseError(`Block ${block.num} references undefined RF event ${block.rfId}`);
+ }
+ for (const [channel, gradId] of [["GX", block.gxId], ["GY", block.gyId], ["GZ", block.gzId]]) {
+ if (gradId > 0 && !seq.arbitraryGrads.has(gradId) && !seq.trapGrads.has(gradId)) {
+ parseError(`Block ${block.num} references undefined ${channel} gradient event ${gradId}`);
+ }
+ }
+ if (block.adcId > 0 && !seq.adcs.has(block.adcId)) {
+ parseError(`Block ${block.num} references undefined ADC event ${block.adcId}`);
+ }
+ if (block.extId > 0 && !seq.extensions.has(block.extId)) {
+ parseError(`Block ${block.num} references undefined extension list ${block.extId}`);
+ }
+ }
+ for (const ext of seq.extensions.values()) {
+ if (ext.nextId > 0 && !seq.extensions.has(ext.nextId)) {
+ parseError(`Extension list ${ext.id} references undefined next extension ${ext.nextId}`);
+ }
+ const type = seq.extensionTypes.get(ext.type) ?? 999 /* EXT_UNKNOWN */;
+ if (type === 999 /* EXT_UNKNOWN */) continue;
+ if (!extensionPayloadExists(seq, type, ext.ref)) {
+ const name = seq.extensionNames.get(ext.type) ?? `type ${ext.type}`;
+ parseError(`Extension list ${ext.id} references undefined ${name} payload ${ext.ref}`);
+ }
+ }
+ }
+ function requireNumericDefinition(seq, name) {
+ const value = seq.definitions.get(name);
+ if (!value || value.length === 0 || !Number.isFinite(value[0])) {
+ parseError(`Required definition ${name} is not present in the file`);
+ }
+ }
+ function extensionPayloadExists(seq, type, ref) {
+ switch (type) {
+ case 1 /* EXT_TRIGGER */:
+ return seq.triggers.some((v) => v.id === ref);
+ case 2 /* EXT_ROTATION */:
+ return seq.rotations.some((v) => v.id === ref);
+ case 3 /* EXT_LABELSET */:
+ return seq.labelSets.some((v) => v.id === ref);
+ case 4 /* EXT_LABELINC */:
+ return seq.labelIncs.some((v) => v.id === ref);
+ case 5 /* EXT_DELAY */:
+ return seq.softDelays.some((v) => v.id === ref);
+ case 6 /* EXT_RF_SHIM */:
+ return seq.rfShims.some((v) => v.id === ref);
+ case 100 /* EXT_NCO */:
+ return seq.ncos.some((v) => v.id === ref);
+ default:
+ return false;
+ }
+ }
+
+ // src/pulseq/decoder.ts
+ var GAMMA_HZ_T = 42576e3;
+ var DEFAULT_B0_T = 3;
+ function getB0(seq) {
+ const raw = seq.definitions.get("B0");
+ if (raw && Array.isArray(raw) && raw.length > 0) return +raw[0];
+ const raw2 = seq.definitions.get("b0") ?? seq.definitions.get("b_0");
+ if (raw2 && Array.isArray(raw2) && raw2.length > 0) return +raw2[0];
+ return DEFAULT_B0_T;
+ }
+ function effFreqOff(freqOffset, freqPPM, b0) {
+ return freqOffset + freqPPM * 1e-6 * GAMMA_HZ_T * b0;
+ }
+ function effPhaseOff(phaseOffset, phasePPM, b0) {
+ return phaseOffset + phasePPM * 1e-6 * GAMMA_HZ_T * b0;
+ }
+ function decodeAllBlocks(seq) {
+ return decodeBlockRange(seq, 0, seq.blocks.length);
+ }
+ function decodeBlockRange(seq, startBlockIdx, endBlockIdx) {
+ _trigCache.clear();
+ _ncoCache.clear();
+ const totalBlocks = seq.blocks.length;
+ const s = Math.max(0, Math.min(startBlockIdx, totalBlocks));
+ const e = Math.max(s, Math.min(endBlockIdx, totalBlocks));
+ if (s >= e) return [];
+ let cumulative = 0;
+ for (let i = 0; i < Math.min(s, totalBlocks); i++) {
+ cumulative += blockDurationSeconds(seq, seq.blocks[i]);
+ }
+ const decoded = [];
+ for (let i = s; i < e; i++) {
+ const block = seq.blocks[i];
+ const dur = blockDurationSeconds(seq, block);
+ const db = { index: block.num, duration: dur, startTime: cumulative };
+ if (block.rfId > 0) {
+ const rf = seq.rfs.get(block.rfId);
+ if (rf) db.rf = decodeRF(seq, rf, cumulative, dur);
+ }
+ db.gx = decodeGradient(seq, block.gxId, cumulative, dur, "gx");
+ db.gy = decodeGradient(seq, block.gyId, cumulative, dur, "gy");
+ db.gz = decodeGradient(seq, block.gzId, cumulative, dur, "gz");
+ if (block.adcId > 0) {
+ const adc = seq.adcs.get(block.adcId);
+ if (adc) db.adc = decodeADC(adc, cumulative, seq);
+ }
+ if (block.extId > 0) {
+ const ext = seq.extensions.get(block.extId);
+ if (ext) decodeExtensions(seq, ext, db, cumulative);
+ }
+ decoded.push(db);
+ cumulative += dur;
+ }
+ return decoded;
+ }
+ function getTotalDuration(seq) {
+ let total = 0;
+ for (const block of seq.blocks) {
+ total += blockDurationSeconds(seq, block);
+ }
+ return total;
+ }
+ function blockDurationSeconds(seq, block) {
+ if (seq.versionCombined < VER_PRE_14) return block.dur * 1e-6;
+ return block.dur * seq.rasterTimes.blockDurationRaster;
+ }
+ function decodeRF(seq, rf, blockStart, _blockDur) {
+ const raster = seq.rasterTimes.rfRaster;
+ const rfDelay = rf.delay * 1e-6;
+ const rfStart = blockStart + rfDelay;
+ const b0 = getB0(seq);
+ const freqFull = effFreqOff(rf.freqOffset, rf.freqPPM, b0);
+ const phaseFull = effPhaseOff(rf.phaseOffset, rf.phasePPM, b0);
+ const magShape = seq.shapes.get(rf.magShapeId);
+ const nSamples = magShape?.numSamples ?? Math.max(2, Math.round(_blockDur / raster));
+ const mag = magShape ? new Float64Array(magShape.samples) : makeConstant(nSamples, 1);
+ const phShape = seq.shapes.get(rf.phaseShapeId);
+ const ph = phShape ? new Float64Array(phShape.samples) : new Float64Array(mag.length);
+ const timeShape = rf.timeShapeId > 0 ? seq.shapes.get(rf.timeShapeId)?.samples ?? null : null;
+ const n = Math.min(mag.length, ph.length);
+ const t = new Float64Array(n);
+ const amp = new Float64Array(n);
+ const phase = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ t[i] = timeShape ? rfStart + timeShape[i] * raster : rfStart + (i + 0.5) * raster;
+ amp[i] = rf.amplitude * mag[i];
+ const dt = t[i] - rfStart;
+ phase[i] = 2 * Math.PI * ph[i] + phaseFull + 2 * Math.PI * freqFull * dt;
+ }
+ const duration = n > 0 ? t[n - 1] - rfStart + raster : 0;
+ const centerTime = rf.center >= 0 ? blockStart + rfDelay + rf.center * 1e-6 : estimateRfPeakTime(t, amp, rfStart, duration);
+ let use = rf.use || "";
+ if (!use || use === "u") {
+ let faDeg = 0;
+ for (let i = 1; i < n; i++) {
+ const dt = t[i] - t[i - 1];
+ faDeg += 360 * (amp[i] + amp[i - 1]) * 0.5 * dt;
+ }
+ use = faDeg >= 120 ? "r" : "e";
+ }
+ return {
+ blockIndex: rf.id,
+ startTime: rfStart,
+ centerTime,
+ duration,
+ timePoints: t,
+ magnitude: amp,
+ phase,
+ amplitude: rf.amplitude,
+ freqOffset: freqFull,
+ phaseOffset: phaseFull,
+ use
+ };
+ }
+ function decodeGradient(seq, gradId, blockStart, blockDur, channel) {
+ if (gradId <= 0) return zeroGradient(blockStart, blockDur, channel);
+ const trap = seq.trapGrads.get(gradId);
+ if (trap) return decodeTrap(trap, blockStart, channel);
+ const arb = seq.arbitraryGrads.get(gradId);
+ if (arb) return decodeArb(seq, arb, blockStart, channel);
+ return zeroGradient(blockStart, blockDur, channel);
+ }
+ function zeroGradient(t0, dur, ch) {
+ return {
+ blockIndex: 0,
+ startTime: t0,
+ duration: dur,
+ timePoints: new Float64Array([t0, t0 + dur]),
+ waveform: new Float64Array([0, 0]),
+ amplitude: 0,
+ type: "none",
+ channel: ch
+ };
+ }
+ function decodeTrap(trap, blockStart, ch) {
+ const rise = trap.rise * 1e-6;
+ const flat = trap.flat * 1e-6;
+ const fall = trap.fall * 1e-6;
+ const delay = trap.delay * 1e-6;
+ const gradStart = blockStart + delay;
+ const tRel = [0, rise, rise + flat, rise + flat + fall];
+ const wfRel = [0, trap.amplitude, trap.amplitude, 0];
+ if (delay > 0) {
+ const tp2 = new Float64Array(5);
+ const wf2 = new Float64Array(5);
+ tp2[0] = blockStart;
+ wf2[0] = 0;
+ for (let i = 0; i < 4; i++) {
+ tp2[i + 1] = gradStart + tRel[i];
+ wf2[i + 1] = wfRel[i];
+ }
+ return {
+ blockIndex: trap.id,
+ startTime: blockStart,
+ duration: delay + rise + flat + fall,
+ timePoints: tp2,
+ waveform: wf2,
+ amplitude: trap.amplitude,
+ type: "trap",
+ channel: ch
+ };
+ }
+ const tp = new Float64Array(4);
+ const wf = new Float64Array(4);
+ for (let i = 0; i < 4; i++) {
+ tp[i] = gradStart + tRel[i];
+ wf[i] = wfRel[i];
+ }
+ return {
+ blockIndex: trap.id,
+ startTime: blockStart,
+ duration: rise + flat + fall,
+ timePoints: tp,
+ waveform: wf,
+ amplitude: trap.amplitude,
+ type: "trap",
+ channel: ch
+ };
+ }
+ function decodeArb(seq, arb, blockStart, ch) {
+ const shape = seq.shapes.get(arb.shapeId);
+ if (!shape) return zeroGradient(blockStart, 0, ch);
+ const raster = seq.rasterTimes.gradientRaster;
+ const delay = arb.delay * 1e-6;
+ const gradStart = blockStart + delay;
+ const n = shape.numSamples;
+ const oversampled = arb.timeId === -1;
+ const timeShape = arb.timeId > 0 ? seq.shapes.get(arb.timeId)?.samples ?? null : null;
+ if (timeShape) {
+ const tp2 = new Float64Array(n);
+ const wf2 = new Float64Array(n);
+ for (let i = 0; i < n; i++) {
+ tp2[i] = gradStart + timeShape[i] * raster;
+ wf2[i] = arb.amplitude * shape.samples[i];
+ }
+ const dur2 = n > 0 ? tp2[n - 1] - blockStart + raster : delay;
+ return {
+ blockIndex: arb.id,
+ startTime: blockStart,
+ duration: dur2,
+ timePoints: tp2,
+ waveform: wf2,
+ amplitude: arb.amplitude,
+ type: "arb",
+ channel: ch
+ };
+ }
+ const tp = new Float64Array(n + 2);
+ const wf = new Float64Array(n + 2);
+ tp[0] = gradStart;
+ wf[0] = edgeAmplitude(arb.first, arb.amplitude, shape.samples, true);
+ if (oversampled) {
+ const dt = raster * 0.5;
+ for (let i = 0; i < n; i++) {
+ tp[i + 1] = gradStart + (i + 1) * dt;
+ wf[i + 1] = arb.amplitude * shape.samples[i];
+ }
+ tp[n + 1] = gradStart + (n + 1) * dt;
+ } else {
+ for (let i = 0; i < n; i++) {
+ tp[i + 1] = gradStart + (i + 0.5) * raster;
+ wf[i + 1] = arb.amplitude * shape.samples[i];
+ }
+ tp[n + 1] = gradStart + n * raster;
+ }
+ wf[wf.length - 1] = edgeAmplitude(arb.last, arb.amplitude, shape.samples, false);
+ const dur = tp[tp.length - 1] - blockStart;
+ return {
+ blockIndex: arb.id,
+ startTime: blockStart,
+ duration: dur,
+ timePoints: tp,
+ waveform: wf,
+ amplitude: arb.amplitude,
+ type: "arb",
+ channel: ch
+ };
+ }
+ function edgeAmplitude(stored, amplitude, samples, first) {
+ let value;
+ if (Number.isFinite(stored)) {
+ value = stored;
+ if (Math.abs(value) > 1 + 1e-6 && Math.abs(amplitude) > 0) value /= amplitude;
+ } else if (samples.length === 0) {
+ value = 0;
+ } else if (samples.length === 1) {
+ value = samples[0];
+ } else if (first) {
+ value = 0.5 * (3 * samples[0] - samples[1]);
+ } else {
+ value = 0.5 * (3 * samples[samples.length - 1] - samples[samples.length - 2]);
+ }
+ return value * amplitude;
+ }
+ function decodeADC(adc, blockStart, seq) {
+ const b0 = getB0(seq);
+ const freqFull = effFreqOff(adc.freqOffset, adc.freqPPM, b0);
+ const phaseFull = effPhaseOff(adc.phaseOffset, adc.phasePPM, b0);
+ return {
+ blockIndex: adc.id,
+ startTime: blockStart,
+ numSamples: adc.numSamples,
+ dwell: adc.dwell * 1e-9,
+ // ns β s
+ delay: adc.delay * 1e-6,
+ // Β΅s β s
+ freqOffset: freqFull,
+ phaseOffset: phaseFull
+ };
+ }
+ var _trigCache = /* @__PURE__ */ new Map();
+ var _ncoCache = /* @__PURE__ */ new Map();
+ function decodeExtensions(seq, ext, db, blockStart) {
+ const visited = /* @__PURE__ */ new Set();
+ let cur = ext;
+ while (cur && !visited.has(cur.id)) {
+ visited.add(cur.id);
+ const type = seq.extensionTypes.get(cur.type) ?? 999 /* EXT_UNKNOWN */;
+ if (type === 1 /* EXT_TRIGGER */) {
+ let cached = _trigCache.get(cur.id);
+ if (!cached) {
+ const trigger = findById(seq.triggers, cur.ref);
+ if (trigger) {
+ cached = {
+ blockIndex: trigger.id,
+ startTime: 0,
+ channel: trigger.channel,
+ delay: trigger.delay * 1e-6,
+ duration: trigger.duration * 1e-6
+ };
+ _trigCache.set(cur.id, cached);
+ }
+ }
+ if (cached) {
+ if (!db.triggers) db.triggers = [];
+ db.triggers.push({ ...cached, startTime: blockStart });
+ }
+ } else if (type === 100 /* EXT_NCO */) {
+ let cached = _ncoCache.get(cur.id);
+ if (!cached) {
+ const nco = findById(seq.ncos, cur.ref);
+ if (nco) {
+ cached = {
+ blockIndex: nco.id,
+ startTime: 0,
+ channel: nco.channel,
+ frequency: nco.frequency,
+ phase: nco.phase,
+ delay: nco.delay * 1e-6,
+ duration: nco.duration * 1e-6
+ };
+ _ncoCache.set(cur.id, cached);
+ }
+ }
+ if (cached) {
+ if (!db.nco) db.nco = [];
+ db.nco.push({ ...cached, startTime: blockStart });
+ }
+ } else if (type === 2 /* EXT_ROTATION */) {
+ const rotation = findById(seq.rotations, cur.ref);
+ if (rotation) db.rotation = { id: rotation.id, values: [...rotation.values] };
+ } else if (type === 3 /* EXT_LABELSET */) {
+ const label = findById(seq.labelSets, cur.ref);
+ if (label) {
+ if (!db.labelSets) db.labelSets = [];
+ db.labelSets.push({ ...label });
+ }
+ } else if (type === 4 /* EXT_LABELINC */) {
+ const label = findById(seq.labelIncs, cur.ref);
+ if (label) {
+ if (!db.labelIncs) db.labelIncs = [];
+ db.labelIncs.push({ ...label });
+ }
+ } else if (type === 5 /* EXT_DELAY */) {
+ const delay = findById(seq.softDelays, cur.ref);
+ if (delay) db.softDelay = { ...delay };
+ } else if (type === 6 /* EXT_RF_SHIM */) {
+ const shim = findById(seq.rfShims, cur.ref);
+ if (shim) {
+ db.rfShim = {
+ id: shim.id,
+ nChannels: shim.nChannels,
+ amplitudes: [...shim.amplitudes],
+ phases: [...shim.phases]
+ };
+ }
+ }
+ cur = cur.nextId > 0 ? seq.extensions.get(cur.nextId) : void 0;
+ }
+ }
+ function makeConstant(n, value) {
+ const a = new Float64Array(Math.max(n, 2));
+ a.fill(value);
+ return a;
+ }
+ function estimateRfPeakTime(timePoints, magnitude, startTime, duration) {
+ if (!timePoints.length || !magnitude.length) return startTime + duration * 0.5;
+ let peak = Math.abs(magnitude[0]);
+ for (let i = 1; i < magnitude.length; i++) {
+ const v = Math.abs(magnitude[i]);
+ if (v > peak) peak = v;
+ }
+ const threshold = Math.abs(peak) * 0.99999;
+ let firstPeak = -1;
+ let lastPeak = -1;
+ for (let i = 0; i < magnitude.length; i++) {
+ if (Math.abs(magnitude[i]) >= threshold) {
+ if (firstPeak < 0) firstPeak = i;
+ lastPeak = i;
+ }
+ }
+ if (firstPeak < 0 || lastPeak < 0) return startTime + duration * 0.5;
+ return 0.5 * (timePoints[Math.min(firstPeak, timePoints.length - 1)] + timePoints[Math.min(lastPeak, timePoints.length - 1)]);
+ }
+ function findById(items, id) {
+ return items.find((item) => item.id === id);
+ }
+
+ // src/pulseq/kspace.ts
+ function calculateKspace(blocks, gradientRaster, totalDuration, trajectoryDelay = 0, _options) {
+ if (!blocks.length || !gradientRaster || gradientRaster <= 0) return null;
+ const GR = gradientRaster;
+ const RF = _options?.rfRaster && _options.rfRaster > 0 ? _options.rfRaster : 1e-6;
+ const tacc = 1e-10;
+ const gradientSupport = _options?.gradientSupport ?? "endpoints";
+ const excT = [], refT = [];
+ const gradTimes = [];
+ let totalAdcSamples = 0;
+ for (const b of blocks) {
+ if (b.adc) totalAdcSamples += b.adc.numSamples;
+ }
+ const adcT = new Float64Array(totalAdcSamples);
+ let adcIdx = 0;
+ for (const b of blocks) {
+ collectGradientSupport(b.gx, gradTimes, gradientSupport);
+ collectGradientSupport(b.gy, gradTimes, gradientSupport);
+ collectGradientSupport(b.gz, gradTimes, gradientSupport);
+ if (b.rf) {
+ const iso = Number.isFinite(b.rf.centerTime) ? b.rf.centerTime : b.rf.startTime + b.rf.duration * 0.5;
+ const u = b.rf.use || "";
+ if (u === "e" || u === "" || u === "u") excT.push(iso);
+ else if (u === "r") refT.push(iso);
+ }
+ if (b.adc) {
+ const t0 = b.adc.startTime + b.adc.delay;
+ const dwell = b.adc.dwell;
+ const nSamp = b.adc.numSamples;
+ for (let s = 0; s < nSamp; s++)
+ adcT[adcIdx++] = t0 + (s + 0.5) * dwell + trajectoryDelay;
+ }
+ }
+ const cand = [];
+ const pushC = (t) => {
+ if (isFinite(t) && t >= -tacc) cand.push(Math.max(0, tacc * Math.round(t / tacc)));
+ };
+ for (const t of gradTimes) pushC(t);
+ for (const t of excT) {
+ pushC(t);
+ pushC(t - RF);
+ pushC(t - 2 * RF);
+ }
+ for (const t of refT) {
+ pushC(t);
+ pushC(t - RF);
+ }
+ for (const t of adcT) pushC(t);
+ pushC(0);
+ pushC(totalDuration);
+ if (totalDuration > 0) {
+ const nS = Math.max(1, Math.round(totalDuration / GR));
+ for (let i = 0; i <= nS; i++) pushC(i * GR);
+ }
+ if (cand.length === 0) return null;
+ cand.sort((a, b) => a - b);
+ const grid = [];
+ for (let i = 0; i < cand.length; i++) {
+ if (i === 0 || cand[i] - cand[i - 1] > tacc * 0.5) grid.push(cand[i]);
+ }
+ const N = grid.length;
+ if (N < 2) return null;
+ if (_options?.maxGridPoints && N > _options.maxGridPoints) return null;
+ const gx = new Float64Array(N), gy = new Float64Array(N), gz = new Float64Array(N);
+ const edges = [0];
+ let cum = 0;
+ for (const b of blocks) {
+ cum += b.duration;
+ edges.push(cum);
+ }
+ for (let i = 0; i < N; i++) {
+ const t = grid[i];
+ const bi = blockIdx(t, edges);
+ if (bi >= 0 && bi < blocks.length) {
+ const block = blocks[bi];
+ const localX = gradVal(block.gx, t);
+ const localY = gradVal(block.gy, t);
+ const localZ = gradVal(block.gz, t);
+ const rotated = rotateGradient(block, localX, localY, localZ);
+ gx[i] = rotated[0];
+ gy[i] = rotated[1];
+ gz[i] = rotated[2];
+ }
+ }
+ const kx = new Float64Array(N), ky = new Float64Array(N), kz = new Float64Array(N);
+ for (let i = 1; i < N; i++) {
+ const dt = grid[i] - grid[i - 1];
+ if (dt <= 0) {
+ kx[i] = kx[i - 1];
+ ky[i] = ky[i - 1];
+ kz[i] = kz[i - 1];
+ continue;
+ }
+ const gxm = 0.5 * (gx[i - 1] + gx[i]), gym = 0.5 * (gy[i - 1] + gy[i]), gzm = 0.5 * (gz[i - 1] + gz[i]);
+ kx[i] = kx[i - 1] + gxm * dt;
+ ky[i] = ky[i - 1] + gym * dt;
+ kz[i] = kz[i - 1] + gzm * dt;
+ }
+ const eIdx = [], rIdx = [];
+ for (const t of excT) {
+ const i = timeIdx(t, grid);
+ if (i >= 0) eIdx.push(i);
+ }
+ for (const t of refT) {
+ const i = timeIdx(t, grid);
+ if (i >= 0) rIdx.push(i);
+ }
+ eIdx.sort((a, b) => a - b);
+ rIdx.sort((a, b) => a - b);
+ const bounds = [0];
+ for (const i of eIdx) bounds.push(i);
+ for (const i of rIdx) bounds.push(i);
+ bounds.push(N - 1);
+ bounds.sort((a, b) => a - b);
+ const bUniq = [bounds[0]];
+ for (let i = 1; i < bounds.length; i++) if (bounds[i] !== bUniq[bUniq.length - 1]) bUniq.push(bounds[i]);
+ let dkX = -kx[0], dkY = -ky[0], dkZ = -kz[0];
+ let pE = 0, pR = 0;
+ for (let s = 0; s < bUniq.length - 1; s++) {
+ const st = bUniq[s], en = bUniq[s + 1];
+ if (pE < eIdx.length && eIdx[pE] === st) {
+ dkX = -kx[st];
+ dkY = -ky[st];
+ dkZ = -kz[st];
+ pE++;
+ } else if (pR < rIdx.length && rIdx[pR] === st) {
+ dkX = -2 * kx[st] - dkX;
+ dkY = -2 * ky[st] - dkY;
+ dkZ = -2 * kz[st] - dkZ;
+ pR++;
+ }
+ for (let j = st; j < en; j++) {
+ kx[j] += dkX;
+ ky[j] += dkY;
+ kz[j] += dkZ;
+ }
+ }
+ kx[N - 1] += dkX;
+ ky[N - 1] += dkY;
+ kz[N - 1] += dkZ;
+ const kxP = new Float64Array(kx), kyP = new Float64Array(ky), kzP = new Float64Array(kz);
+ for (const i of eIdx) {
+ if (i > 0) {
+ kxP[i - 1] = NaN;
+ kyP[i - 1] = NaN;
+ kzP[i - 1] = NaN;
+ }
+ }
+ const nA = adcT.length;
+ const kxA = new Float64Array(nA), kyA = new Float64Array(nA), kzA = new Float64Array(nA);
+ for (let a = 0; a < nA; a++) {
+ kxA[a] = interp(kx, grid, adcT[a]);
+ kyA[a] = interp(ky, grid, adcT[a]);
+ kzA[a] = interp(kz, grid, adcT[a]);
+ }
+ return { ktraj: [kxP, kyP, kzP], t_ktraj: new Float64Array(grid), ktraj_adc: [kxA, kyA, kzA], t_adc: new Float64Array(adcT) };
+ }
+ function collectGradientSupport(g, support, mode) {
+ if (!g || g.type === "none" || !g.timePoints || g.timePoints.length < 2) return;
+ if (mode === "all") {
+ for (let i = 0; i < g.timePoints.length; i++) support.push(g.timePoints[i]);
+ return;
+ }
+ support.push(g.timePoints[0], g.timePoints[g.timePoints.length - 1]);
+ }
+ function gradVal(g, t) {
+ if (!g || g.type === "none") return 0;
+ const tp = g.timePoints, wf = g.waveform;
+ if (!tp || tp.length < 2) return 0;
+ if (t < tp[0] || t > tp[tp.length - 1]) return 0;
+ let lo = 0, hi = tp.length - 1;
+ while (hi - lo > 1) {
+ const m = lo + hi >> 1;
+ if (tp[m] <= t) lo = m;
+ else hi = m;
+ }
+ const s = tp[hi] - tp[lo];
+ if (s <= 0) return wf[lo];
+ return wf[lo] + (wf[hi] - wf[lo]) * (t - tp[lo]) / s;
+ }
+ function blockIdx(t, edges) {
+ let lo = 0, hi = edges.length - 1;
+ while (lo < hi) {
+ const m = lo + hi >> 1;
+ if (edges[m] <= t + 1e-12) lo = m + 1;
+ else hi = m;
+ }
+ return Math.max(0, lo - 1);
+ }
+ function timeIdx(t, g) {
+ let lo = 0, hi = g.length;
+ while (lo < hi) {
+ const m = lo + hi >> 1;
+ if (g[m] < t - 1e-12) lo = m + 1;
+ else hi = m;
+ }
+ return lo < g.length ? lo : -1;
+ }
+ function interp(d, g, t) {
+ const n = g.length;
+ if (n === 0) return 0;
+ let lo = 0, hi = n;
+ while (lo < hi) {
+ const m = lo + hi >> 1;
+ if (g[m] < t) lo = m + 1;
+ else hi = m;
+ }
+ if (lo === 0) return d[0];
+ if (lo >= n) return d[n - 1];
+ if (Math.abs(g[lo] - t) < 1e-12) return d[lo];
+ const i0 = lo - 1, i1 = lo, dt = g[i1] - g[i0];
+ if (dt <= 0) return d[i1];
+ return d[i0] + (d[i1] - d[i0]) * (t - g[i0]) / dt;
+ }
+ function rotateGradient(block, gx, gy, gz) {
+ const values = block.rotation?.values;
+ if (!values) return [gx, gy, gz];
+ if (values.length === 4) {
+ const [w, x, y, z] = values;
+ const r00 = 1 - 2 * y * y - 2 * z * z;
+ const r01 = 2 * x * y - 2 * w * z;
+ const r02 = 2 * x * z + 2 * w * y;
+ const r10 = 2 * x * y + 2 * w * z;
+ const r11 = 1 - 2 * x * x - 2 * z * z;
+ const r12 = 2 * y * z - 2 * w * x;
+ const r20 = 2 * x * z - 2 * w * y;
+ const r21 = 2 * y * z + 2 * w * x;
+ const r22 = 1 - 2 * x * x - 2 * y * y;
+ return [
+ r00 * gx + r01 * gy + r02 * gz,
+ r10 * gx + r11 * gy + r12 * gz,
+ r20 * gx + r21 * gy + r22 * gz
+ ];
+ }
+ if (values.length === 9) {
+ return [
+ values[0] * gx + values[1] * gy + values[2] * gz,
+ values[3] * gx + values[4] * gy + values[5] * gz,
+ values[6] * gx + values[7] * gy + values[8] * gz
+ ];
+ }
+ return [gx, gy, gz];
+ }
+
+ // src/pulseq/m1.ts
+ var TIME_EPS = 1e-15;
+ function calculateM1(blocks, gradientRaster, options = {}) {
+ const referenceMode = normalizeReferenceMode(options.referenceMode);
+ if (!blocks.length) {
+ return invalidM1("Empty or invalid block list.", referenceMode);
+ }
+ const gx = collectGradientSeries(blocks, "gx");
+ const gy = collectGradientSeries(blocks, "gy");
+ const gz = collectGradientSeries(blocks, "gz");
+ const ranges = [gx, gy, gz].filter((series) => series.time.length > 0).map((series) => [series.time[0], series.time[series.time.length - 1]]);
+ if (!ranges.length) {
+ return invalidM1("No gradient waveform available for M1.", referenceMode);
+ }
+ const tMin = Math.min(...ranges.map((range) => range[0]));
+ const tMax = Math.max(...ranges.map((range) => range[1]));
+ if (!Number.isFinite(tMin) || !Number.isFinite(tMax) || tMax < tMin) {
+ return invalidM1("Invalid gradient time range for M1.", referenceMode);
+ }
+ const warnings = [];
+ const rfEvents = collectRfEvents(blocks, warnings);
+ const excitationTimes = rfEvents.filter((rf) => rf.use === "e").map((rf) => rf.tSec);
+ const refocusingTimes = rfEvents.filter((rf) => rf.use === "r").map((rf) => rf.tSec);
+ const events = buildWalkerEvents(rfEvents);
+ let recentExcCount = 0;
+ let lastExcT = -1e9;
+ for (const rf of rfEvents) {
+ if (rf.use === "e") {
+ if (rf.tSec - lastExcT < 0.1) recentExcCount++;
+ lastExcT = rf.tSec;
+ }
+ }
+ if (recentExcCount > 8) {
+ warnings.push(
+ `Sequence shows ${recentExcCount} closely-spaced (<100 ms) excitation events. This pattern is consistent with a steady-state sequence for which the simplified reset/flip bookkeeping does NOT model coherent pathway interference. Treat the M1 curve as advisory only.`
+ );
+ }
+ if (!excitationTimes.length) {
+ warnings.push(`No excitation RF events found in sequence. M1 will be integrated from t=${tMin.toFixed(6)} s with no signal basis.`);
+ }
+ const rasterSec = gradientRaster > 0 ? gradientRaster : 1e-5;
+ if (rasterSec <= 0) {
+ return invalidM1("gradientRaster must be positive.", referenceMode);
+ }
+ const samples = buildSampleTimes(tMin, tMax, rasterSec);
+ const x = walkM1(gx, samples, events, excitationTimes, tMin, referenceMode);
+ const y = walkM1(gy, samples, events, excitationTimes, tMin, referenceMode);
+ const z = walkM1(gz, samples, events, excitationTimes, tMin, referenceMode);
+ if (x.t.length !== y.t.length || x.t.length !== z.t.length) {
+ warnings.push(`Internal warning: per-axis M1 output sizes disagree (${x.t.length}, ${y.t.length}, ${z.t.length}). Plot may be inconsistent.`);
+ }
+ return {
+ valid: true,
+ ok: true,
+ referenceMode,
+ tSec: new Float64Array(x.t),
+ m1x: new Float64Array(x.m1),
+ m1y: new Float64Array(y.m1),
+ m1z: new Float64Array(z.m1),
+ warnings,
+ excitationTimesSec: new Float64Array(excitationTimes),
+ refocusingTimesSec: new Float64Array(refocusingTimes)
+ };
+ }
+ function invalidM1(error, referenceMode = "rfCenter") {
+ return {
+ valid: false,
+ ok: false,
+ referenceMode,
+ error,
+ tSec: new Float64Array(),
+ m1x: new Float64Array(),
+ m1y: new Float64Array(),
+ m1z: new Float64Array(),
+ warnings: [],
+ excitationTimesSec: new Float64Array(),
+ refocusingTimesSec: new Float64Array()
+ };
+ }
+ function normalizeReferenceMode(mode) {
+ return mode === "observationTime" ? "observationTime" : "rfCenter";
+ }
+ function collectGradientSeries(blocks, channel) {
+ const time = [];
+ const value = [];
+ for (const block of blocks) {
+ const grad = block[channel];
+ if (!grad?.timePoints || !grad.waveform) continue;
+ const n = Math.min(grad.timePoints.length, grad.waveform.length);
+ for (let i = 0; i < n; i++) {
+ time.push(grad.timePoints[i]);
+ value.push(grad.waveform[i]);
+ }
+ }
+ return sanitizeGradientSeries(time, value);
+ }
+ function sanitizeGradientSeries(time, value) {
+ const pairs = [];
+ const n = Math.min(time.length, value.length);
+ for (let i = 0; i < n; i++) {
+ const t = time[i];
+ const v = value[i];
+ if (Number.isFinite(t) && Number.isFinite(v)) pairs.push([t, v]);
+ }
+ pairs.sort((a, b) => a[0] - b[0]);
+ const outT = [];
+ const outV = [];
+ for (const [t, v] of pairs) {
+ const last = outT.length - 1;
+ if (last >= 0 && Math.abs(t - outT[last]) <= TIME_EPS) {
+ outV[last] = 0.5 * (outV[last] + v);
+ continue;
+ }
+ outT.push(t);
+ outV.push(v);
+ }
+ return { time: outT, value: outV };
+ }
+ function collectRfEvents(blocks, warnings) {
+ const events = [];
+ for (const block of blocks) {
+ if (!block.rf) continue;
+ const use = classifyRfUse(block.rf.use);
+ if (!use) continue;
+ const rec = { tSec: block.rf.centerTime, use };
+ events.push(rec);
+ if (use === "u") {
+ warnings.push(`Unknown RF use 'u' at t=${rec.tSec.toFixed(6)} s; M1 bookkeeping treats it as no-op.`);
+ } else if (use === "p") {
+ warnings.push(
+ `Preparation module 'p' at t=${rec.tSec.toFixed(6)} s; treated as M1 reset (simplified handling; prep modules that preserve phase encoding will give wrong results).`
+ );
+ }
+ }
+ events.sort((a, b) => a.tSec - b.tSec);
+ return events;
+ }
+ function classifyRfUse(raw) {
+ const c = (raw || "u").toLowerCase();
+ if (c === "e" || c === "r" || c === "s" || c === "i" || c === "p") return c;
+ return "u";
+ }
+ function buildWalkerEvents(rfs) {
+ const events = [];
+ for (const rf of rfs) {
+ if (rf.use === "i" || rf.use === "u") continue;
+ events.push({
+ tSec: rf.tSec,
+ kind: rf.use === "r" ? "flip" : "reset"
+ });
+ }
+ events.sort((a, b) => {
+ if (a.tSec !== b.tSec) return a.tSec - b.tSec;
+ return a.kind === "reset" && b.kind === "flip" ? -1 : 1;
+ });
+ return events;
+ }
+ function buildSampleTimes(tMin, tMax, rasterSec) {
+ const samples = [];
+ const nSamples = Math.floor((tMax - tMin) / rasterSec) + 1;
+ for (let i = 0; i < nSamples; i++) samples.push(tMin + i * rasterSec);
+ if (!samples.length || samples[samples.length - 1] < tMax - TIME_EPS) samples.push(tMax);
+ return samples;
+ }
+ function walkM1(gradient, samples, events, excitationTimes, tMin, referenceMode) {
+ const outT = [];
+ const outM1 = [];
+ let sign = 1;
+ let tReset = excitationTimes.length ? excitationTimes[0] : tMin;
+ if (samples.length && samples[0] < tReset) tReset = samples[0];
+ let currentT = tReset;
+ let unsignedM0 = 0;
+ let unsignedM1 = 0;
+ const reportedM1At = (t) => {
+ if (referenceMode === "observationTime") return sign * (unsignedM1 - (t - tReset) * unsignedM0);
+ return sign * unsignedM1;
+ };
+ const advanceTo = (targetT) => {
+ if (!(targetT > currentT + TIME_EPS)) return;
+ while (currentT < targetT - TIME_EPS) {
+ let nextT = nextGradientBreakpoint(gradient.time, currentT, targetT);
+ if (!(nextT > currentT)) nextT = targetT;
+ const ga = sampleGradientAt(gradient, currentT);
+ const gb = sampleGradientAt(gradient, nextT);
+ const [m0Seg, m1Seg] = integrateLinearSegment(currentT, nextT, tReset, ga, gb);
+ unsignedM0 += m0Seg;
+ unsignedM1 += m1Seg;
+ currentT = nextT;
+ }
+ };
+ let ei = 0;
+ let si = 0;
+ while (ei < events.length || si < samples.length) {
+ const nextEvtT = ei < events.length ? events[ei].tSec : Number.POSITIVE_INFINITY;
+ const nextSampT = si < samples.length ? samples[si] : Number.POSITIVE_INFINITY;
+ if (nextEvtT <= nextSampT) {
+ advanceTo(nextEvtT);
+ if (events[ei].kind === "reset") {
+ if (!outT.length || outT[outT.length - 1] < nextEvtT - TIME_EPS) {
+ outT.push(nextEvtT);
+ outM1.push(0);
+ } else {
+ outT[outT.length - 1] = nextEvtT;
+ outM1[outM1.length - 1] = 0;
+ }
+ sign = 1;
+ tReset = nextEvtT;
+ currentT = nextEvtT;
+ unsignedM0 = 0;
+ unsignedM1 = 0;
+ } else {
+ outT.push(nextEvtT);
+ outM1.push(reportedM1At(nextEvtT));
+ sign = -sign;
+ }
+ ei++;
+ } else {
+ advanceTo(nextSampT);
+ outT.push(nextSampT);
+ outM1.push(reportedM1At(nextSampT));
+ si++;
+ }
+ }
+ return { t: outT, m1: outM1 };
+ }
+ function sampleGradientAt(gradient, t) {
+ const n = gradient.time.length;
+ if (n <= 0 || t < gradient.time[0] || t > gradient.time[n - 1]) return 0;
+ if (n === 1 || t <= gradient.time[0]) return gradient.value[0];
+ if (t >= gradient.time[n - 1]) return gradient.value[n - 1];
+ let lo = 0;
+ let hi = n - 1;
+ while (hi - lo > 1) {
+ const mid = lo + hi >> 1;
+ if (gradient.time[mid] <= t) lo = mid;
+ else hi = mid;
+ }
+ const t0 = gradient.time[lo];
+ const t1 = gradient.time[hi];
+ if (!(t1 > t0)) return gradient.value[lo];
+ const alpha = (t - t0) / (t1 - t0);
+ return gradient.value[lo] + alpha * (gradient.value[hi] - gradient.value[lo]);
+ }
+ function nextGradientBreakpoint(times, t, target) {
+ if (times.length <= 1 || t >= times[times.length - 1]) return target;
+ let lo = 0;
+ let hi = times.length;
+ const threshold = t + TIME_EPS;
+ while (lo < hi) {
+ const mid = lo + hi >> 1;
+ if (times[mid] <= threshold) lo = mid + 1;
+ else hi = mid;
+ }
+ return lo < times.length ? Math.min(target, times[lo]) : target;
+ }
+ function integrateLinearSegment(a, b, tRef, ga, gb) {
+ const h = b - a;
+ if (!(h > 0)) return [0, 0];
+ const slope = (gb - ga) / h;
+ const aRel = a - tRef;
+ const m0 = ga * h + 0.5 * slope * h * h;
+ const m1 = ga * (aRel * h + 0.5 * h * h) + slope * (0.5 * aRel * h * h + h * h * h / 3);
+ return [m0, m1];
+ }
+
+ // src/pulseq/pns.ts
+ var GAMMA_HZ_PER_T = 42576e3;
+ var TIME_EPS2 = 1e-15;
+ function parsePnsHardwareAsc(text) {
+ const asc = parseAscText(text);
+ const prefix = resolvePnsPrefix(asc);
+ const x = getAxisHardware(
+ asc,
+ `${prefix}flGSWDTauX`,
+ `${prefix}flGSWDAX`,
+ `${prefix}flGSWDStimulationLimitX`,
+ `${prefix}flGSWDStimulationThresholdX`,
+ [
+ "asGPAParameters[0].sGCParameters.flGScaleFactorX",
+ "asGPAParameters.sGCParameters.flGScaleFactorX",
+ "flGScaleFactorX",
+ "flGCGScaleFactorX",
+ "GScaleFactorX"
+ ]
+ );
+ const y = getAxisHardware(
+ asc,
+ `${prefix}flGSWDTauY`,
+ `${prefix}flGSWDAY`,
+ `${prefix}flGSWDStimulationLimitY`,
+ `${prefix}flGSWDStimulationThresholdY`,
+ [
+ "asGPAParameters[0].sGCParameters.flGScaleFactorY",
+ "asGPAParameters.sGCParameters.flGScaleFactorY",
+ "flGScaleFactorY",
+ "flGCGScaleFactorY",
+ "GScaleFactorY"
+ ]
+ );
+ const z = getAxisHardware(
+ asc,
+ `${prefix}flGSWDTauZ`,
+ `${prefix}flGSWDAZ`,
+ `${prefix}flGSWDStimulationLimitZ`,
+ `${prefix}flGSWDStimulationThresholdZ`,
+ [
+ "asGPAParameters[0].sGCParameters.flGScaleFactorZ",
+ "asGPAParameters.sGCParameters.flGScaleFactorZ",
+ "flGScaleFactorZ",
+ "flGCGScaleFactorZ",
+ "GScaleFactorZ"
+ ]
+ );
+ if (!hasValidWeights(x) || !hasValidWeights(y) || !hasValidWeights(z)) {
+ throw new Error("ASC hardware coefficients are invalid (a1+a2+a3 or stim limit).");
+ }
+ return { x, y, z, valid: true };
+ }
+ function calculatePns(blocks, gradientRaster, hardware, gammaHzPerT = GAMMA_HZ_PER_T) {
+ if (!hardware.valid) return invalidPns("PNS hardware is not initialized.");
+ if (!blocks.length) return invalidPns("No sequence loaded.");
+ if (gradientRaster <= 0 || gammaHzPerT <= 0) return invalidPns("Missing GradientRasterTime or gamma.");
+ const dtSec = gradientRaster;
+ const waves = [
+ collectGradientSeries2(blocks, "gx"),
+ collectGradientSeries2(blocks, "gy"),
+ collectGradientSeries2(blocks, "gz")
+ ];
+ const nonEmpty = waves.filter((wave) => wave.time.length > 0);
+ if (!nonEmpty.length) return invalidPns("No gradient waveform available for PNS.");
+ const tFirst = Math.min(...nonEmpty.map((wave) => wave.time[0]));
+ const tLast = Math.max(...nonEmpty.map((wave) => wave.time[wave.time.length - 1]));
+ if (!Number.isFinite(tFirst) || !Number.isFinite(tLast) || tLast <= tFirst) {
+ return invalidPns("No gradient waveform available for PNS.");
+ }
+ let ntMin = Math.floor(tFirst / dtSec + Number.EPSILON) + 0.5;
+ const ntMax = Math.ceil(tLast / dtSec - Number.EPSILON) - 0.5;
+ if (ntMin < 0.5) ntMin = 0.5;
+ if (ntMax < ntMin) return invalidPns("Unable to build regular PNS raster.");
+ const nSamples = Math.floor(ntMax - ntMin + 1);
+ if (nSamples < 2) return invalidPns("Too few samples for PNS computation.");
+ const tAxis = new Float64Array(nSamples);
+ const gxTpm = new Float64Array(nSamples);
+ const gyTpm = new Float64Array(nSamples);
+ const gzTpm = new Float64Array(nSamples);
+ for (let i = 0; i < nSamples; i++) {
+ const tSec = (ntMin + i) * dtSec;
+ tAxis[i] = tSec;
+ gxTpm[i] = interpLinearZero(waves[0], tSec) / gammaHzPerT;
+ gyTpm[i] = interpLinearZero(waves[1], tSec) / gammaHzPerT;
+ gzTpm[i] = interpLinearZero(waves[2], tSec) / gammaHzPerT;
+ }
+ const longestTauMs = Math.max(
+ hardware.x.tau1Ms,
+ hardware.x.tau2Ms,
+ hardware.x.tau3Ms,
+ hardware.y.tau1Ms,
+ hardware.y.tau2Ms,
+ hardware.y.tau3Ms,
+ hardware.z.tau1Ms,
+ hardware.z.tau2Ms,
+ hardware.z.tau3Ms
+ );
+ const zptSec = longestTauMs * 4 / 1e3;
+ const preCount = Math.max(0, Math.round(zptSec / (4 * dtSec)));
+ const postCount = Math.max(0, Math.round(zptSec / dtSec));
+ const gxPadded = padSamples(gxTpm, preCount, postCount);
+ const gyPadded = padSamples(gyTpm, preCount, postCount);
+ const gzPadded = padSamples(gzTpm, preCount, postCount);
+ const stimX = safePnsModel(diff(gxPadded, dtSec), dtSec, hardware.x);
+ const stimY = safePnsModel(diff(gyPadded, dtSec), dtSec, hardware.y);
+ const stimZ = safePnsModel(diff(gzPadded, dtSec), dtSec, hardware.z);
+ const hasAnyNonTrap = blocks.some((block) => block.gx?.type === "arb" || block.gy?.type === "arb" || block.gz?.type === "arb");
+ const hasAnyLabelExt = blocks.some((block) => !!(block.labelSets?.length || block.labelIncs?.length));
+ const shift = hasAnyNonTrap || hasAnyLabelExt ? 1 : 0;
+ const selectedX = [];
+ const selectedY = [];
+ const selectedZ = [];
+ const selectedT = [];
+ for (let origIdx = 0; origIdx < nSamples; origIdx++) {
+ const paddedIdx = preCount + origIdx;
+ let stimIdx = paddedIdx - shift;
+ if (shift > 0 && hasAnyLabelExt && origIdx === tAxis.length - 1) {
+ stimIdx = Math.min(paddedIdx, stimX.length - 1);
+ }
+ if (stimIdx < 0 || stimIdx >= stimX.length || stimIdx >= stimY.length || stimIdx >= stimZ.length) continue;
+ selectedX.push(stimX[stimIdx]);
+ selectedY.push(stimY[stimIdx]);
+ selectedZ.push(stimZ[stimIdx]);
+ selectedT.push(tAxis[origIdx]);
+ }
+ const timeSec = new Float64Array(selectedX.length);
+ const pnsX = new Float64Array(selectedX.length);
+ const pnsY = new Float64Array(selectedX.length);
+ const pnsZ = new Float64Array(selectedX.length);
+ const pnsNorm = new Float64Array(selectedX.length);
+ let ok = true;
+ for (let i = 0; i < selectedX.length; i++) {
+ const xNorm = 0.01 * selectedX[i];
+ const yNorm = 0.01 * selectedY[i];
+ const zNorm = 0.01 * selectedZ[i];
+ const norm = Math.sqrt(xNorm * xNorm + yNorm * yNorm + zNorm * zNorm);
+ timeSec[i] = selectedT[i];
+ pnsX[i] = xNorm;
+ pnsY[i] = yNorm;
+ pnsZ[i] = zNorm;
+ pnsNorm[i] = norm;
+ if (norm >= 1) ok = false;
+ }
+ return { valid: true, ok, timeSec, pnsX, pnsY, pnsZ, pnsNorm };
+ }
+ function safePnsModel(dgdt, dtSec, hw) {
+ const absDgdt = new Float64Array(dgdt.length);
+ for (let i = 0; i < dgdt.length; i++) absDgdt[i] = Math.abs(dgdt[i]);
+ const dtMs = dtSec * 1e3;
+ const lp1 = lowpassTau(dgdt, hw.tau1Ms, dtMs);
+ const lp2 = lowpassTau(absDgdt, hw.tau2Ms, dtMs);
+ const lp3 = lowpassTau(dgdt, hw.tau3Ms, dtMs);
+ const stim = new Float64Array(dgdt.length);
+ const denom = hw.stimLimit > 0 ? hw.stimLimit : 1;
+ for (let i = 0; i < dgdt.length; i++) {
+ const s1 = hw.a1 * Math.abs(lp1[i]);
+ const s2 = hw.a2 * lp2[i];
+ const s3 = hw.a3 * Math.abs(lp3[i]);
+ stim[i] = (s1 + s2 + s3) / denom * hw.gScale * 100;
+ }
+ return stim;
+ }
+ function invalidPns(error) {
+ return {
+ valid: false,
+ ok: false,
+ error,
+ timeSec: new Float64Array(),
+ pnsX: new Float64Array(),
+ pnsY: new Float64Array(),
+ pnsZ: new Float64Array(),
+ pnsNorm: new Float64Array()
+ };
+ }
+ function parseAscText(text) {
+ const scalar = /* @__PURE__ */ new Map();
+ const array = /* @__PURE__ */ new Map();
+ const re = /^\s*([A-Za-z0-9_.[\]]+?)(?:\[(\d+)])?\s*=\s*([-+]?\d*\.?\d+(?:[eE][-+]?\d+)?)\s*$/;
+ for (const rawLine of text.split(/\r?\n/)) {
+ const line = rawLine.trim();
+ if (!line || line.startsWith("#") || line.startsWith("###")) continue;
+ if (/^\$include\b/i.test(line)) {
+ throw new Error("ASC contains $include directives. Use a combined ASC profile in the web viewer, or open it through the VS Code extension so companion ASC files can be resolved.");
+ }
+ const match = re.exec(line);
+ if (!match) continue;
+ const key = match[1].trim();
+ const index = match[2] === void 0 ? -1 : Number.parseInt(match[2], 10);
+ const value = Number(match[3]);
+ if (!Number.isFinite(value)) continue;
+ if (index >= 0) {
+ const values = array.get(key) ?? [];
+ values[index] = value;
+ array.set(key, values);
+ } else {
+ scalar.set(key, value);
+ }
+ }
+ return { scalar, array };
+ }
+ function resolvePnsPrefix(asc) {
+ if (asc.array.has("flGSWDTauX")) return "";
+ if (asc.array.has("GradPatSup.Phys.PNS.flGSWDTauX")) return "GradPatSup.Phys.PNS.";
+ const candidates = [...asc.array.keys()].filter((key) => key.endsWith("flGSWDTauX") && !key.toLowerCase().includes(".carns.")).sort();
+ if (candidates.length) return candidates[0].slice(0, -"flGSWDTauX".length);
+ return "GradPatSup.Phys.PNS.";
+ }
+ function getAxisHardware(asc, tauKey, aKey, stimLimitKey, stimThreshKey, gScaleKeys) {
+ const tau = findArray(asc, tauKey);
+ const weights = findArray(asc, aKey);
+ if (!tau || !weights) throw new Error(`Missing ASC arrays for ${tauKey} or ${aKey}`);
+ if (tau.length < 3 || weights.length < 3) throw new Error(`ASC arrays ${tauKey}/${aKey} require at least 3 values`);
+ const stimLimit = findScalar(asc, stimLimitKey);
+ const stimThreshold = findScalar(asc, stimThreshKey);
+ if (stimLimit === void 0 || stimThreshold === void 0) {
+ throw new Error(`Missing ASC scalar ${stimLimitKey} or ${stimThreshKey}`);
+ }
+ let gScale;
+ for (const key of gScaleKeys) {
+ gScale = findScalar(asc, key);
+ if (gScale !== void 0) break;
+ }
+ if (gScale === void 0) {
+ throw new Error("ASC is missing g_scale factors (X/Y/Z). Select a full ASC (e.g. *_twoFilesCombined.asc).");
+ }
+ return {
+ tau1Ms: tau[0],
+ tau2Ms: tau[1],
+ tau3Ms: tau[2],
+ a1: weights[0],
+ a2: weights[1],
+ a3: weights[2],
+ stimLimit,
+ stimThreshold,
+ gScale
+ };
+ }
+ function findArray(asc, key) {
+ const exact = asc.array.get(key);
+ if (exact) return exact;
+ const keyNorm = normalizeAscKey(key);
+ const chosen = [...asc.array.keys()].filter((candidate) => normalizeAscKey(candidate) === keyNorm && !candidate.toLowerCase().includes(".carns.")).sort()[0];
+ return chosen ? asc.array.get(chosen) : void 0;
+ }
+ function findScalar(asc, key) {
+ const exact = asc.scalar.get(key);
+ if (exact !== void 0) return exact;
+ const keyNorm = normalizeAscKey(key);
+ const chosen = [...asc.scalar.keys()].filter((candidate) => normalizeAscKey(candidate) === keyNorm && !candidate.toLowerCase().includes(".carns.")).sort()[0];
+ return chosen ? asc.scalar.get(chosen) : void 0;
+ }
+ function normalizeAscKey(key) {
+ return key.trim().replace(/\[\d+]/g, "");
+ }
+ function hasValidWeights(hw) {
+ return Math.abs(hw.a1 + hw.a2 + hw.a3 - 1) <= 0.01 && hw.stimLimit > 0;
+ }
+ function lowpassTau(input, tauMs, dtMs) {
+ const out = new Float64Array(input.length);
+ if (!input.length) return out;
+ if (tauMs <= 0 || dtMs <= 0) {
+ out.set(input);
+ return out;
+ }
+ const alpha = dtMs / (tauMs + dtMs);
+ out[0] = alpha * input[0];
+ for (let i = 1; i < input.length; i++) out[i] = alpha * input[i] + (1 - alpha) * out[i - 1];
+ return out;
+ }
+ function collectGradientSeries2(blocks, channel) {
+ const time = [];
+ const value = [];
+ for (const block of blocks) {
+ const grad = block[channel];
+ if (!grad?.timePoints || !grad.waveform) continue;
+ const n = Math.min(grad.timePoints.length, grad.waveform.length);
+ for (let i = 0; i < n; i++) {
+ time.push(grad.timePoints[i]);
+ value.push(grad.waveform[i]);
+ }
+ }
+ return sanitizeGradientSeries2(time, value);
+ }
+ function sanitizeGradientSeries2(time, value) {
+ const pairs = [];
+ const n = Math.min(time.length, value.length);
+ for (let i = 0; i < n; i++) {
+ if (Number.isFinite(time[i]) && Number.isFinite(value[i])) pairs.push([time[i], value[i]]);
+ }
+ pairs.sort((a, b) => a[0] - b[0]);
+ const outT = [];
+ const outV = [];
+ for (const [t, v] of pairs) {
+ const last = outT.length - 1;
+ if (last >= 0 && Math.abs(t - outT[last]) <= TIME_EPS2) {
+ outV[last] = 0.5 * (outV[last] + v);
+ continue;
+ }
+ outT.push(t);
+ outV.push(v);
+ }
+ return { time: outT, value: outV };
+ }
+ function interpLinearZero(series, t) {
+ const n = series.time.length;
+ if (!n || t < series.time[0] || t > series.time[n - 1]) return 0;
+ if (n === 1 || t <= series.time[0]) return series.value[0];
+ if (t >= series.time[n - 1]) return series.value[n - 1];
+ let lo = 0;
+ let hi = n - 1;
+ while (hi - lo > 1) {
+ const mid = lo + hi >> 1;
+ if (series.time[mid] <= t) lo = mid;
+ else hi = mid;
+ }
+ const t0 = series.time[lo];
+ const t1 = series.time[hi];
+ if (!(t1 > t0)) return series.value[lo];
+ const alpha = (t - t0) / (t1 - t0);
+ return series.value[lo] + alpha * (series.value[hi] - series.value[lo]);
+ }
+ function padSamples(input, preCount, postCount) {
+ const out = new Float64Array(preCount + input.length + postCount);
+ out.set(input, preCount);
+ return out;
+ }
+ function diff(input, dtSec) {
+ const out = new Float64Array(Math.max(0, input.length - 1));
+ for (let i = 0; i < out.length; i++) out[i] = (input[i + 1] - input[i]) / dtSec;
+ return out;
+ }
+
+ // src/pulseq/trdetect.ts
+ var GAMMA_HZ_T2 = 42576e3;
+ var DEFAULT_B0_T2 = 3;
+ function detectSequenceTiming(seq) {
+ const b0 = getB02(seq);
+ const supportsRfUse = seq.versionCombined >= 1005e3;
+ let teTimeSec = 0;
+ let hasExplicitTE = false;
+ const teDef = seq.definitions.get("EchoTime") ?? seq.definitions.get("TE");
+ if (teDef && teDef.length > 0) {
+ teTimeSec = teDef[0];
+ hasExplicitTE = true;
+ }
+ let trTimeSec = 0;
+ let hasExplicitTR = false;
+ const trDef = seq.definitions.get("RepetitionTime") ?? seq.definitions.get("TR");
+ if (trDef && trDef.length > 0) {
+ trTimeSec = trDef[0];
+ hasExplicitTR = true;
+ }
+ const rfUsePerBlock = [];
+ const excitationTimesSec = [];
+ let rfUseGuessed = false;
+ const blockStartTimes = computeCumulativeTimes(seq);
+ for (let i = 0; i < seq.blocks.length; i++) {
+ const blk = seq.blocks[i];
+ if (blk.rfId <= 0) {
+ rfUsePerBlock.push(0);
+ continue;
+ }
+ const rf = seq.rfs.get(blk.rfId);
+ if (!rf) {
+ rfUsePerBlock.push(0);
+ continue;
+ }
+ const useChar = classifyRfUse2(rf, seq, supportsRfUse, b0);
+ const useCode = useChar.charCodeAt(0);
+ rfUsePerBlock.push(useCode);
+ if (useChar === "e") {
+ const center = rf.center >= 0 ? rf.center * 1e-6 : estimateRfCenter(rf, seq);
+ const excTime = blockStartTimes[i] + rf.delay * 1e-6 + center;
+ excitationTimesSec.push(excTime);
+ }
+ if (!supportsRfUse && useChar !== "u") rfUseGuessed = true;
+ }
+ let trCount = 0;
+ const trStartBlocks = [];
+ if (!hasExplicitTR && excitationTimesSec.length >= 2) {
+ trTimeSec = estimateTRFromExcitations(excitationTimesSec);
+ hasExplicitTR = false;
+ }
+ if (trTimeSec > 0) {
+ const totalDuration = blockStartTimes.length > 0 ? blockStartTimes[blockStartTimes.length - 1] + blockDurationSeconds2(seq, seq.blocks[seq.blocks.length - 1]) : 0;
+ trCount = Math.max(1, Math.ceil(totalDuration / trTimeSec));
+ const tol = trTimeSec * 0.3;
+ let trIdx = 0;
+ for (let i = 0; i < seq.blocks.length; i++) {
+ const blkStart = blockStartTimes[i];
+ const expected = trIdx * trTimeSec;
+ if (blkStart >= expected - tol && trIdx < trCount) {
+ trStartBlocks.push(i);
+ trIdx++;
+ }
+ }
+ trStartBlocks.push(seq.blocks.length);
+ trCount = trStartBlocks.length - 1;
+ } else {
+ trCount = 0;
+ for (let i = 0; i < seq.blocks.length; i++) {
+ if (seq.blocks[i].adcId > 0) {
+ trStartBlocks.push(i);
+ trCount++;
+ }
+ }
+ trStartBlocks.push(seq.blocks.length);
+ }
+ return {
+ teTimeSec,
+ hasExplicitTE,
+ trTimeSec,
+ hasExplicitTR,
+ trCount,
+ trStartBlocks,
+ excitationTimesSec,
+ rfUseGuessed,
+ rfUsePerBlock
+ };
+ }
+ function getB02(seq) {
+ const raw = seq.definitions.get("B0") ?? seq.definitions.get("b0") ?? seq.definitions.get("b_0");
+ if (raw && Array.isArray(raw) && raw.length > 0) return +raw[0];
+ return DEFAULT_B0_T2;
+ }
+ function computeCumulativeTimes(seq) {
+ const times = [];
+ let cum = 0;
+ for (const blk of seq.blocks) {
+ times.push(cum);
+ cum += blockDurationSeconds2(seq, blk);
+ }
+ return times;
+ }
+ function blockDurationSeconds2(seq, block) {
+ if (seq.versionCombined < VER_PRE_14) return block.dur * 1e-6;
+ return block.dur * seq.rasterTimes.blockDurationRaster;
+ }
+ function classifyRfUse2(rf, seq, supportsMetadata, b0Tesla) {
+ if (supportsMetadata && rf.use && rf.use !== "u" && rf.use !== "U") {
+ return rf.use.toLowerCase();
+ }
+ const faDeg = estimateFlipAngleDeg(rf, seq);
+ if (faDeg < 90.01) return "e";
+ const freqPPM = rf.freqPPM !== 0 ? rf.freqPPM : b0Tesla > 0 ? 1e6 * rf.freqOffset / (GAMMA_HZ_T2 * b0Tesla) : 0;
+ const durEst = estimateRfDuration(rf, seq);
+ if (durEst > 6e-3 && freqPPM >= -4.5 && freqPPM <= -3) return "s";
+ return "r";
+ }
+ function estimateFlipAngleDeg(rf, seq) {
+ const magShape = seq.shapes.get(rf.magShapeId);
+ if (magShape && magShape.numSamples > 0) {
+ const raster = seq.rasterTimes.rfRaster;
+ const timeShape = rf.timeShapeId > 0 ? seq.shapes.get(rf.timeShapeId)?.samples : void 0;
+ let area = 0;
+ let prevT = timeShape ? timeShape[0] * raster : 0.5 * raster;
+ let prevAmp = Math.abs(rf.amplitude * magShape.samples[0]);
+ for (let i = 1; i < magShape.numSamples; i++) {
+ const t = timeShape ? timeShape[i] * raster : (i + 0.5) * raster;
+ const amp = Math.abs(rf.amplitude * magShape.samples[i]);
+ const dt = t - prevT;
+ if (dt > 0) area += 0.5 * (prevAmp + amp) * dt;
+ prevT = t;
+ prevAmp = amp;
+ }
+ return 360 * area;
+ }
+ const absAmp = Math.abs(rf.amplitude);
+ if (absAmp > 3e3) return 180;
+ if (absAmp > 1500) return 120;
+ return 90;
+ }
+ function estimateRfCenter(rf, _seq) {
+ const magShape = _seq.shapes.get(rf.magShapeId);
+ if (!magShape || magShape.numSamples <= 0) return 0;
+ let peakIdx = 0;
+ let peak = Math.abs(magShape.samples[0]);
+ for (let i = 1; i < magShape.numSamples; i++) {
+ const v = Math.abs(magShape.samples[i]);
+ if (v > peak) {
+ peak = v;
+ peakIdx = i;
+ }
+ }
+ const raster = _seq.rasterTimes.rfRaster;
+ const timeShape = rf.timeShapeId > 0 ? _seq.shapes.get(rf.timeShapeId)?.samples : void 0;
+ return timeShape ? (timeShape[peakIdx] ?? 0) * raster : (peakIdx + 0.5) * raster;
+ }
+ function estimateRfDuration(rf, seq) {
+ const magShape = seq.shapes.get(rf.magShapeId);
+ if (!magShape || magShape.numSamples <= 0) return 0;
+ const raster = seq.rasterTimes.rfRaster;
+ const timeShape = rf.timeShapeId > 0 ? seq.shapes.get(rf.timeShapeId)?.samples : void 0;
+ if (timeShape && timeShape.length > 0) {
+ return timeShape[timeShape.length - 1] * raster + raster;
+ }
+ return magShape.numSamples * raster;
+ }
+ function estimateTRFromExcitations(excTimesSec) {
+ if (excTimesSec.length < 2) return 0;
+ const intervals = [];
+ for (let i = 1; i < excTimesSec.length; i++) {
+ const dt = excTimesSec[i] - excTimesSec[i - 1];
+ if (dt > 1e-9) intervals.push(dt);
+ }
+ if (intervals.length === 0) return 0;
+ intervals.sort((a, b) => a - b);
+ const median = intervals[Math.floor(intervals.length / 2)];
+ const niceMs = niceRound(median * 1e3, 10);
+ return niceMs * 1e-3;
+ }
+ function niceRound(value, base) {
+ return Math.round(value / base) * base;
+ }
+
+ // src/pulseq/kspaceExportArtifacts.ts
+ function exportKspaceArtifacts(sequenceText, sequenceName, options = {}) {
+ const seq = parseSequenceText(sequenceText);
+ const decoded = decodeAllBlocks(seq);
+ const totalDuration = getTotalDuration(seq);
+ const gradientSupport = options.gradientSupport ?? "all";
+ const kspace = calculateKspace(
+ decoded,
+ seq.rasterTimes.gradientRaster,
+ totalDuration,
+ 0,
+ { maxGridPoints: options.maxGridPoints, rfRaster: seq.rasterTimes.rfRaster, gradientSupport }
+ );
+ if (!kspace) {
+ throw new Error("Unable to calculate k-space trajectory for sequence");
+ }
+ const metadata = createMetadata(
+ seq,
+ kspace,
+ sequenceName,
+ options.sequenceSha256 ?? "unknown",
+ options.packageVersion ?? "unknown",
+ !!options.includeFullTrajectory,
+ totalDuration,
+ gradientSupport
+ );
+ return {
+ ktrajAdcText: formatTrajectoryText(kspace.ktraj_adc),
+ ktrajText: options.includeFullTrajectory ? formatTrajectoryText(kspace.ktraj) : void 0,
+ metadata
+ };
+ }
+ function formatTrajectoryText(series) {
+ assertThreeEqualLengthSeries(series);
+ const n = series[0].length;
+ if (n === 0) return "";
+ const rows = [];
+ for (let i = 0; i < n; i++) {
+ rows.push(`${formatFloat(series[0][i])} ${formatFloat(series[1][i])} ${formatFloat(series[2][i])}`);
+ }
+ return `${rows.join("\n")}
+`;
+ }
+ function formatFloat(value) {
+ if (Number.isNaN(value)) return "NaN";
+ if (!Number.isFinite(value)) return value > 0 ? "Infinity" : "-Infinity";
+ const normalized = Object.is(value, -0) ? 0 : value;
+ return normalized.toExponential(12).replace(/e([+-])(\d+)$/, (_match, sign, exponent) => `e${sign}${exponent.padStart(2, "0")}`);
+ }
+ function createMetadata(seq, kspace, sequenceName, sequenceSha256, packageVersion, includeFullTrajectory, totalDurationSec, gradientSupport) {
+ return {
+ schemaVersion: 1,
+ sequenceName,
+ sequenceSha256,
+ packageVersion,
+ pulseqVersion: {
+ major: seq.version.major,
+ minor: seq.version.minor,
+ revision: seq.version.revision,
+ combined: seq.versionCombined
+ },
+ blockCount: seq.blocks.length,
+ rasterTimes: {
+ blockDuration: seq.rasterTimes.blockDurationRaster,
+ gradient: seq.rasterTimes.gradientRaster,
+ rf: seq.rasterTimes.rfRaster,
+ adc: seq.rasterTimes.adcRaster
+ },
+ totalDurationSec,
+ adcSampleCount: kspace.t_adc.length,
+ trajectorySampleCount: kspace.t_ktraj.length,
+ units: {
+ trajectory: "1/m",
+ time: "s",
+ gradient: "Hz/m",
+ convention: "Pulseq gradient integral without 2*pi factor"
+ },
+ calculation: {
+ gradientSupport
+ },
+ files: includeFullTrajectory ? { ktrajAdc: "ktraj_adc.txt", ktraj: "ktraj.txt" } : { ktrajAdc: "ktraj_adc.txt" }
+ };
+ }
+ function assertThreeEqualLengthSeries(series) {
+ if (series.length !== 3) {
+ throw new Error(`Expected three trajectory axes, received ${series.length}`);
+ }
+ const n = series[0].length;
+ if (series[1].length !== n || series[2].length !== n) {
+ throw new Error("Trajectory axes have mismatched sample counts");
+ }
+ }
+
+ // web/pulseq-browser.ts
+ var PACKAGE_VERSION = version;
+ return __toCommonJS(pulseq_browser_exports);
+})();
diff --git a/python/src/seqeyes/resources/viewer.html b/python/src/seqeyes/resources/viewer.html
new file mode 100644
index 0000000..3d24061
--- /dev/null
+++ b/python/src/seqeyes/resources/viewer.html
@@ -0,0 +1,1589 @@
+
+
+
+
+
+SeqEyes β Pulseq MRI Sequence Viewer
+
+
+
+
+
+
+
+
+
+
+
+
+
Theme:
+
+
Time:
+
Grad:
+
+
+
+
+
β hover for time
+
+
+
+
+
+
Loading sequenceβ¦
+
+
0%
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/python/tests/conftest.py b/python/tests/conftest.py
new file mode 100644
index 0000000..b124138
--- /dev/null
+++ b/python/tests/conftest.py
@@ -0,0 +1,9 @@
+# SeqEyes Python β Test Configuration
+
+import pytest
+
+
+def pytest_configure(config):
+ config.addinivalue_line(
+ "markers", "slow: marks tests as slow (deselect with '-m \"not slow\"')"
+ )
diff --git a/python/tests/debug_full.html b/python/tests/debug_full.html
new file mode 100644
index 0000000..b43eac7
--- /dev/null
+++ b/python/tests/debug_full.html
@@ -0,0 +1,3530 @@
+
+
+
+
+
+SeqEyes β Pulseq MRI Sequence Viewer
+
+
+
+
+
+
+
+
+
+
+
Theme:
+
+
Time:
+
Grad:
+
+
+
+
+
β hover for time
+
+
+
+
+
+
Loading sequenceβ¦
+
+
0%
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/python/tests/debug_minimal.html b/python/tests/debug_minimal.html
new file mode 100644
index 0000000..3704266
--- /dev/null
+++ b/python/tests/debug_minimal.html
@@ -0,0 +1,281 @@
+
+
+
+
+
+SeqEyes Debug β Hardcoded Sequence
+
+
+
+SeqEyes Debug | β hover for time | Blocks: 0 | TD: 0
+
+
+
+
+
+
diff --git a/python/tests/debug_viewer.html b/python/tests/debug_viewer.html
new file mode 100644
index 0000000..f59bdf6
--- /dev/null
+++ b/python/tests/debug_viewer.html
@@ -0,0 +1,3530 @@
+
+
+
+
+
+SeqEyes β Pulseq MRI Sequence Viewer
+
+
+
+
+
+
+
+
+
+
+
Theme:
+
+
Time:
+
Grad:
+
+
+
+
+
β hover for time
+
+
+
+
+
+
Loading sequenceβ¦
+
+
0%
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/python/tests/debug_viewer.py b/python/tests/debug_viewer.py
new file mode 100644
index 0000000..d55cdda
--- /dev/null
+++ b/python/tests/debug_viewer.py
@@ -0,0 +1,50 @@
+"""Debug: Generate the viewer HTML and check its integrity."""
+import os, re
+from seqeyes._renderer import _build_html
+
+# Use a known-good .seq file
+seq_path = os.path.join(os.path.dirname(__file__), '..', '..', 'test', 'seqeyes_demo_seq_files', 'writeFid.seq')
+seq_text = open(seq_path, 'r').read()
+
+html = _build_html(seq_text, theme='dark')
+path = os.path.join(os.path.dirname(__file__), 'debug_viewer.html')
+with open(path, 'w', encoding='utf-8') as f:
+ f.write(html)
+
+# Check key JS components
+checks = [
+ ('convertBlock function', 'function convertBlock(blk)'),
+ ('loadSequenceText function', 'function loadSequenceText(rawText)'),
+ ('SEQEYES_RAW_B64 assignment', 'window.SEQEYES_RAW_B64'),
+ ('Pulseq bundle: parseSequenceText', 'parseSequenceText'),
+ ('Pulseq bundle: decodeAllBlocks', 'decodeAllBlocks'),
+ ('convertBlock: index mapping', 'b.i = blk.index'),
+ ('convertBlock: startTime mapping', 'b.s = blk.startTime'),
+ ('convertBlock: RF amplitude', 'b.rf.a'),
+ ('convertBlock: RF timePoints', 'b.rf.t'),
+ ('convertBlock: Grad type mapping', 'b.gx.ty'),
+ ('convertBlock: ADC start', 'b.adc.s'),
+ ('convertBlock: Triggers', 'b.trg'),
+ ('draw function', 'function draw()'),
+ ('computeGlobalMax', 'function computeGlobalMax()'),
+ ('fit function', 'function fit()'),
+ ('drawBlocks function', 'function drawBlocks'),
+ ('viewerDrawFrame', 'viewerDrawFrame'),
+]
+for label, pattern in checks:
+ found = pattern in html
+ print(f' {"OK" if found else "MISSING"} {label}')
+
+# Check that the sequence data is base64-encoded
+import base64
+seq_b64 = base64.b64encode(seq_text.encode()).decode()
+data_in_html = seq_b64[:50] in html or seq_b64[-50:] in html
+print(f' {"OK" if data_in_html else "MISSING"} Sequence data in HTML')
+
+# Count script tags
+script_count = html.count('