Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 22 additions & 8 deletions CLAUDE.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,16 @@ Pre-print: <https://www.biorxiv.org/content/10.1101/2023.11.28.568744v1>.
## Layout

- `warpkit/` — Python package. Public entry points: `warpkit.distortion.medic`,
`warpkit.utilities.*`. CLI scripts in `warpkit/scripts/` are auto-discovered
by `setup.py` and exposed as `medic` and `extract_field_from_maps`.
`warpkit.utilities.*`. CLI scripts live in `warpkit/scripts/` and are
registered explicitly in `[project.scripts]` in `pyproject.toml`. All CLIs
ship with a `wk-` prefix to avoid colliding with same-named tools from
FSL/ANTs/AFNI/etc.: `wk-medic`, `wk-unwrap-phase`, `wk-compute-fieldmap`,
`wk-apply-warp`, `wk-convert-warp`, `wk-convert-fieldmap`,
`wk-compute-jacobian`. Adding a new CLI means adding a new file under
`warpkit/scripts/` *and* a new line to `[project.scripts]` — there is
no longer any auto-discovery. Shared IO helpers used by the
`wk-convert-*` and `wk-compute-jacobian` scripts live in
`warpkit/scripts/_warp_io.py` (private; not a CLI).
- `warpkit/warpkit_cpp.pyi` + `warpkit/py.typed` — type info for the compiled
extension, shipped via `MANIFEST.in` and `[tool.setuptools.package-data]`.
Regenerate the stub after pybind11 binding changes via the wrapper script —
Expand All @@ -40,6 +48,9 @@ uv sync --group dev --config-setting editable_mode=strict
# tests
uv run pytest -q

# coverage (matches the CI `coverage` job)
uv run coverage run && uv run coverage report -m

# lint + types (matches pre-commit; never bypass with --no-verify)
uvx ruff check
uvx ruff format
Expand Down Expand Up @@ -106,9 +117,12 @@ and call out anything CI-relevant (wheel matrix, pybind11 ABI, ITK).

## CI specifics

GitHub Actions builds wheels for Python 3.11–3.14 on `ubuntu-latest` and
`macos-latest` via cibuildwheel. `pyproject.toml`'s `[tool.cibuildwheel]`
skips `*musllinux*` and free-threaded builds (`cp313t-*`, `cp314t-*`);
re-enabling free-threaded support requires auditing the pybind11 + ITK code
paths for the no-GIL ABI. PyPI publish and the GHCR Docker image only run on
a published GitHub release.
GitHub Actions builds wheels for Python 3.11–3.14 on `ubuntu-latest`,
`ubuntu-24.04-arm`, and `macos-latest` via cibuildwheel.
`pyproject.toml`'s `[tool.cibuildwheel]` skips `*musllinux*` and the
`cp314t-*` free-threaded build; re-enabling free-threaded support requires
auditing the pybind11 + ITK code paths for the no-GIL ABI. The sdist job
also runs `uv run coverage run` then `coverage report -m` — keep coverage
healthy when adding code (the `[tool.coverage.report]` config in
`pyproject.toml` omits the test files). PyPI publish and the GHCR Docker
image only run on a published GitHub release.
6 changes: 3 additions & 3 deletions Dockerfile
Original file line number Diff line number Diff line change
Expand Up @@ -21,11 +21,11 @@ ENV UV_PROJECT_ENVIRONMENT=/opt/warpkit/.venv
ENV UV_PYTHON_PREFERENCE=only-managed
RUN cd /opt/warpkit && uv sync --group dev --config-setting editable_mode=strict -v

# put the project venv on PATH so `medic` and friends resolve
# put the project venv on PATH so `wk-medic` and friends resolve
ENV PATH=/opt/warpkit/.venv/bin:${PATH}

# test warpkit
RUN cd /opt/warpkit && uv run pytest -s -v

# set medic script as entrypoint
ENTRYPOINT ["medic"]
# set wk-medic script as entrypoint
ENTRYPOINT ["wk-medic"]
71 changes: 66 additions & 5 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ Pre-built wheels are published for Linux (x86_64) and macOS (universal2). If `pi
docker run -it --rm ghcr.io/vanandrew/warpkit:latest --help
```

The image's entrypoint is the `medic` CLI.
The image's entrypoint is the `wk-medic` CLI.

### From source

Expand Down Expand Up @@ -99,12 +99,28 @@ displacement_field = displacement_map_to_field(

### CLI

A `medic` script is installed on `PATH`. Acquisition parameters can come either from BIDS sidecars or from the command line — pick one.
All warpkit CLIs are installed on `PATH` with a `wk-` prefix to avoid colliding
with same-named tools from FSL/ANTs/AFNI/etc.:

| Command | Purpose |
| ---------------------- | ------------------------------------------------------------------------------------- |
| `wk-medic` | End-to-end MEDIC pipeline: phase + magnitude → field maps + displacement maps. |
| `wk-unwrap-phase` | Stage 1: ROMEO multi-echo phase unwrapping → unwrapped phase + per-frame masks. |
| `wk-compute-fieldmap` | Stage 2: take stage-1 outputs → native + displacement + undistorted-space field maps. |
| `wk-apply-warp` | Resample an image through a displacement map / field (single or per-frame series). |
| `wk-convert-warp` | Convert between maps ↔ fields and between ITK / FSL / ANTs / AFNI; `--invert` warps. |
| `wk-convert-fieldmap` | Convert between mm displacement maps/fields and Hz B0 field maps. |
| `wk-compute-jacobian` | Per-voxel Jacobian determinant (1 = no change, <1 = compression, >1 = expansion). |

`wk-medic` runs the full pipeline; `wk-unwrap-phase` + `wk-compute-fieldmap`
run the same thing in two stages so you can inspect/reuse the unwrapped
phase. Acquisition parameters can come either from BIDS sidecars or from the
command line — pick one.

From BIDS sidecars:

```bash
medic \
wk-medic \
--magnitude mag_e1.nii.gz mag_e2.nii.gz mag_e3.nii.gz \
--phase phase_e1.nii.gz phase_e2.nii.gz phase_e3.nii.gz \
--metadata mag_e1.json mag_e2.json mag_e3.json \
Expand All @@ -116,7 +132,7 @@ medic \
Or by passing acquisition parameters directly:

```bash
medic \
wk-medic \
--magnitude mag_e1.nii.gz mag_e2.nii.gz mag_e3.nii.gz \
--phase phase_e1.nii.gz phase_e2.nii.gz phase_e3.nii.gz \
--TEs 14.2 38.93 63.66 \
Expand All @@ -127,7 +143,52 @@ medic \

`--TEs` is in **milliseconds**, `--total-readout-time` in **seconds**, and `--phase-encoding-direction` is one of `i, j, k, i-, j-, k-, x, y, z, x-, y-, z-`.

Run `medic --help` for the full option list (noise-frame trimming, CPU count, debug mode, etc.).
Run any `wk-*` CLI with `--help` for the full option list.

#### Common follow-on workflows

Apply a MEDIC displacement-map series to the BOLD it was derived from
(per-frame distortion correction):

```bash
wk-apply-warp \
--input bold.nii.gz \
--transform sub-01_run-01_displacementmaps.nii \
--phase-encoding-axis j \
--output bold_corrected.nii.gz
```

Convert MEDIC's per-frame displacement maps into per-frame ANTs-format
displacement fields:

```bash
wk-convert-warp \
--input sub-01_run-01_displacementmaps.nii \
--to field --axis j --to-format ants \
--output field_{0..14}.nii.gz
```

Compute the per-frame Jacobian determinant (volume-change map) of those
displacement maps:

```bash
wk-compute-jacobian \
--input sub-01_run-01_displacementmaps.nii \
--axis j \
--output sub-01_run-01_jacobian.nii
```

Convert MEDIC's mm displacement maps to a Hz B0 field map (or back —
this CLI handles either direction, with maps or fields on the mm side):

```bash
wk-convert-fieldmap \
--input sub-01_run-01_displacementmaps.nii \
--to fieldmap \
--total-readout-time 0.0501 \
--phase-encoding-direction j- \
--output sub-01_run-01_fieldmap.nii
```

## Authors

Expand Down
119 changes: 40 additions & 79 deletions include/romeo/romeo.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,67 +20,28 @@ namespace py = pybind11;

namespace romeo {

// Python-facing facade over the pure-C++ ROMEO implementation.
// Python-facing entry points over the pure-C++ ROMEO implementation.
//
// Method names retain the `romeo_*` prefix from the original Julia-backed
// pybind class so call sites in warpkit/unwrap.py read naturally.
template <typename T>
class Romeo {
public:
Romeo() = default;
~Romeo() = default;

// Port of ROMEO.jl `calculateweights(phase; weights=:romeo, ...)` — the
// only weight preset we support. Exposed primarily so Python tests can
// validate the internal machinery against literal goldens from
// ROMEO.jl test/specialcases.jl. Not used by warpkit itself.
//
// `phase` is a column-major (nx, ny, nz) array. `mag`, `phase2`, `mask` may
// be 0-sized / empty arrays to indicate "not provided"; `TEs` is required
// only when `phase2` is provided (length 2: [te_phase, te_phase2]).
// Returns a (3, nx, ny, nz) uint8 array.
py::array_t<std::uint8_t, py::array::f_style> calculate_weights(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> mag,
py::array_t<T, py::array::f_style> phase2,
py::array_t<T, py::array::f_style> TEs,
py::array_t<bool, py::array::f_style> mask);

py::array_t<T, py::array::f_style> romeo_voxelquality(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> TEs,
py::array_t<T, py::array::f_style> mag);

py::array_t<T, py::array::f_style> romeo_unwrap3D(py::array_t<T, py::array::f_style> phase,
std::string weights,
py::array_t<T, py::array::f_style> mag,
py::array_t<bool, py::array::f_style> mask,
bool correctglobal = false,
int maxseeds = 1,
bool merge_regions = false,
bool correct_regions = false);

py::array_t<T, py::array::f_style> romeo_unwrap4D(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> TEs,
std::string weights,
py::array_t<T, py::array::f_style> mag,
py::array_t<bool, py::array::f_style> mask,
bool correctglobal = false,
int maxseeds = 1,
bool merge_regions = false,
bool correct_regions = false);
};

// ----------------------------------------------------------------------------
// Method implementations (Phase 1: all throw; wired to the public entry points
// under include/romeo/*.h so Phase 2+ can fill them in without touching this
// file's signatures).
// ----------------------------------------------------------------------------

// Function names retain the `romeo_*` prefix from the original Julia-backed
// API so call sites in warpkit/unwrap.py read naturally. They are stateless
// free functions; pybind11 binds them as module-level functions in
// `warpkit_cpp` (no class wrapper).

// Port of ROMEO.jl `calculateweights(phase; weights=:romeo, ...)` — the
// only weight preset we support. Exposed primarily so Python tests can
// validate the internal machinery against literal goldens from
// ROMEO.jl test/specialcases.jl. Not used by warpkit itself.
//
// `phase` is a column-major (nx, ny, nz) array. `mag`, `phase2`, `mask` may
// be 0-sized / empty arrays to indicate "not provided"; `TEs` is required
// only when `phase2` is provided (length 2: [te_phase, te_phase2]).
// Returns a (3, nx, ny, nz) uint8 array.
template <typename T>
py::array_t<std::uint8_t, py::array::f_style> Romeo<T>::calculate_weights(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> mag,
py::array_t<T, py::array::f_style> phase2,
py::array_t<T, py::array::f_style> TEs,
py::array_t<bool, py::array::f_style> mask) {
py::array_t<std::uint8_t, py::array::f_style> calculate_weights(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> mag,
py::array_t<T, py::array::f_style> phase2,
py::array_t<T, py::array::f_style> TEs,
py::array_t<bool, py::array::f_style> mask) {
if (phase.ndim() != 3) throw std::invalid_argument("calculate_weights: phase must be 3D");
const auto nx = static_cast<std::size_t>(phase.shape(0));
const auto ny = static_cast<std::size_t>(phase.shape(1));
Expand Down Expand Up @@ -118,9 +79,9 @@ py::array_t<std::uint8_t, py::array::f_style> Romeo<T>::calculate_weights(py::ar
}

template <typename T>
py::array_t<T, py::array::f_style> Romeo<T>::romeo_voxelquality(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> TEs,
py::array_t<T, py::array::f_style> mag) {
py::array_t<T, py::array::f_style> romeo_voxelquality(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> TEs,
py::array_t<T, py::array::f_style> mag) {
if (phase.ndim() != 4) throw std::invalid_argument("romeo_voxelquality: phase must be 4D");
const auto nx = static_cast<std::size_t>(phase.shape(0));
const auto ny = static_cast<std::size_t>(phase.shape(1));
Expand Down Expand Up @@ -150,14 +111,14 @@ py::array_t<T, py::array::f_style> Romeo<T>::romeo_voxelquality(py::array_t<T, p
}

template <typename T>
py::array_t<T, py::array::f_style> Romeo<T>::romeo_unwrap3D(py::array_t<T, py::array::f_style> phase,
std::string weights,
py::array_t<T, py::array::f_style> mag,
py::array_t<bool, py::array::f_style> mask,
bool correctglobal,
int maxseeds,
bool merge_regions,
bool correct_regions) {
py::array_t<T, py::array::f_style> romeo_unwrap3D(py::array_t<T, py::array::f_style> phase,
std::string weights,
py::array_t<T, py::array::f_style> mag,
py::array_t<bool, py::array::f_style> mask,
bool correctglobal = false,
int maxseeds = 1,
bool merge_regions = false,
bool correct_regions = false) {
if (weights != "romeo")
throw std::invalid_argument("romeo_unwrap3D: only the \"romeo\" weight preset is supported.");
if (merge_regions || correct_regions)
Expand Down Expand Up @@ -198,15 +159,15 @@ py::array_t<T, py::array::f_style> Romeo<T>::romeo_unwrap3D(py::array_t<T, py::a
}

template <typename T>
py::array_t<T, py::array::f_style> Romeo<T>::romeo_unwrap4D(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> TEs,
std::string weights,
py::array_t<T, py::array::f_style> mag,
py::array_t<bool, py::array::f_style> mask,
bool correctglobal,
int maxseeds,
bool merge_regions,
bool correct_regions) {
py::array_t<T, py::array::f_style> romeo_unwrap4D(py::array_t<T, py::array::f_style> phase,
py::array_t<T, py::array::f_style> TEs,
std::string weights,
py::array_t<T, py::array::f_style> mag,
py::array_t<bool, py::array::f_style> mask,
bool correctglobal = false,
int maxseeds = 1,
bool merge_regions = false,
bool correct_regions = false) {
if (weights != "romeo")
throw std::invalid_argument("romeo_unwrap4D: only the \"romeo\" weight preset is supported.");
if (merge_regions || correct_regions)
Expand Down
11 changes: 10 additions & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ authors = [{ name = "Andrew Van", email = "vanandrew77@gmail.com" }]
keywords = ["neuroimaging"]
classifiers = ["Programming Language :: Python :: 3"]
urls = { github = "https://github.com/vanandrew/warpkit" }
dynamic = ["version", "scripts"]
dynamic = ["version"]
dependencies = [
"nibabel >= 4.0.2",
"numpy >= 1.23.3",
Expand All @@ -17,6 +17,15 @@ dependencies = [
"indexed-gzip >= 1.7.0",
]

[project.scripts]
"wk-medic" = "warpkit.scripts.medic:main"
"wk-unwrap-phase" = "warpkit.scripts.unwrap_phase:main"
"wk-compute-fieldmap" = "warpkit.scripts.compute_fieldmap:main"
"wk-apply-warp" = "warpkit.scripts.apply_warp:main"
"wk-convert-warp" = "warpkit.scripts.convert_warp:main"
"wk-convert-fieldmap" = "warpkit.scripts.convert_fieldmap:main"
"wk-compute-jacobian" = "warpkit.scripts.compute_jacobian:main"

[dependency-groups]
dev = [
"coverage[toml] >= 5.5",
Expand Down
8 changes: 0 additions & 8 deletions setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@

THISDIR = Path(__file__).parent
cmake_dir = (Path(THISDIR)).absolute().as_posix()
scripts_path = THISDIR / "warpkit" / "scripts"

IS_CIBUILDWHEEL = os.environ.get("CIBUILDWHEEL", "0") == "1"
IS_MACOS = os.environ.get("RUNNER_OS", "0") == "macOS"
Expand All @@ -32,11 +31,4 @@
)
],
cmdclass={"build_ext": BuildExtension},
entry_points={
"console_scripts": [
f"{f.stem}=warpkit.scripts.{f.stem}:main"
for f in scripts_path.glob("*.py")
if f.name not in "__init__.py"
]
},
)
Loading