From 634fcbc0b9358856be767550b78b0512535fcfa2 Mon Sep 17 00:00:00 2001 From: Andrew Van Date: Sat, 25 Apr 2026 17:02:27 -0500 Subject: [PATCH 1/2] :recycle: Streamline ROMEO weight + voxel_quality kernels Refactor calculate_weights_romeo_impl into a single (k, j, i) pass over voxels: skip boundary edges by construction (i+1 --- include/romeo/voxel_quality.h | 47 +++---------- include/romeo/weights.h | 124 ++++++++++++++-------------------- 2 files changed, 59 insertions(+), 112 deletions(-) diff --git a/include/romeo/voxel_quality.h b/include/romeo/voxel_quality.h index 7533352..95534ce 100644 --- a/include/romeo/voxel_quality.h +++ b/include/romeo/voxel_quality.h @@ -48,50 +48,23 @@ inline std::vector voxel_quality(const T* phase4d, std::size_t nx, std::size_ std::vector 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(nx); - const std::ptrdiff_t sz = static_cast(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(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(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(sz); - qmap[idx] += w[2u + 3u * prev]; - } - } - } - - for (std::size_t i = 0; i < vol; ++i) qmap[i] /= T(6); return qmap; } diff --git a/include/romeo/weights.h b/include/romeo/weights.h index 1a87773..b68f283 100644 --- a/include/romeo/weights.h +++ b/include/romeo/weights.h @@ -119,91 +119,65 @@ inline std::vector 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(nx), static_cast(nx * ny)}; const std::ptrdiff_t signed_n = static_cast(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 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(nx); + const std::ptrdiff_t sz = static_cast(nx * ny); std::vector 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(dim) + 3u * static_cast(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(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)); + } + } } } From 7f8eddf42ab82bf86174dc208321d467ad1b3464 Mon Sep 17 00:00:00 2001 From: Andrew Van Date: Sat, 25 Apr 2026 17:02:36 -0500 Subject: [PATCH 2/2] :wrench: Ship ROMEO MIT notice with wheels and standalone bundles The ROMEO C++ port under include/romeo/ is a derivative of MIT-licensed code, so the notice has to travel with binary distributions, not just the source tree. - pyproject.toml: declare license-files = [LICENSE, include/romeo/LICENSE] (PEP 639). Verified via `uv build --sdist` that PKG-INFO emits both files as License-File entries, so wheels also carry them under *.dist-info/licenses/. - packaging/pyinstaller/build_bundle.py: drop both files into a LICENSES/ directory inside each standalone bundle. - bundle_README.md: point users at LICENSES/. The Dockerfile already copies the full source tree (ADD . /opt/warpkit) so both notices ship in the image untouched. Co-Authored-By: Claude Opus 4.7 (1M context) --- packaging/pyinstaller/build_bundle.py | 15 +++++++++++++++ packaging/pyinstaller/bundle_README.md | 6 ++++++ pyproject.toml | 4 ++++ 3 files changed, 25 insertions(+) diff --git a/packaging/pyinstaller/build_bundle.py b/packaging/pyinstaller/build_bundle.py index 0e877a0..f79c2c1 100644 --- a/packaging/pyinstaller/build_bundle.py +++ b/packaging/pyinstaller/build_bundle.py @@ -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("")) @@ -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) diff --git a/packaging/pyinstaller/bundle_README.md b/packaging/pyinstaller/bundle_README.md index 0be6e8d..c7bb27c 100644 --- a/packaging/pyinstaller/bundle_README.md +++ b/packaging/pyinstaller/bundle_README.md @@ -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. diff --git a/pyproject.toml b/pyproject.toml index 20e7ac5..b538789 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -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",