From 12140bc3c3fcb84aab44e4a68d0b304cc6b921c4 Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Thu, 14 May 2026 11:36:52 -0500 Subject: [PATCH 1/7] Add Aqua.jl quality checks; move RegisterMismatch to test-only dep MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit RegisterMismatch was in [deps] but never imported in src/ — moved to [extras]. Added compat entries for all test extras (Aqua, ForwardDiff, ImageMagick, Random, RegisterCore, RegisterMismatch, Test, TestImages). All Aqua checks pass. Co-Authored-By: Claude Sonnet 4.6 --- .claude/freshen-package-status | 2 +- Project.toml | 13 +++++++++++-- test/runtests.jl | 5 +++++ 3 files changed, 17 insertions(+), 3 deletions(-) diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status index 663dd37..bdb428a 100644 --- a/.claude/freshen-package-status +++ b/.claude/freshen-package-status @@ -2,7 +2,7 @@ DONE: design review DONE: API review (plan → API_REVIEW_PLAN.md) DONE: update .gitignore DONE: format with runic -TODO: add Aqua.jl +DONE: add Aqua.jl TODO: remove deprecations TODO: add ExplicitImports.jl TODO: limit struct mutability diff --git a/Project.toml b/Project.toml index 4d714c0..3ac980a 100644 --- a/Project.toml +++ b/Project.toml @@ -9,27 +9,36 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" RegisterDeformation = "c19381b7-cf49-59d7-881c-50dfbd227eaf" -RegisterMismatch = "3c0dd727-6833-5f5d-a1e8-c0d421935c74" RegisterPenalty = "464fa2a9-b19c-5c59-8698-f58c971f971e" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] +Aqua = "0.8" +ForwardDiff = "0.10, 1" ImageCore = "0.8.1, 0.9, 0.10" +ImageMagick = "1" Interpolations = "0.15, 0.16" OffsetArrays = "0.10, 0.11, 1" ProgressMeter = "0.7, 0.8, 0.9, 1" +RegisterCore = "1" RegisterDeformation = "1" +Random = "1" RegisterMismatch = "1" RegisterPenalty = "1" StaticArrays = "0.10, 0.11, 0.12, 1" +Test = "1" +TestImages = "1" julia = "1.10" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RegisterCore = "67712758-55e7-5c3c-8e85-dda1d7758434" +RegisterMismatch = "3c0dd727-6833-5f5d-a1e8-c0d421935c74" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" [targets] -test = ["Test", "ForwardDiff", "ImageMagick", "RegisterCore", "TestImages"] +test = ["Aqua", "Test", "ForwardDiff", "ImageMagick", "Random", "RegisterCore", "RegisterMismatch", "TestImages"] diff --git a/test/runtests.jl b/test/runtests.jl index eff3749..123ae57 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,4 @@ +using Aqua using TestImages @static((Sys.islinux() || Sys.iswindows()) && using ImageMagick) # https://github.com/JuliaImages/ImageView.jl/pull/156#issuecomment-418200062 using Interpolations, RegisterMismatch, RegisterPenalty, RegisterDeformation @@ -12,6 +13,10 @@ if !isdefined(@__MODULE__, :test_hindsight) include("utils.jl") end +@testset "Aqua" begin + Aqua.test_all(RegisterHindsight) +end + @testset "1-dimensional" begin fixed = sin.(range(0, stop = 4π, length = 40)) # Create a circumstance where we should be able to get the answer almost exactly From be85eca475b1af4c29d5f0851ffcbb568bccf67b Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Thu, 14 May 2026 11:54:45 -0500 Subject: [PATCH 2/7] Add ExplicitImports.jl; make all implicit imports explicit in src/ Co-Authored-By: Claude Sonnet 4.6 --- .claude/freshen-package-status | 4 ++-- Project.toml | 8 +++++--- src/RegisterHindsight.jl | 12 ++++++++++-- test/runtests.jl | 11 +++++++++++ 4 files changed, 28 insertions(+), 7 deletions(-) diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status index bdb428a..89a5ed7 100644 --- a/.claude/freshen-package-status +++ b/.claude/freshen-package-status @@ -3,8 +3,8 @@ DONE: API review (plan → API_REVIEW_PLAN.md) DONE: update .gitignore DONE: format with runic DONE: add Aqua.jl -TODO: remove deprecations -TODO: add ExplicitImports.jl +DONE: remove deprecations +DONE: add ExplicitImports.jl TODO: limit struct mutability TODO: improve test coverage TODO: add and improve docstrings diff --git a/Project.toml b/Project.toml index 3ac980a..441d375 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RegisterHindsight" uuid = "73fda715-febb-5522-a0fc-086e91958f61" -authors = ["Tim Holy "] version = "0.3.0" +authors = ["Tim Holy "] [deps] ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" @@ -14,15 +14,16 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] Aqua = "0.8" +ExplicitImports = "1" ForwardDiff = "0.10, 1" ImageCore = "0.8.1, 0.9, 0.10" ImageMagick = "1" Interpolations = "0.15, 0.16" OffsetArrays = "0.10, 0.11, 1" ProgressMeter = "0.7, 0.8, 0.9, 1" +Random = "1" RegisterCore = "1" RegisterDeformation = "1" -Random = "1" RegisterMismatch = "1" RegisterPenalty = "1" StaticArrays = "0.10, 0.11, 0.12, 1" @@ -32,6 +33,7 @@ julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -41,4 +43,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" [targets] -test = ["Aqua", "Test", "ForwardDiff", "ImageMagick", "Random", "RegisterCore", "RegisterMismatch", "TestImages"] +test = ["Aqua", "ExplicitImports", "Test", "ForwardDiff", "ImageMagick", "Random", "RegisterCore", "RegisterMismatch", "TestImages"] diff --git a/src/RegisterHindsight.jl b/src/RegisterHindsight.jl index de1b821..42ba592 100644 --- a/src/RegisterHindsight.jl +++ b/src/RegisterHindsight.jl @@ -1,8 +1,16 @@ module RegisterHindsight -using ImageCore, StaticArrays, OffsetArrays, ProgressMeter -using Interpolations, RegisterDeformation, RegisterPenalty +using ImageCore: ImageCore, float64, gray +using StaticArrays: StaticArrays +using OffsetArrays: OffsetArrays +using ProgressMeter: ProgressMeter, ProgressUnknown +using Interpolations: Interpolations, AbstractExtrapolation, AbstractInterpolation, + BSpline, Linear, ScaledInterpolation +# Interpolations internals needed to manually evaluate weighted interpolation sums; +# workaround for Interpolations v0.15 removing WeightedArbIndex as a direct array index. using Interpolations: tcollect, itpflag, value_weights, coefficients, indextuple, weights +using RegisterDeformation: RegisterDeformation, GridDeformation, extrapolate, interpolate +using RegisterPenalty: RegisterPenalty, DeformationPenalty, penalty! # optimize! is deliberately unexported because it conflicts with RegisterOptimize.optimize! diff --git a/test/runtests.jl b/test/runtests.jl index 123ae57..bf19749 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using Aqua +using ExplicitImports using TestImages @static((Sys.islinux() || Sys.iswindows()) && using ImageMagick) # https://github.com/JuliaImages/ImageView.jl/pull/156#issuecomment-418200062 using Interpolations, RegisterMismatch, RegisterPenalty, RegisterDeformation @@ -17,6 +18,16 @@ end Aqua.test_all(RegisterHindsight) end +@testset "ExplicitImports" begin + # Interpolations internals (tcollect, itpflag, etc.) and Base.Slice / Interpolations.dimension_wis + # are intentional non-public accesses required by the WeightedArbIndex workaround. + ExplicitImports.test_explicit_imports( + RegisterHindsight; + ignore = (:tcollect, :itpflag, :value_weights, :coefficients, :indextuple, + :weights, :Slice, :dimension_wis), + ) +end + @testset "1-dimensional" begin fixed = sin.(range(0, stop = 4π, length = 40)) # Create a circumstance where we should be able to get the answer almost exactly From dfeb3b8a301f24800a2c48ff3fedb1ed122e3347 Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Thu, 14 May 2026 12:30:00 -0500 Subject: [PATCH 3/7] =?UTF-8?q?Add=20:gradient=20to=20ExplicitImports=20ig?= =?UTF-8?q?nore=20list=20=E2=80=94=20not=20public=20in=20Interpolations=20?= =?UTF-8?q?v0.16+?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Sonnet 4.6 --- .claude/freshen-package-status | 4 ++-- test/runtests.jl | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status index 89a5ed7..74c4427 100644 --- a/.claude/freshen-package-status +++ b/.claude/freshen-package-status @@ -5,7 +5,7 @@ DONE: format with runic DONE: add Aqua.jl DONE: remove deprecations DONE: add ExplicitImports.jl -TODO: limit struct mutability -TODO: improve test coverage +DONE: limit struct mutability +DONE: improve test coverage TODO: add and improve docstrings TODO: add or improve documentation diff --git a/test/runtests.jl b/test/runtests.jl index bf19749..4c46221 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,10 +21,11 @@ end @testset "ExplicitImports" begin # Interpolations internals (tcollect, itpflag, etc.) and Base.Slice / Interpolations.dimension_wis # are intentional non-public accesses required by the WeightedArbIndex workaround. + # gradient is not public in Interpolations v0.16+ but has no public replacement. ExplicitImports.test_explicit_imports( RegisterHindsight; ignore = (:tcollect, :itpflag, :value_weights, :coefficients, :indextuple, - :weights, :Slice, :dimension_wis), + :weights, :Slice, :dimension_wis, :gradient), ) end From 6072c86c083133fc25a049daa4da2de9b5a174ee Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Thu, 14 May 2026 13:16:38 -0500 Subject: [PATCH 4/7] Add docstrings for all public symbols; add Documenter doctests MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Module docstring for RegisterHindsight explaining purpose and entry points - InterpolatingDeformation: documents ScaledInterpolation constraint and InPlace BC requirement (previously a hidden runtime error) - All penalty_hindsight* functions: new docstrings covering both single-ϕ and dual-ϕ overloads, gradient (!) variants, and equal-footing semantics - optimize!: corrected stepsize default (1 → 1.0), documented all three moving overloads, added # Returns section and a jldoctest example - prepare_value_axes: replaced broken Interpolations v0.14-era indexing example with accurate description of current internals - Add Documenter as test dependency; run doctest(RegisterHindsight; manual=false) in runtests.jl Co-Authored-By: Claude Sonnet 4.6 --- .claude/freshen-package-status | 2 +- Project.toml | 4 +- src/RegisterHindsight.jl | 190 ++++++++++++++++++++++++++++++--- test/runtests.jl | 6 ++ 4 files changed, 185 insertions(+), 17 deletions(-) diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status index 74c4427..78b560c 100644 --- a/.claude/freshen-package-status +++ b/.claude/freshen-package-status @@ -7,5 +7,5 @@ DONE: remove deprecations DONE: add ExplicitImports.jl DONE: limit struct mutability DONE: improve test coverage -TODO: add and improve docstrings +DONE: add and improve docstrings TODO: add or improve documentation diff --git a/Project.toml b/Project.toml index 441d375..2a7a40c 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] Aqua = "0.8" +Documenter = "1" ExplicitImports = "1" ForwardDiff = "0.10, 1" ImageCore = "0.8.1, 0.9, 0.10" @@ -33,6 +34,7 @@ julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" @@ -43,4 +45,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" [targets] -test = ["Aqua", "ExplicitImports", "Test", "ForwardDiff", "ImageMagick", "Random", "RegisterCore", "RegisterMismatch", "TestImages"] +test = ["Aqua", "Documenter", "ExplicitImports", "Test", "ForwardDiff", "ImageMagick", "Random", "RegisterCore", "RegisterMismatch", "TestImages"] diff --git a/src/RegisterHindsight.jl b/src/RegisterHindsight.jl index 42ba592..ce00fc6 100644 --- a/src/RegisterHindsight.jl +++ b/src/RegisterHindsight.jl @@ -1,3 +1,23 @@ +""" + RegisterHindsight + +Hindsight-based image registration refinement. + +This module provides a gradient-descent optimizer that refines a deformation +field `ϕ` by minimizing the mean-square intensity error between a `fixed` image +and a warped `moving` image, subject to a deformation penalty. + +The primary entry point is [`optimize!`](@ref RegisterHindsight.optimize!). The +lower-level penalty functions (`penalty_hindsight`, `penalty_hindsight_data`, +`penalty_hindsight_reg`, and their in-place gradient variants) are also +accessible as `RegisterHindsight.penalty_hindsight` etc. + +The deformation `ϕ` must be an [`InterpolatingDeformation`](@ref), i.e., a +`GridDeformation` whose displacement field is backed by a `ScaledInterpolation` +with `InPlace` boundary conditions. Construct one with `interpolate!(copy(ϕ0))` +from a plain `GridDeformation` `ϕ0` (note: `interpolate(ϕ0)` produces a +different, incompatible type). +""" module RegisterHindsight using ImageCore: ImageCore, float64, gray @@ -14,6 +34,24 @@ using RegisterPenalty: RegisterPenalty, DeformationPenalty, penalty! # optimize! is deliberately unexported because it conflicts with RegisterOptimize.optimize! +""" + InterpolatingDeformation{T, N, A <: ScaledInterpolation} + +Type alias for a `GridDeformation` whose displacement field is backed by a +`ScaledInterpolation`. All functions in this module operate on deformations of +this type. + +Construct one with `interpolate!(copy(ϕ0))` from a plain `GridDeformation` +`ϕ0`. Note that `interpolate(ϕ0)` (without `!`) produces a different, +incompatible type. + +!!! note + The regularization penalty (`penalty_hindsight_reg!`) applies to the + interpolation coefficient array and therefore requires `InPlace` boundary + conditions so that the coefficient array has no padding. If the deformation + was not constructed with `InPlace` boundary conditions, a runtime error is + thrown. +""" const InterpolatingDeformation{T, N, A <: ScaledInterpolation} = GridDeformation{T, N, A} # Interpolations v0.15 no longer accepts WeightedArbIndex as a direct array index. @@ -30,10 +68,32 @@ function _coefs_at(coefs, wI) return result end +""" + penalty_hindsight(ϕ, dp, fixed, moving) → scalar + +Return the total hindsight penalty for deformation `ϕ`: the sum of the data +penalty (mean-square intensity error) and the regularization penalty from `dp`. + +`moving` must be an `AbstractInterpolation`. + +See also [`penalty_hindsight_data`](@ref RegisterHindsight.penalty_hindsight_data), +[`penalty_hindsight_reg`](@ref RegisterHindsight.penalty_hindsight_reg), +[`penalty_hindsight!`](@ref RegisterHindsight.penalty_hindsight!). +""" function penalty_hindsight(ϕ::InterpolatingDeformation, dp::DeformationPenalty, fixed, moving) return penalty_hindsight_reg(ϕ, dp) + penalty_hindsight_data(ϕ, fixed, moving) end +""" + penalty_hindsight!(g, ϕ, dp, fixed, moving) → scalar + +Compute the total hindsight penalty for `ϕ` and write its gradient to `g`. +Returns the same scalar value as +[`penalty_hindsight`](@ref RegisterHindsight.penalty_hindsight)`(ϕ, dp, fixed, moving)`. + +`g` must be allocated as `similar(ϕ.u.itp.coefs)` and is written in-place +(replacing any previous contents). +""" function penalty_hindsight!(g, ϕ::InterpolatingDeformation, dp::DeformationPenalty, fixed, moving) gd = similar(g) ret = penalty_hindsight_reg!(g, ϕ, dp) + penalty_hindsight_data!(gd, ϕ, fixed, moving) @@ -43,7 +103,16 @@ function penalty_hindsight!(g, ϕ::InterpolatingDeformation, dp::DeformationPena return ret end -# For comparison of two deformations +""" + penalty_hindsight(ϕ1, ϕ2, dp, fixed, moving) → (scalar, scalar) + +Return the total hindsight penalties for two deformations `ϕ1` and `ϕ2` as a +tuple `(p1, p2)`, evaluated on the same set of voxels. A voxel contributes to +both penalties only if both `ϕ1(x)` and `ϕ2(x)` map to in-bounds (finite) +positions in `moving`, ensuring the two values are directly comparable. + +`ϕ1` and `ϕ2` must have the same `axes` and `nodes`. +""" function penalty_hindsight( ϕ1::InterpolatingDeformation, ϕ2::InterpolatingDeformation, @@ -55,10 +124,33 @@ function penalty_hindsight( return rp1 + dp1, rp2 + dp2 end +""" + penalty_hindsight_reg(ϕ, dp) → scalar + +Return the regularization component of the hindsight penalty for deformation +`ϕ`, computed by `RegisterPenalty.penalty!` applied to the interpolation +coefficient array of `ϕ`. + +!!! note + Requires `InPlace` boundary conditions. See [`InterpolatingDeformation`](@ref). +""" function penalty_hindsight_reg(ϕ::InterpolatingDeformation, dp) return penalty_hindsight_reg!(nothing, ϕ, dp) end +""" + penalty_hindsight_reg!(g, ϕ, dp) → scalar + +Compute the regularization penalty for `ϕ` and write its gradient to `g`. +Returns the same scalar value as +[`penalty_hindsight_reg`](@ref RegisterHindsight.penalty_hindsight_reg)`(ϕ, dp)`. + +`g` must be `nothing` (to skip gradient computation) or allocated as +`similar(ϕ.u.itp.coefs)` and is written in-place. + +Throws an error if the deformation does not use `InPlace` boundary conditions +(i.e., if `axes(ϕ.u.itp.coefs) ≠ ϕ.u.itp.parentaxes`). +""" function penalty_hindsight_reg!(g, ϕ::InterpolatingDeformation, dp) # The regularization penalty. We apply this to the interpolation # coefficients rather than the on-grid values. This may be @@ -68,6 +160,16 @@ function penalty_hindsight_reg!(g, ϕ::InterpolatingDeformation, dp) return penalty!(g, dp, itp.coefs) end +""" + penalty_hindsight_data(ϕ, fixed, moving) → scalar + +Return the data penalty for deformation `ϕ`: the mean-square intensity error +between `fixed` and `moving` evaluated at the deformed positions `ϕ(x)`, +averaged over all voxels where both `fixed[x]` and `moving(ϕ(x))` are finite. + +`moving` must be an `AbstractInterpolation` (typically also an extrapolation +returning `NaN` for out-of-bounds coordinates). +""" function penalty_hindsight_data( ϕ::InterpolatingDeformation{T, N, A}, fixed::AbstractArray{T1, N}, @@ -93,6 +195,20 @@ function penalty_hindsight_data( return mm / valid end +""" + penalty_hindsight_data!(g, ϕ, fixed, moving) → scalar + +Compute the data penalty for `ϕ` and write its gradient with respect to the +interpolation coefficients of `ϕ` to `g`. Returns the same scalar value as +[`penalty_hindsight_data`](@ref RegisterHindsight.penalty_hindsight_data)`(ϕ, fixed, moving)`. + +`g` must be allocated as `similar(ϕ.u.itp.coefs)` and is zeroed and written +in-place. Its element type must be an `SVector` matching the spatial +dimensionality of `ϕ` (e.g., `SVector{2, Float64}` for 2-D images). + +`moving` must be an `AbstractInterpolation` (not merely an extrapolation) so +that `Interpolations.gradient` can be evaluated on it. +""" # This re-uses the work of computing the weights for both the value and the gradient function penalty_hindsight_data!( g, @@ -141,12 +257,17 @@ function penalty_hindsight_data!( end """ - coefs, wI = prepare_value_axes(ϕ::InterpolatingDeformation) + coefs, wis = prepare_value_axes(ϕ::InterpolatingDeformation) + +Return the interpolation coefficient array `coefs` and a tuple `wis` of +per-dimension `WeightedIndex` arrays for evaluating `ϕ` at each position in +its node grid. `wis[d][i]` is the `WeightedIndex` for the `i`-th position +along dimension `d`. -Return the coefficients and `Interpolations.WeightedIndex` values needed for evaluating -`ϕ` at each position in the range grid. `coefs` is an array and `wI` a tuple, `wI[d]` -corresponding to axis `d` of the range grid. In particular, at location `I = (I1, ..., In)` -the value is `coefs[wI[1][I1], ..., wI[n][In]]`. +This is a low-level helper used internally by +[`penalty_hindsight_data`](@ref RegisterHindsight.penalty_hindsight_data) and +[`penalty_hindsight_data!`](@ref RegisterHindsight.penalty_hindsight_data!). +The coefficient-weighted sum at a given node position is computed by `_coefs_at`. """ function prepare_value_axes(ϕ::InterpolatingDeformation) itp = ϕ.u.itp @@ -157,6 +278,16 @@ function prepare_value_axes(ϕ::InterpolatingDeformation) return coefficients(itp), wis end +""" + penalty_hindsight_data(ϕ1, ϕ2, fixed, moving) → (scalar, scalar) + +Return the data penalties for two deformations `ϕ1` and `ϕ2` as a tuple +`(p1, p2)`, evaluated on the same set of voxels. A voxel contributes to both +penalties only if `moving` returns a finite value at both `ϕ1(x)` and +`ϕ2(x)`, ensuring the two values are directly comparable. + +`ϕ1` and `ϕ2` must have the same `nodes`. +""" # This implementation allows comparing ϕ1 and ϕ2 on equal footing, meaning using the same # voxels of the image data. A voxel gets included only if both ϕ1 and ϕ2 are in-bounds. # This cannot be guaranteed for the single-ϕ implementation of penalty_hindsight_data. @@ -193,18 +324,47 @@ function penalty_hindsight_data( end """ - result = optimize!(ϕ, dp, fixed, moving; stepsize=1, itermax=1000) + result = optimize!(ϕ, dp, fixed, moving; stepsize=1.0, itermax=1000) + +Minimize the mean-square error between `fixed[x...]` and `moving(ϕ(x)...)` by +updating `ϕ` in place via gradient descent on the interpolation coefficients. +`dp` is a deformation penalty (e.g., `AffinePenalty`). `stepsize` is the +maximum per-iteration update in pixels applied to elements of `ϕ`. + +`moving` may be: +- an `AbstractExtrapolation` (used directly), +- an `AbstractInterpolation` (wrapped automatically with `NaN` extrapolation), or +- a plain `AbstractArray` (interpolated with `BSpline(Linear())` then extrapolated). + +Terminates when a gradient-descent step would increase the penalty value, or +after `itermax` iterations. + +# Returns + +A named tuple `(; final, initial)` where `final` is the penalty after +optimization and `initial` is the penalty before optimization. Tuple +destructuring also works: `final, initial = optimize!(...)`. + +# Example + +```jldoctest +julia> using RegisterDeformation, RegisterPenalty, Interpolations + +julia> fixed = sin.(range(0, stop=4π, length=40)); + +julia> nodes = (range(1, stop=40, length=5),); + +julia> ϕ = interpolate!(GridDeformation(zeros(1, 5), nodes)); + +julia> ap = AffinePenalty(ϕ.nodes, 0.01); -Improve `ϕ` by minimizing the mean-square error between `fixed[x...]` -and `moving(ϕ(x)...)`. `dp` is a deformation penalty, e.g., -an `AffinePenalty`. `stepsize` specifies the maximum update, in pixels, of -elements of `ϕ`. +julia> moving = sin.(range(0.2, stop=4π + 0.2, length=40)); -The operation terminates when the step increases the value of the penalty, -or when more than `itermax` iterations have occurred. +julia> result = RegisterHindsight.optimize!(ϕ, ap, fixed, moving; stepsize=0.1); -Returns a named tuple `(; final, initial)` with the penalty after optimization -and the penalty before optimization. Tuple destructuring `final, initial = optimize!(...)` also works. +julia> result.final < result.initial +true +``` """ function optimize!(ϕ::InterpolatingDeformation, dp::DeformationPenalty, fixed, moving::AbstractExtrapolation; stepsize = 1.0, itermax = 1000) # Optimize the interpolation coefficients, rather than the values diff --git a/test/runtests.jl b/test/runtests.jl index 4c46221..e55fbe4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using Aqua +using Documenter using ExplicitImports using TestImages @static((Sys.islinux() || Sys.iswindows()) && using ImageMagick) # https://github.com/JuliaImages/ImageView.jl/pull/156#issuecomment-418200062 @@ -14,6 +15,11 @@ if !isdefined(@__MODULE__, :test_hindsight) include("utils.jl") end +DocMeta.setdocmeta!(RegisterHindsight, :DocTestSetup, :(using RegisterHindsight); recursive = true) +@testset "Doctests" begin + doctest(RegisterHindsight; manual = false) +end + @testset "Aqua" begin Aqua.test_all(RegisterHindsight) end From 9dc7d3e291b2b7a10a4978c6c246ada2af51c889 Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Thu, 14 May 2026 13:24:45 -0500 Subject: [PATCH 5/7] Add README with badges, installation, concept, and usage example Co-Authored-By: Claude Sonnet 4.6 --- .claude/freshen-package-status | 11 ----- README.md | 89 ++++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 11 deletions(-) delete mode 100644 .claude/freshen-package-status create mode 100644 README.md diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status deleted file mode 100644 index 78b560c..0000000 --- a/.claude/freshen-package-status +++ /dev/null @@ -1,11 +0,0 @@ -DONE: design review -DONE: API review (plan → API_REVIEW_PLAN.md) -DONE: update .gitignore -DONE: format with runic -DONE: add Aqua.jl -DONE: remove deprecations -DONE: add ExplicitImports.jl -DONE: limit struct mutability -DONE: improve test coverage -DONE: add and improve docstrings -TODO: add or improve documentation diff --git a/README.md b/README.md new file mode 100644 index 0000000..307723a --- /dev/null +++ b/README.md @@ -0,0 +1,89 @@ +# RegisterHindsight + +[![CI](https://github.com/HolyLab/RegisterHindsight.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/RegisterHindsight.jl/actions/workflows/CI.yml) +[![codecov](https://codecov.io/gh/HolyLab/RegisterHindsight.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HolyLab/RegisterHindsight.jl) + +RegisterHindsight refines image-registration deformation fields using gradient +descent. Given a `fixed` image and a `moving` image, it adjusts the +displacements of a `GridDeformation` to minimize the mean-square intensity +error subject to a deformation-smoothness penalty. + +It is part of the +[HolyLab registration pipeline](https://github.com/HolyLab/HolyLabRegistry.git) +and is designed to be called after a coarse registration step (e.g., +[RegisterOptimize](https://github.com/HolyLab/RegisterOptimize.jl)) to squeeze +out residual misalignment. + +## Installation + +RegisterHindsight is registered in the +[HolyLab registry](https://github.com/HolyLab/HolyLabRegistry). Add the +registry once, then install normally: + +```julia +using Pkg +pkg"registry add https://github.com/HolyLab/HolyLabRegistry.git" +Pkg.add("RegisterHindsight") +``` + +## Concept + +"Hindsight" refers to the optimization strategy: rather than working with the +raw displacement values on the deformation grid, the optimizer adjusts the +*interpolation coefficients* that back the deformation field. This exposes a +smooth, differentiable objective that can be descended with a simple line +search. + +The deformation `ϕ` must be an **interpolating deformation** — a +`GridDeformation` whose displacement field is backed by a +`ScaledInterpolation` with `InPlace` boundary conditions. Construct one with: + +```julia +ϕ = interpolate!(copy(ϕ0)) # note: interpolate! (with !), not interpolate +``` + +`interpolate(ϕ0)` (without `!`) produces a different, incompatible type. + +## Usage + +```julia +using RegisterDeformation, RegisterPenalty, Interpolations +using RegisterHindsight + +# 1-D example: shift a sine wave by ~0.05 radians +fixed = sin.(range(0, stop=4π, length=40)) +moving = sin.(range(0.2, stop=4π + 0.2, length=40)) + +nodes = (range(1, stop=40, length=5),) +ϕ = interpolate!(GridDeformation(zeros(1, 5), nodes)) + +ap = AffinePenalty(ϕ.nodes, 0.01) # smoothness penalty weight 0.01 + +result = RegisterHindsight.optimize!(ϕ, ap, fixed, moving; stepsize=0.1) + +result.final < result.initial # true — penalty decreased +``` + +`optimize!` returns a named tuple `(; final, initial)` with the penalty +before and after optimization, so you can check how much improvement was +achieved. + +### Penalty functions + +The penalty and its gradient can also be called directly, which is useful for +diagnostics or for building a custom optimizer: + +| Function | Description | +|---|---| +| `RegisterHindsight.penalty_hindsight(ϕ, dp, fixed, moving)` | Total penalty (data + regularization) | +| `RegisterHindsight.penalty_hindsight!(g, ϕ, dp, fixed, moving)` | Total penalty + gradient in `g` | +| `RegisterHindsight.penalty_hindsight_data(ϕ, fixed, moving)` | Data (intensity error) term only | +| `RegisterHindsight.penalty_hindsight_reg(ϕ, dp)` | Regularization term only | + +A two-deformation overload `penalty_hindsight(ϕ1, ϕ2, dp, fixed, moving)` +evaluates both deformations on the same set of valid voxels so the two +penalties are directly comparable. + +> **Note:** `optimize!` is intentionally not exported because it would +> conflict with `RegisterOptimize.optimize!`. Call it as +> `RegisterHindsight.optimize!(...)`. From 7b44b622db34d36134c935ed96beb7c154af2d5b Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Thu, 14 May 2026 13:25:22 -0500 Subject: [PATCH 6/7] Bump version to 1.0.0 Co-Authored-By: Claude Sonnet 4.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2a7a40c..7ff0ec5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "RegisterHindsight" uuid = "73fda715-febb-5522-a0fc-086e91958f61" -version = "0.3.0" +version = "1.0.0" authors = ["Tim Holy "] [deps] From 41365f80141b472d22dd0c72574b59487eb5f1f2 Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Thu, 14 May 2026 13:32:42 -0500 Subject: [PATCH 7/7] Fix ProgressUnknown deprecation warning that broke doctests Co-Authored-By: Claude Sonnet 4.6 --- src/RegisterHindsight.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RegisterHindsight.jl b/src/RegisterHindsight.jl index ce00fc6..8ca237f 100644 --- a/src/RegisterHindsight.jl +++ b/src/RegisterHindsight.jl @@ -376,7 +376,7 @@ function optimize!(ϕ::InterpolatingDeformation, dp::DeformationPenalty, fixed, ∇objective!(g, ϕ) = penalty_hindsight!(g, ϕ, dp, fixed, moving) pold = p0 = objective(ϕ) iter = 0 - prog = ProgressUnknown("Performing descent:") + prog = ProgressUnknown(; desc="Performing descent:") while iter < itermax iter += 1 ∇objective!(g, ϕ)