Skip to content

♻️ Streamline ROMEO kernels and ship MIT notice with binaries #22

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
vanandrew merged 2 commits into from
Apr 25, 2026
Merged

♻️ Streamline ROMEO kernels and ship MIT notice with binaries #22

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
634fcbc
:recycle: Streamline ROMEO weight + voxel_quality kernels
vanandrew Apr 25, 2026
7f8eddf
:wrench: Ship ROMEO MIT notice with wheels and standalone bundles
vanandrew Apr 25, 2026
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
47 changes: 10 additions & 37 deletions include/romeo/voxel_quality.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,50 +48,23 @@ inline std::vector<T> voxel_quality(const T* phase4d, std::size_t nx, std::size_

std::vector<T> qmap(vol, T(0));

// Forward-edge sum: qmap[v] = Σ_d w[d, v] for the 3 dims.
for (std::size_t idx = 0; idx < vol; ++idx) {
qmap[idx] = w[0u + 3u * idx] + w[1u + 3u * idx] + w[2u + 3u * idx];
}

// Backward-edge contributions. Julia:
// qmap[2:end,:,:] .+= weights[1, 1:end-1, :, :]
// qmap[:,2:end,:] .+= weights[2, :, 1:end-1, :]
// qmap[:,:,2:end] .+= weights[3, :, :, 1:end-1]
// In 0-based: each interior voxel adds the forward-edge weight of its
// previous neighbor along that dim.
const std::ptrdiff_t sx = 1;
const std::ptrdiff_t sy = static_cast<std::ptrdiff_t>(nx);
const std::ptrdiff_t sz = static_cast<std::ptrdiff_t>(nx * ny);

// Single fused pass: each voxel sums its three forward edges plus the
// forward-edge weight of its previous neighbor along each axis (= the
// backward-edge contribution from Julia's three .+= expressions).
// Associativity matches the original four-pass form: forward triplet
// first, then back-x, then back-y, then back-z.
for (std::size_t k = 0; k < nz; ++k) {
for (std::size_t j = 0; j < ny; ++j) {
for (std::size_t i = 1; i < nx; ++i) {
const std::size_t idx = i + nx * (j + ny * k);
const std::size_t prev = idx - static_cast<std::size_t>(sx);
qmap[idx] += w[0u + 3u * prev];
}
}
}
for (std::size_t k = 0; k < nz; ++k) {
for (std::size_t j = 1; j < ny; ++j) {
for (std::size_t i = 0; i < nx; ++i) {
const std::size_t idx = i + nx * (j + ny * k);
const std::size_t prev = idx - static_cast<std::size_t>(sy);
qmap[idx] += w[1u + 3u * prev];
T q = w[0u + 3u * idx] + w[1u + 3u * idx] + w[2u + 3u * idx];
if (i > 0) q += w[0u + 3u * (idx - 1)];
if (j > 0) q += w[1u + 3u * (idx - nx)];
if (k > 0) q += w[2u + 3u * (idx - nx * ny)];
qmap[idx] = q / T(6);
}
}
}
for (std::size_t k = 1; k < nz; ++k) {
for (std::size_t j = 0; j < ny; ++j) {
for (std::size_t i = 0; i < nx; ++i) {
const std::size_t idx = i + nx * (j + ny * k);
const std::size_t prev = idx - static_cast<std::size_t>(sz);
qmap[idx] += w[2u + 3u * prev];
}
}
}

for (std::size_t i = 0; i < vol; ++i) qmap[i] /= T(6);
return qmap;
}

Expand Down
124 changes: 49 additions & 75 deletions include/romeo/weights.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,91 +119,65 @@ inline std::vector<OutT> calculate_weights_romeo_impl(const T* phase, std::size_
if (phase2 == nullptr || TEs == nullptr) flags.phase_gradient_coherence = false;

const std::size_t n = nx * ny * nz;
const std::ptrdiff_t stride[3] = {1, static_cast<std::ptrdiff_t>(nx), static_cast<std::ptrdiff_t>(nx * ny)};
const std::ptrdiff_t signed_n = static_cast<std::ptrdiff_t>(n);

// ROMEO's `parsekwargs` multiplies mag by mask when both are present. We
// do the same pre-pass so `maxmag` is computed over the masked region.
std::vector<T> masked_mag;
const T* mag_in_use = mag;
if (mag != nullptr && mask != nullptr) {
masked_mag.assign(mag, mag + n);
for (std::size_t i = 0; i < n; ++i) {
if (!mask[i]) masked_mag[i] = T(0);
}
mag_in_use = masked_mag.data();
}

// maxmag: `maximum(mag[isfinite.(mag)])` — infinities and NaNs excluded.
T maxmag = T(0);
bool has_finite_mag = false;
if (mag_in_use != nullptr) {
for (std::size_t i = 0; i < n; ++i) {
T v = mag_in_use[i];
if (std::isfinite(v) && (!has_finite_mag || v > maxmag)) {
maxmag = v;
has_finite_mag = true;
}
}
}
(void)maxmag; // only needed once magweight/magweight2 (flags 5-6) are ported.
(void)has_finite_mag;
const std::ptrdiff_t sx = 1;
const std::ptrdiff_t sy = static_cast<std::ptrdiff_t>(nx);
const std::ptrdiff_t sz = static_cast<std::ptrdiff_t>(nx * ny);

std::vector<OutT> weights(3 * n, zero_value);

for (int dim = 0; dim < 3; ++dim) {
const std::ptrdiff_t step = stride[dim];
for (std::ptrdiff_t idx = 0; idx < signed_n; ++idx) {
const std::ptrdiff_t jdx = idx + step;
if (jdx >= signed_n) continue;
if (mask != nullptr && !mask[idx]) continue;
// Mirrors ROMEO's `mag .* mask` pre-pass without materializing a copy: any
// masked-out voxel reads as 0 magnitude.
auto mag_at = [&](std::ptrdiff_t i) -> T {
if (mask != nullptr && !mask[i]) return T(0);
return mag[i];
};

T w = T(1);
if (flags.phase_coherence) {
w *= T(0.1) + T(0.9) * detail::phase_coherence(phase[idx], phase[jdx]);
}
if (flags.phase_gradient_coherence) {
w *= T(0.1) +
T(0.9) * detail::phase_gradient_coherence(phase[idx], phase[jdx], phase2[idx], phase2[jdx],
TEs[0], TEs[1]);
}
if (flags.phase_linearity) {
w *= T(0.1) + T(0.9) * detail::phase_linearity(phase, idx, jdx, step, signed_n);
}
if (mag_in_use != nullptr) {
T mi = mag_in_use[idx];
T mj = mag_in_use[jdx];
T small = std::min(mi, mj);
T big = std::max(mi, mj);
if (flags.mag_coherence) {
w *= T(0.1) + T(0.9) * detail::mag_coherence(small, big);
}
// flags 5-6 (magweight, magweight2) not ported — warpkit uses :romeo, not :romeo6.
auto edge_weight = [&](std::ptrdiff_t idx, std::ptrdiff_t jdx, std::ptrdiff_t stride) -> T {
T w = T(1);
if (flags.phase_coherence) {
w *= T(0.1) + T(0.9) * detail::phase_coherence(phase[idx], phase[jdx]);
}
if (flags.phase_gradient_coherence) {
w *= T(0.1) + T(0.9) * detail::phase_gradient_coherence(phase[idx], phase[jdx], phase2[idx], phase2[jdx],
TEs[0], TEs[1]);
}
if (flags.phase_linearity) {
w *= T(0.1) + T(0.9) * detail::phase_linearity(phase, idx, jdx, stride, signed_n);
}
if (mag != nullptr) {
T mi = mag_at(idx);
T mj = mag_at(jdx);
T small = std::min(mi, mj);
T big = std::max(mi, mj);
if (flags.mag_coherence) {
w *= T(0.1) + T(0.9) * detail::mag_coherence(small, big);
}

// Column-major (3, nx, ny, nz): index = dim + 3*idx
weights[static_cast<std::size_t>(dim) + 3u * static_cast<std::size_t>(idx)] = rescale_fn(w);
// flags 5-6 (magweight, magweight2) not ported — warpkit uses :romeo, not :romeo6.
}
}
return w;
};

// Zero out the boundary edges that don't actually exist (Julia:
// weights[1,end,:,:] .= 0 / [2,:,end,:] .= 0 / [3,:,:,end] .= 0).
// One pass over voxels in (k, j, i) order. Boundary edges are skipped by
// construction via the `+1 < n*` guards, so the trailing zero-fill pass
// from the Julia source is unnecessary — the value-init of `weights`
// already leaves those slots at `zero_value`.
for (std::size_t k = 0; k < nz; ++k) {
for (std::size_t j = 0; j < ny; ++j) {
std::size_t idx = (nx - 1) + nx * (j + ny * k);
weights[0u + 3u * idx] = zero_value;
}
}
for (std::size_t k = 0; k < nz; ++k) {
for (std::size_t i = 0; i < nx; ++i) {
std::size_t idx = i + nx * ((ny - 1) + ny * k);
weights[1u + 3u * idx] = zero_value;
}
}
for (std::size_t j = 0; j < ny; ++j) {
for (std::size_t i = 0; i < nx; ++i) {
std::size_t idx = i + nx * (j + ny * (nz - 1));
weights[2u + 3u * idx] = zero_value;
for (std::size_t i = 0; i < nx; ++i) {
const std::size_t idx = i + nx * (j + ny * k);
if (mask != nullptr && !mask[idx]) continue;
const std::ptrdiff_t sidx = static_cast<std::ptrdiff_t>(idx);
if (i + 1 < nx) {
weights[0u + 3u * idx] = rescale_fn(edge_weight(sidx, sidx + sx, sx));
}
if (j + 1 < ny) {
weights[1u + 3u * idx] = rescale_fn(edge_weight(sidx, sidx + sy, sy));
}
if (k + 1 < nz) {
weights[2u + 3u * idx] = rescale_fn(edge_weight(sidx, sidx + sz, sz));
}
}
}
}

Expand Down
15 changes: 15 additions & 0 deletions packaging/pyinstaller/build_bundle.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,20 @@ def write_readme(bundle: Path, version: str, target: str) -> None:
(bundle / "README.md").write_text(body, encoding="utf-8")


def write_licenses(bundle: Path) -> None:
# Bundle warpkit's own LICENSE plus the upstream ROMEO MIT notice. The
# ROMEO notice is required: the wk-* binaries link in compiled code derived
# from include/romeo/, which the MIT terms require we ship the notice with.
repo_root = Path(__file__).resolve().parents[2]
licenses_dir = bundle / "LICENSES"
licenses_dir.mkdir(exist_ok=True)
shutil.copy2(repo_root / "LICENSE", licenses_dir / "warpkit.LICENSE")
shutil.copy2(
repo_root / "include" / "romeo" / "LICENSE",
licenses_dir / "ROMEO.LICENSE",
)


def make_zip(bundle: Path, out_zip: Path) -> None:
# shutil.make_archive's base_dir keeps a tidy top-level folder inside the zip.
base_name = str(out_zip.with_suffix(""))
Expand Down Expand Up @@ -144,6 +158,7 @@ def main() -> int:

adhoc_sign_macos(versioned)
write_readme(versioned, version, target)
write_licenses(versioned)

out_zip = dist / f"warpkit-{version}-{target}.zip"
make_zip(versioned, out_zip)
Expand Down
6 changes: 6 additions & 0 deletions packaging/pyinstaller/bundle_README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,3 +47,9 @@ xattr -r -d com.apple.quarantine /opt/warpkit-@VERSION@
- `wk-compute-jacobian` — compute the Jacobian determinant of a warp

Run any of them with `--help` for usage.

## Licenses

`LICENSES/warpkit.LICENSE` is the warpkit license. `LICENSES/ROMEO.LICENSE` is
the upstream MIT notice for the ROMEO phase-unwrapping algorithms compiled
into these binaries.
4 changes: 4 additions & 0 deletions pyproject.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,10 @@ authors = [{ name = "Andrew Van", email = "vanandrew77@gmail.com" }]
keywords = ["neuroimaging"]
classifiers = ["Programming Language :: Python :: 3"]
urls = { github = "https://github.com/vanandrew/warpkit" }
# include/romeo/LICENSE is the upstream MIT notice for the ROMEO C++ port; it
# travels with the headers in source distributions and (via PEP 639
# license-files) gets bundled into wheels under *.dist-info/licenses/.
license-files = ["LICENSE", "include/romeo/LICENSE"]
dynamic = ["version"]
dependencies = [
"nibabel >= 4.0.2",
Expand Down
Loading