From 3eb11540ac00e3aa2c699083c91f0350c4605232 Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Tue, 19 May 2026 14:35:43 -0500
Subject: [PATCH 01/11] Add Aqua.jl quality checks

- Remove stale RegisterMismatch from [deps] (only appeared in a warning string)
- Add missing compat entries for stdlib and test-only deps
- Add Aqua testset with treat_as_own for the intentional RegisterCore.maxshift extension

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 Project.toml              | 15 ++++++++++++---
 test/register_optimize.jl |  5 +++++
 2 files changed, 17 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 51cda05..4e0e593 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "RegisterOptimize"
 uuid = "b981a1b5-4587-53b0-9c3d-454c648f6f8d"
-authors = ["Tim Holy <tim.holy@gmail.com>"]
 version = "1.0.0"
+authors = ["Tim Holy <tim.holy@gmail.com>"]
 
 [deps]
 CachedInterpolations = "b9709bfb-d23d-5560-8276-8c35c4b76823"
@@ -21,33 +21,42 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 RegisterCore = "67712758-55e7-5c3c-8e85-dda1d7758434"
 RegisterDeformation = "c19381b7-cf49-59d7-881c-50dfbd227eaf"
 RegisterFit = "36121b08-3789-5198-aff2-59a3443d9b59"
-RegisterMismatch = "3c0dd727-6833-5f5d-a1e8-c0d421935c74"
 RegisterPenalty = "464fa2a9-b19c-5c59-8698-f58c971f971e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
+Aqua = "0.8"
 CachedInterpolations = "1"
 CenterIndexedArrays = "1"
 CoordinateTransformations = "0.5, 0.6"
 ForwardDiff = "0.10, 1"
 Interpolations = "0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
 Ipopt = "0.5, 0.6, 1"
+ImageMagick = "1"
+Images = "0.26"
 IterativeSolvers = "0.7, 0.8, 0.9"
 JuMP = "1.2"
+LinearAlgebra = "1"
 MathOptInterface = "1"
 NLsolve = "3, 4"
 Optim = "0.18, 0.19, 0.20, 1"
+Printf = "1"
 ProgressMeter = "0.7, 0.8, 0.9, 1"
 RegisterCore = "1"
 RegisterDeformation = "1"
 RegisterFit = "1"
 RegisterPenalty = "1"
+RegisterUtilities = "1"
+Rotations = "1"
 StaticArrays = "0.10, 0.11, 0.12, 1"
 Statistics = "1"
+Test = "1"
+TestImages = "1"
 julia = "1.10"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 RegisterUtilities = "d4862ba2-f42c-5aeb-af4f-96a8884a16c4"
@@ -56,4 +65,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["Test", "Rotations", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
+test = ["Aqua", "Test", "Rotations", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
diff --git a/test/register_optimize.jl b/test/register_optimize.jl
index 3aed28f..879c02b 100644
--- a/test/register_optimize.jl
+++ b/test/register_optimize.jl
@@ -4,6 +4,11 @@ import RegisterOptimize
 using RegisterCore, RegisterPenalty, RegisterDeformation, RegisterMismatch, RegisterFit
 using Images, CoordinateTransformations, Rotations, RegisterOptimize, LinearAlgebra
 using RegisterUtilities
+using Aqua
+
+@testset "Aqua" begin
+    Aqua.test_all(RegisterOptimize; piracies=(; treat_as_own=[RegisterCore.maxshift]))
+end
 
 ###
 ### Global-optimum initial guess

From 65acecfb717b7f6b0f7e9a3c908f54021e864f13 Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Tue, 19 May 2026 14:40:50 -0500
Subject: [PATCH 02/11] Add ExplicitImports.jl and make all module imports
 explicit

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 Project.toml              |  8 +++++---
 src/RegisterOptimize.jl   | 27 ++++++++++++++++++++++-----
 test/register_optimize.jl |  6 ++++++
 3 files changed, 33 insertions(+), 8 deletions(-)

diff --git a/Project.toml b/Project.toml
index 4e0e593..75b009b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,11 +30,12 @@ Aqua = "0.8"
 CachedInterpolations = "1"
 CenterIndexedArrays = "1"
 CoordinateTransformations = "0.5, 0.6"
+ExplicitImports = "1"
 ForwardDiff = "0.10, 1"
-Interpolations = "0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
-Ipopt = "0.5, 0.6, 1"
 ImageMagick = "1"
 Images = "0.26"
+Interpolations = "0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
+Ipopt = "0.5, 0.6, 1"
 IterativeSolvers = "0.7, 0.8, 0.9"
 JuMP = "1.2"
 LinearAlgebra = "1"
@@ -57,6 +58,7 @@ julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 RegisterUtilities = "d4862ba2-f42c-5aeb-af4f-96a8884a16c4"
@@ -65,4 +67,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["Aqua", "Test", "Rotations", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
+test = ["Aqua", "ExplicitImports", "Test", "Rotations", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
diff --git a/src/RegisterOptimize.jl b/src/RegisterOptimize.jl
index 92cc61a..dedf4f3 100644
--- a/src/RegisterOptimize.jl
+++ b/src/RegisterOptimize.jl
@@ -1,11 +1,28 @@
 module RegisterOptimize
 
 import MathOptInterface as MOI
-using JuMP: JuMP, Model, optimizer_with_attributes, @variable, @objective, @operator, @constraint, termination_status, LOCALLY_SOLVED
-using Ipopt, Optim, Interpolations, ForwardDiff, StaticArrays, IterativeSolvers, ProgressMeter
-using RegisterCore, RegisterDeformation, RegisterPenalty, RegisterFit, CachedInterpolations, CenterIndexedArrays
-using Printf, LinearAlgebra, Statistics, CoordinateTransformations
+using JuMP: JuMP, Model, optimizer_with_attributes, @variable, @objective, @operator, termination_status, LOCALLY_SOLVED
+using CachedInterpolations: CachedInterpolations, CachedInterpolation, cachedinterpolators
+using CenterIndexedArrays: CenterIndexedArrays, CenterIndexedArray
+using CoordinateTransformations: CoordinateTransformations, AffineMap
+using ForwardDiff: ForwardDiff
+using Interpolations: Interpolations, AbstractExtrapolation, AbstractInterpolation,
+                      BSpline, InPlace, Linear, OnCell, Quadratic
+using Ipopt: Ipopt
+using IterativeSolvers: IterativeSolvers, cg
+using LinearAlgebra: LinearAlgebra, I, dot, mul!, tr
+using Optim: Optim, Fminbox, LBFGS, OnceDifferentiable
+using Printf: Printf, @printf
+using ProgressMeter: ProgressMeter, @showprogress
+using RegisterCore: RegisterCore, ColonFun, NumDenom, maxshift, ratio
+using RegisterDeformation: RegisterDeformation, GridDeformation, extrapolate,
+                           interpolate!, rotation2, rotation3, rotationparameters,
+                           tformeye, tformrotate, tformtranslate, transform
 using RegisterDeformation: convert_to_fixed, convert_from_fixed
+using RegisterFit: RegisterFit, mms2fit, uclamp!
+using RegisterPenalty: RegisterPenalty, AffinePenalty, DeformationPenalty, penalty!
+using StaticArrays: StaticArrays, SArray, SMatrix, SVector, Size, StaticMatrix, similar_type
+using Statistics: Statistics, mean
 using Base: tail
 
 import Base: *
@@ -677,7 +694,7 @@ function optimize!(ϕs::Vector{<:GridDeformation}, ϕs_old, dp::AffinePenalty, m
     ub1 = T[mxs...] .- T(RegisterFit.register_half)
     ub = repeat(ub1, outer = [div(length(uvec), length(ub1))])
     results = Optim.optimize(df, -ub, ub, uvec, Fminbox(LBFGS()), Optim.Options(x_tol = 1.0e-4, kwargs...))
-    return _copy!(ϕs, Optim.minimizer(results)), Optim.minimum(results)
+    return _copy!(ϕs, Optim.minimizer(results)), minimum(results)
 end
 
 function optimize!(ϕs::Vector{<:GridDeformation}, ϕs_old, dp::AffinePenalty{T, N}, mmis::Array{Tf}; λt = nothing, kwargs...) where {Tf <: Number, T, N}
diff --git a/test/register_optimize.jl b/test/register_optimize.jl
index 879c02b..2c056cf 100644
--- a/test/register_optimize.jl
+++ b/test/register_optimize.jl
@@ -5,11 +5,17 @@ using RegisterCore, RegisterPenalty, RegisterDeformation, RegisterMismatch, Regi
 using Images, CoordinateTransformations, Rotations, RegisterOptimize, LinearAlgebra
 using RegisterUtilities
 using Aqua
+using ExplicitImports
 
 @testset "Aqua" begin
     Aqua.test_all(RegisterOptimize; piracies=(; treat_as_own=[RegisterCore.maxshift]))
 end
 
+@testset "ExplicitImports" begin
+    @test ExplicitImports.check_no_implicit_imports(RegisterOptimize) === nothing
+    @test ExplicitImports.check_no_stale_explicit_imports(RegisterOptimize) === nothing
+end
+
 ###
 ### Global-optimum initial guess
 ###

From 16e61d0957f867e7996eefca5a7361dfc0af6d5b Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Tue, 19 May 2026 14:44:41 -0500
Subject: [PATCH 03/11] Mark all mutable struct fields const
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

No fields are ever reassigned after construction — only their internal
contents are mutated in-place (e.g. ap.λ, ϕ.u). Adding const enforces
this invariant and helps the compiler.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 src/RegisterOptimize.jl | 44 ++++++++++++++++++++---------------------
 1 file changed, 22 insertions(+), 22 deletions(-)

diff --git a/src/RegisterOptimize.jl b/src/RegisterOptimize.jl
index dedf4f3..8f4b448 100644
--- a/src/RegisterOptimize.jl
+++ b/src/RegisterOptimize.jl
@@ -270,11 +270,11 @@ to_float(::Type{T}, A, B) where {T} = convert(Array{Float32}, A), convert(Array{
 ### Rigid registration from raw images, MathProg interface
 ###
 mutable struct RigidValue{N, A <: AbstractArray, I <: AbstractExtrapolation, SDT} <: MOI.AbstractNLPEvaluator
-    fixed::A
-    wfixed::A
-    moving::I
-    SD::SDT
-    thresh
+    const fixed::A
+    const wfixed::A
+    const moving::I
+    const SD::SDT
+    const thresh
 end
 
 function RigidValue(fixed::AbstractArray, moving::AbstractArray{T}, SD, thresh) where {T <: Real}
@@ -301,8 +301,8 @@ function (d::RigidValue)(x)
 end
 
 mutable struct RigidOpt{RV <: RigidValue, G} <: GradOnlyBoundsOnly
-    rv::RV
-    g::G
+    const rv::RV
+    const g::G
 end
 
 function RigidOpt(fixed, moving, SD, thresh)
@@ -426,9 +426,9 @@ end
 
 # A type for computing multiplication by the linear operator
 mutable struct AffineQHessian{AP <: AffinePenalty, M <: StaticMatrix, N, Φ}
-    ap::AP
-    Qs::Array{M, N}
-    ϕ_old::Φ
+    const ap::AP
+    const Qs::Array{M, N}
+    const ϕ_old::Φ
 end
 
 function AffineQHessian(ap::AffinePenalty{T}, Qs::AbstractArray{TQ, N}, ϕ_old) where {T, TQ, N}
@@ -649,10 +649,10 @@ end
 
 
 mutable struct DeformOpt{D, Dold, DP, M} <: GradOnlyBoundsOnly
-    ϕ::D
-    ϕ_old::Dold
-    dp::DP
-    mmis::M
+    const ϕ::D
+    const ϕ_old::Dold
+    const dp::DP
+    const mmis::M
 end
 # (d::DeformOpt)(x) = MOI.eval_objective(d, x)
 
@@ -705,11 +705,11 @@ function optimize!(ϕs::Vector{<:GridDeformation}, ϕs_old, dp::AffinePenalty{T,
 end
 
 mutable struct DeformTseriesOpt{D, Dsold, DP, T, M} <: GradOnlyBoundsOnly
-    ϕs::Vector{D}
-    ϕs_old::Dsold
-    dp::DP
-    λt::T
-    mmis::M
+    const ϕs::Vector{D}
+    const ϕs_old::Dsold
+    const dp::DP
+    const λt::T
+    const mmis::M
 end
 
 # Using MOI is a legacy of the old Ipopt interface, but
@@ -1107,9 +1107,9 @@ end
 
 
 mutable struct SigmoidOpt{G, H} <: BoundsOnly
-    data::Vector{Float64}
-    g::G
-    h::H
+    const data::Vector{Float64}
+    const g::G
+    const h::H
 end
 
 SigmoidOpt(data::Vector{Float64}) = SigmoidOpt(data, y -> ForwardDiff.gradient(x -> sigpenalty(x, data), y), y -> ForwardDiff.hessian(x -> sigpenalty(x, data), y))

From e0314fda4c94c789bacc10ee75ff9d136d04326f Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 09:46:38 -0500
Subject: [PATCH 04/11] Improve test coverage: rigid registration, to_full,
 fit_sigmoid, grid_rotations

Fix _optimize! dispatch for BSpline{Linear{...}} (Interpolations v0.15+ changed
the type from BSpline{Linear} to BSpline{Linear{Throw{OnGrid}}}), so the
subgradient path was silently bypassed. Add new testsets for to_float, p2rigid,
to_full, grid_rotations, fit_sigmoid/sigpenalty, and optimize_rigid. Also declare
RegisterMismatch as an explicit test dependency (it was already used but unlisted).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 Project.toml              |  4 +-
 src/RegisterOptimize.jl   |  2 +-
 test/register_optimize.jl | 86 +++++++++++++++++++++++++++++++++++++++
 3 files changed, 90 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 75b009b..20aa69e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -47,6 +47,7 @@ ProgressMeter = "0.7, 0.8, 0.9, 1"
 RegisterCore = "1"
 RegisterDeformation = "1"
 RegisterFit = "1"
+RegisterMismatch = "1"
 RegisterPenalty = "1"
 RegisterUtilities = "1"
 Rotations = "1"
@@ -61,10 +62,11 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+RegisterMismatch = "3c0dd727-6833-5f5d-a1e8-c0d421935c74"
 RegisterUtilities = "d4862ba2-f42c-5aeb-af4f-96a8884a16c4"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["Aqua", "ExplicitImports", "Test", "Rotations", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
+test = ["Aqua", "ExplicitImports", "Test", "Rotations", "RegisterMismatch", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
diff --git a/src/RegisterOptimize.jl b/src/RegisterOptimize.jl
index 8f4b448..98f7a67 100644
--- a/src/RegisterOptimize.jl
+++ b/src/RegisterOptimize.jl
@@ -540,7 +540,7 @@ end
 # algorithms. In that case we use a subgradient method, using gradient
 # descent with a "constant" step length (using an L1 measure of
 # length). See https://en.wikipedia.org/wiki/Subgradient_method.
-function _optimize!(ϕ, ϕ_old, dp::DeformationPenalty, mmis, ::Type{BSpline{Linear}}; stepsize = 1.0, kwargs...)
+function _optimize!(ϕ, ϕ_old, dp::DeformationPenalty, mmis, ::Type{<:BSpline{<:Interpolations.Linear}}; stepsize = 1.0, kwargs...)
     mxs = maxshift(first(mmis))
     g = similar(ϕ.u)
     gview = convert_from_fixed(g)
diff --git a/test/register_optimize.jl b/test/register_optimize.jl
index 2c056cf..39c003b 100644
--- a/test/register_optimize.jl
+++ b/test/register_optimize.jl
@@ -332,3 +332,89 @@ end
         @test -1.01 <= ϕ.u[i][1] <= -0.99
     end
 end
+
+@testset "to_float" begin
+    A, B = RegisterOptimize.to_float([1, 2, 3], [4, 5, 6])
+    @test eltype(A) == Float32
+    @test eltype(B) == Float32
+
+    A, B = RegisterOptimize.to_float([1.0, 2.0], [3.0, 4.0])
+    @test eltype(A) == Float64
+    @test eltype(B) == Float64
+
+    A, B = RegisterOptimize.to_float(Float32[1, 2], [3.0, 4.0])
+    @test eltype(A) == Float64
+    @test eltype(B) == Float64
+end
+
+@testset "p2rigid" begin
+    SD1 = reshape([1.0], 1, 1)
+    r1 = RegisterOptimize.p2rigid([5.0], SD1)
+    @test r1.translation ≈ [5.0]
+
+    SD = Matrix{Float64}(I, 2, 2)
+    angle = π / 6
+    r2 = RegisterOptimize.p2rigid([angle, 0.0, 0.0], SD)
+    @test r2.linear ≈ [cos(angle) -sin(angle); sin(angle) cos(angle)]
+    @test r2.translation ≈ [0.0, 0.0]
+
+    SD3 = Matrix{Float64}(I, 3, 3)
+    r3 = RegisterOptimize.p2rigid([0.0, 0.0, π / 6, 1.0, 2.0, 3.0], SD3)
+    @test r3.translation ≈ [1.0, 2.0, 3.0]
+
+    @test_throws ErrorException RegisterOptimize.p2rigid([1.0, 2.0], SD)
+end
+
+@testset "to_full" begin
+    nodes = (range(1, stop = 20, length = 4), range(1, stop = 15, length = 4))
+    ap = AffinePenalty(nodes, 1.0)
+    gridsize = map(length, nodes)
+    Qs = [Matrix{Float64}(I, 2, 2) for _ in CartesianIndices(gridsize)]
+    A = RegisterOptimize.to_full(ap, Qs)
+    @test size(A) == (32, 32)
+    @test A ≈ A'
+    @test all(>=(0), diag(A))
+end
+
+@testset "grid_rotations 2D" begin
+    SD = Matrix{Float64}(I, 2, 2)
+    rots = RegisterOptimize.grid_rotations(π / 4, 3, SD)
+    @test length(rots) == 3
+    @test rots[2].linear ≈ Matrix{Float64}(I, 2, 2)
+    @test rots[2].translation ≈ zeros(2)
+
+    # Even grid size gets rounded up to next odd integer
+    rots4 = RegisterOptimize.grid_rotations(π / 4, 4, SD)
+    @test length(rots4) == 5
+end
+
+@testset "fit_sigmoid" begin
+    n = 20
+    bot, top, ctr, wid = 1.0, 10.0, 5.0, 2.0
+    data = bot .+ (top - bot) ./ (1 .+ exp.(-(collect(1:n) .- ctr) ./ wid))
+
+    @test RegisterOptimize.sigpenalty([bot, top, ctr, wid], data) < 1e-20
+
+    b, t, c, w, fval = RegisterOptimize.fit_sigmoid(data, bot, top, ctr, wid)
+    @test isfinite(fval)
+    @test fval >= 0
+
+    b2, t2, c2, w2, fval2 = RegisterOptimize.fit_sigmoid(data)
+    @test isfinite(fval2)
+
+    @test_throws ErrorException RegisterOptimize.fit_sigmoid([1.0, 2.0, 3.0])
+    @test_throws ErrorException RegisterOptimize.fit_sigmoid([1.0, 2.0, 3.0], 1.0, 3.0, 2.0, 1.0)
+end
+
+@testset "optimize_rigid" begin
+    img = zeros(30, 30)
+    img[8:22, 8:22] .= 1.0
+    angle = 0.1
+    tfm_in = AffineMap(rotation2(angle), zeros(2))
+    img_rot = transform(img, tfm_in)
+    SD = Matrix{Float64}(I, 2, 2)
+    result_tfm, fval = RegisterOptimize.optimize_rigid(img, img_rot, tfm_in, [5, 5]; SD)
+    @test fval < 0.05
+    # recovered rotation is close to rotation2(-angle)
+    @test result_tfm.linear[1, 2] ≈ sin(angle) atol = 0.05
+end

From 47d92e7936d1e8e21040651da2550b19108b4543 Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 10:11:50 -0500
Subject: [PATCH 05/11] Rewrite all public docstrings and add Documenter
 doctests
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- Module: remove nonexistent qd_rigid, list unexported helpers separately
- auto_λ: fix 6-element return signature (was 5, silently misassigned λ_all
  to datapenalty); two-signature header; Returns: bullet list
- auto_λt: fix "fix"→"fit"; document λts/datapenalty return types
- fixed_λ: clarify mmis must be "produced by interpolate_mm!" (not "already
  interpolating"); document ϕ and penalty types
- initial_deformation: qualify "globally-optimal" as w.r.t. the quadratic
  approximation; document converged=false meaning and fallback behavior; add
  jldoctest example
- fit_sigmoid: replace misleading [optional-arg] notation with two explicit
  signatures; add return description including objective_value; add jldoctest;
  imperative mood
- Add Documenter as test dependency; run doctests via doctest(; manual=false)

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 Project.toml              |   4 +-
 src/RegisterOptimize.jl   | 213 +++++++++++++++++++++++---------------
 test/register_optimize.jl |  10 ++
 3 files changed, 145 insertions(+), 82 deletions(-)

diff --git a/Project.toml b/Project.toml
index 20aa69e..d163cf4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,6 +27,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Aqua = "0.8"
+Documenter = "1"
 CachedInterpolations = "1"
 CenterIndexedArrays = "1"
 CoordinateTransformations = "0.5, 0.6"
@@ -59,6 +60,7 @@ julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
@@ -69,4 +71,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [targets]
-test = ["Aqua", "ExplicitImports", "Test", "Rotations", "RegisterMismatch", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
+test = ["Aqua", "Documenter", "ExplicitImports", "Test", "Rotations", "RegisterMismatch", "RegisterUtilities", "ImageMagick", "Images", "TestImages"]
diff --git a/src/RegisterOptimize.jl b/src/RegisterOptimize.jl
index 98f7a67..7d3c314 100644
--- a/src/RegisterOptimize.jl
+++ b/src/RegisterOptimize.jl
@@ -35,18 +35,22 @@ export
     initial_deformation
 
 """
-This module provides convenience functions for minimizing the mismatch
-between images. It supports both rigid registration and deformable
-registration.
-
-The main functions are:
-
-- `optimize_rigid`: iteratively improve a rigid transformation, given raw images
-- `rotation_gridsearch`: brute-force search a grid of possible rotations and shifts to align raw images
-- `qd_rigid`: find a rotation and shift to align raw images using the QuadDIRECT algorithm
-- `initial_deformation`: provide an initial guess based on mismatch quadratic fits
-- `RegisterOptimize.optimize!`: iteratively improve a deformation, given mismatch data
-- `fixed_λ` and `auto_λ`: "complete" optimizers that generate initial guesses and then find the minimum
+Convenience functions for minimizing the mismatch between images,
+supporting both rigid and deformable registration.
+
+**Exported functions:**
+
+- `initial_deformation`: initial deformation from quadratic mismatch fits
+- `fixed_λ`: full deformable registration at a fixed regularization strength
+- `auto_λ`: full deformable registration with automatic regularization selection
+- `auto_λt`: estimate optimal temporal regularization strength for image sequences
+- `fit_sigmoid`: fit data to a logistic function (used internally by `auto_λ`)
+
+**Internal helpers (not exported):**
+
+- `RegisterOptimize.optimize_rigid`: iteratively refine a rigid transformation
+- `RegisterOptimize.rotation_gridsearch`: brute-force grid search over rotations
+- `RegisterOptimize.optimize!`: iteratively improve a deformation given mismatch data
 """
 RegisterOptimize
 
@@ -320,19 +324,46 @@ MOI.eval_objective_gradient(d::RigidOpt, grad_f, x) =
 ### quadratic-fit mismatch data
 ###
 """
-`u0, converged = initial_deformation(ap::AffinePenalty, cs, Qs; λt=nothing)`
-prepares a globally-optimal initial guess for a deformation, given a
-quadratic fit to the aperture-wise mismatch data. `cs` and `Qs` must
-be arrays-of-arrays in the shape of the u0-grid, each entry as
-calculated by `qfit`. The initial guess minimizes the function
+    u0, converged = initial_deformation(ap::AffinePenalty, cs, Qs; λt=nothing)
+
+Compute the optimal initial deformation `u0` with respect to the
+quadratic approximation of the mismatch data. `cs` and `Qs` are
+arrays-of-arrays matching the deformation grid shape, with each entry
+calculated by `qfit`. `u0` minimizes
 
 ```
     ap(ϕ(u0)) + ∑_i (u0[i]-cs[i])' * Qs[i] * (u0[i]-cs[i])
 ```
-where `ϕ(u0)` is the deformation associated with `u0`.
 
-Pass `λt` to include a temporal roughness penalty (requires
-`cs::AbstractArray{SVector{…}}` and `Qs::AbstractArray{SMatrix{…}}`).
+Because this objective is a quadratic approximation of the true mismatch,
+`u0` is globally optimal only with respect to that approximation. Pass it
+as the starting point for `fixed_λ` or `RegisterOptimize.optimize!`.
+
+`converged` is `true` if the iterative solver converged. When `converged`
+is `false` the returned `u0` may be inaccurate; `fixed_λ` and `auto_λ`
+fall back to using `cs` directly in that case.
+
+Pass `λt` to include a temporal roughness penalty for an image sequence
+(requires `cs::AbstractArray{SVector{…}}` and `Qs::AbstractArray{SMatrix{…}}`).
+
+# Examples
+```jldoctest
+julia> using RegisterPenalty, StaticArrays
+
+julia> nodes = (range(1, stop=10, length=3),);
+
+julia> ap = AffinePenalty(nodes, 1.0);
+
+julia> cs = [[2.0], [-1.0], [0.5]]; Qs = [reshape([1.0], 1, 1) for _ in 1:3];
+
+julia> u0, converged = initial_deformation(ap, cs, Qs);
+
+julia> converged
+true
+
+julia> length(u0) == length(cs)
+true
+```
 """
 function initial_deformation(ap::AffinePenalty, cs, Qs; λt = nothing)
     return _initial_deformation(ap, cs, Qs)
@@ -736,16 +767,17 @@ function _copy!(ϕs::Vector{D}, x::Array{T}) where {D <: GridDeformation, T <: N
 end
 
 """
-`ϕ, penalty = fixed_λ(cs, Qs, nodes, affinepenalty, mmis; λt=nothing)` computes an
-optimal deformation `ϕ` (or sequence `ϕs` when `λt` is given) and its
-total `penalty` (data penalty + regularization penalty).  `cs` and `Qs`
-come from `qfit`, `nodes` specifies the deformation grid, `affinepenalty`
-the `AffinePenalty` object for that grid, and `mmis` is the
-array-of-mismatch arrays (already interpolating, see `interpolate_mm!`).
+    ϕ, penalty = fixed_λ(cs, Qs, nodes, affinepenalty, mmis; λt=nothing)
 
-Pass `λt` to include a temporal roughness penalty for an image sequence
-(requires SVector/SMatrix-typed `cs`/`Qs`); the return value is then a
-`Vector{GridDeformation}`.
+Compute the optimal deformation `ϕ` and scalar total `penalty` (data +
+regularization) for a fixed regularization strength embedded in
+`affinepenalty`. `cs` and `Qs` come from `qfit`, `nodes` specifies the
+deformation grid, `affinepenalty` is an `AffinePenalty` for that grid,
+and `mmis` is the array of mismatch arrays produced by `interpolate_mm!`.
+
+`ϕ` is a `GridDeformation`. Pass `λt` to add a temporal roughness
+penalty for an image sequence (requires `SVector`/`SMatrix`-typed
+`cs`/`Qs`), in which case `ϕ` is a `Vector{GridDeformation}`.
 
 See also: `auto_λ`.
 """
@@ -804,39 +836,34 @@ end
 ### Set λ automatically
 ###
 """
-`ϕ, penalty, λ, datapenalty, quality = auto_λ(fixed, moving, gridsize, maxshift, (λmin, λmax))` automatically chooses "the best"
-value of `λ` to serve in the spatial regularization penalty. It tests a
-sequence of `λ` values, starting with `λmin` and each successive value
-two-fold larger than the previous; for each such `λ`, it optimizes the
-registration and then evaluates just the "data" portion of the
-penalty.  The "best" value is selected by a sigmoidal fit of the
-impact of `λ` on the data penalty, choosing a value that lies at the
-initial upslope of the sigmoid (indicating that the penalty is large
-enough to begin limiting the form of the deformation, but not yet to
-substantially decrease the quality of the registration).
-
-`ϕ, penalty, λ, datapenalty, quality = auto_λ(cs, Qs, nodes, mmis,
-(λmin, λmax))` is used if you've already computed mismatch data. `cs`
-and `Qs` come from `qfit`, `nodes` specifies the deformation grid, and
-`mmis` is the array-of-mismatch arrays (already interpolating, see
-`interpolate_mm!`).
-
-As a first pass, try setting `λmin=1e-6` and `λmax=100`. You can plot
-the returned `datapenalty` and check that it is approximately
-sigmoidal; if not, you will need to alter the range you supply.
-
-Upon return, `ϕ` is the chosen deformation, `penalty` its total
-penalty (data penalty+regularization penalty), `λ` is the chosen value
-of `λ`, `datapenalty` is a vector containing the data penalty for each
-tested `λ` value, and `quality` an estimate (possibly broken) of the
-fidelity of the sigmoidal fit.
-
-If you have data for an image sequence, pass `stackidx=k` to analyze
-only the `k`-th slice of `cs`, `Qs`, and `mmis` along their last
-dimension.
-
-See also: `fixed_λ`. Because `auto_λ` performs the optimization
-repeatedly for many different `λ`s, it is slower than `fixed_λ`.
+    ϕ, penalty, λ, λ_all, datapenalty, quality = auto_λ(fixed, moving, gridsize, maxshift, (λmin, λmax))
+    ϕ, penalty, λ, λ_all, datapenalty, quality = auto_λ(cs, Qs, nodes, mmis, (λmin, λmax))
+
+Automatically choose the spatial regularization strength `λ` for
+deformable image registration. Tests a geometric sequence of `λ` values
+from `λmin` to `λmax` (each step two-fold larger), optimizes the
+registration at each, and selects the best value via a sigmoidal fit of
+the data penalty vs. `λ` curve, targeting the onset of regularization.
+
+The second form accepts pre-computed mismatch data: `cs` and `Qs` from
+`qfit`, `nodes` specifying the deformation grid, and `mmis` the
+interpolated mismatch arrays (see `interpolate_mm!`).
+
+As a first pass, try `λmin=1e-6` and `λmax=100`. Plot the returned
+`datapenalty` vs `λ_all` to verify an approximately sigmoidal shape;
+if not, widen or shift the range. Pass `stackidx=k` to restrict
+analysis to the `k`-th frame of an image sequence.
+
+Returns:
+- `ϕ`: the `GridDeformation` at the chosen `λ`
+- `penalty`: total (data + regularization) penalty at `ϕ`
+- `λ`: the selected regularization strength
+- `λ_all`: `Vector` of all tested `λ` values
+- `datapenalty`: `Vector` of data-only penalties for each tested `λ`
+- `quality`: scalar estimate of the sigmoidal fit quality (lower is better)
+
+See also: `fixed_λ`. Because `auto_λ` optimizes for many `λ` values, it
+is slower than `fixed_λ`.
 """
 function auto_λ(fixed::AbstractArray{R}, moving::AbstractArray{S}, gridsize::NTuple{N}, maxshift::NTuple{N}, λrange; thresh = (0.5)^ndims(fixed) * length(fixed) / prod(gridsize), kwargs...) where {R <: Real, S <: Real, N}
     T = Float64
@@ -957,20 +984,23 @@ end
 
 # Because of the long run times, here we only use the quadratic approximation
 """
-`λts, datapenalty = auto_λt(Es, cs, Qs, ap, (λtmin, λtmax))` estimates
-the whole-experiment mismatch penalty as a function of `λt`, choosing
-values starting at `λtmin` and increasing two-fold until `λtmax`.
-`Es`, `cs`, and `Qs` come from the quadratic fix of the mismatch, and
-`ap` is the (spatial) affine-residual penalty.  As a first guess, try
-`λtmin=1e-6` and `λtmax=1`.  (Larger values of `λt` are noticeably
-slower to optimize.)
-
-By plotting `datapenalty` vs `λts` (with a log-scale on the x-axis),
-you can find the "kink" at which the value of `λt` begins to constrain
-the optimization.  Good choices for `λt` tend to be near this kink.
-Since only an approximation of the mismatch is used, the value of the
-estimated data penalty will not be terribly accurate, but the hope is
-that its dependence on `λt` will be approximately correct.
+    λts, datapenalty = auto_λt(Es, cs, Qs, ap, (λtmin, λtmax))
+
+Estimate the whole-experiment mismatch penalty as a function of the
+temporal regularization strength `λt`, sampling values from `λtmin` to
+`λtmax` in two-fold steps. `Es`, `cs`, and `Qs` come from the quadratic
+fit of the mismatch (e.g., via `mms2fit`), and `ap` is the (spatial)
+`AffinePenalty`. As a first guess, try `λtmin=1e-6` and `λtmax=1`.
+(Larger values of `λt` are noticeably slower to optimize.)
+
+Returns `λts`, a `Vector` of the tested `λt` values, and `datapenalty`,
+a `Vector` of the estimated data penalty at each value. By plotting
+`datapenalty` vs `λts` with a log-scale on the x-axis, you can find
+the "kink" at which `λt` begins to constrain the optimization. Good
+choices for `λt` tend to be near this kink. Since only a quadratic
+approximation of the mismatch is used, the absolute values of
+`datapenalty` will not be accurate, but their dependence on `λt` should
+be approximately correct.
 """
 function auto_λt(Es, cs, Qs, ap, λtrange)
     ngrid = prod(size(Es)[1:(end - 1)])
@@ -1046,14 +1076,35 @@ end
 # Used in automatically setting λ
 
 """
-`fit_sigmoid(data, [bottom, top, center, width])` fits the y-values in `data` to a logistic function
+    fit_sigmoid(data, bottom, top, center, width)
+    fit_sigmoid(data)
+
+Fit the y-values in `data` to a logistic function
+
 ```
    y = bottom + (top-bottom)./(1 + exp(-(data-center)/width))
 ```
-This is "non-extrapolating": the parameter values are constrained to
-be within the range of the supplied data (i.e., `bottom` and `top`
-between the min and max values of `data`, `center` within `[1,
-length(data)]`, and `0.1 <= width <= length(data)`.)
+
+The five-argument form uses the supplied `bottom`, `top`, `center`, and
+`width` as the initial guess for the optimizer. The one-argument form
+derives the initial guess automatically from the data.
+
+All four parameters are constrained to lie within
+`[minimum(data), maximum(data)]`. Throws an error if `length(data) < 4`.
+
+Returns the 5-tuple `(bottom, top, center, width, objective_value)`,
+where `objective_value` is the residual sum of squares at the optimum
+(smaller values indicate a better fit).
+
+# Examples
+```jldoctest
+julia> n = 20; data = 1.0 .+ 9.0 ./ (1 .+ exp.(-(collect(1:n) .- 5.0) ./ 2.0));
+
+julia> b, t, c, w, obj = fit_sigmoid(data);
+
+julia> isfinite(obj) && obj >= 0
+true
+```
 """
 function fit_sigmoid(data, bottom, top, center, width)
     length(data) >= 4 || error("Too few data points for sigmoidal fit")
diff --git a/test/register_optimize.jl b/test/register_optimize.jl
index 39c003b..b606e2d 100644
--- a/test/register_optimize.jl
+++ b/test/register_optimize.jl
@@ -5,6 +5,7 @@ using RegisterCore, RegisterPenalty, RegisterDeformation, RegisterMismatch, Regi
 using Images, CoordinateTransformations, Rotations, RegisterOptimize, LinearAlgebra
 using RegisterUtilities
 using Aqua
+using Documenter
 using ExplicitImports
 
 @testset "Aqua" begin
@@ -418,3 +419,12 @@ end
     # recovered rotation is close to rotation2(-angle)
     @test result_tfm.linear[1, 2] ≈ sin(angle) atol = 0.05
 end
+
+@testset "Doctests" begin
+    DocMeta.setdocmeta!(
+        RegisterOptimize, :DocTestSetup,
+        :(using RegisterOptimize, RegisterPenalty, StaticArrays);
+        recursive = true,
+    )
+    doctest(RegisterOptimize; manual = false)
+end

From c13bcd381f3b3dbb5704a9f667fb33c60ad96d8a Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 10:16:53 -0500
Subject: [PATCH 06/11] Update freshen-package status: coverage and docstrings
 done

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 .claude/freshen-package-status | 11 +++++++++++
 1 file changed, 11 insertions(+)
 create mode 100644 .claude/freshen-package-status

diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status
new file mode 100644
index 0000000..b2bf205
--- /dev/null
+++ b/.claude/freshen-package-status
@@ -0,0 +1,11 @@
+DONE: design review
+DONE: API review
+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

From 7f00b0c264596db84400337fda4e213eef224e74 Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 10:19:27 -0500
Subject: [PATCH 07/11] Add README with installation, concepts, and usage
 examples
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Covers: HolyLabRegistry installation, ecosystem table, concepts
(mismatch data, regularization, rigid vs deformable), and usage
examples for fixed_λ, auto_λ, optimize_rigid, and temporal sequences.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 .claude/freshen-package-status |   2 +-
 README.md                      | 113 +++++++++++++++++++++++++++++++++
 2 files changed, 114 insertions(+), 1 deletion(-)
 create mode 100644 README.md

diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status
index b2bf205..9e5e151 100644
--- a/.claude/freshen-package-status
+++ b/.claude/freshen-package-status
@@ -8,4 +8,4 @@ DONE: add ExplicitImports.jl
 DONE: limit struct mutability
 DONE: improve test coverage
 DONE: add and improve docstrings
-TODO: add or improve documentation
+DONE: add or improve documentation
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..f2e9abc
--- /dev/null
+++ b/README.md
@@ -0,0 +1,113 @@
+# RegisterOptimize.jl
+
+[![CI](https://github.com/HolyLab/RegisterOptimize.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/RegisterOptimize.jl/actions/workflows/CI.yml)
+[![codecov](https://codecov.io/gh/HolyLab/RegisterOptimize.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HolyLab/RegisterOptimize.jl)
+
+RegisterOptimize provides optimization routines for image registration: given a pair of images, it finds the rigid or deformable transformation that minimizes the mismatch between them. It is part of the [HolyLab image registration ecosystem](#ecosystem).
+
+## Installation
+
+This package is registered in the [HolyLabRegistry](https://github.com/HolyLab/HolyLabRegistry). Add that registry once, then install normally:
+
+```julia
+using Pkg
+pkg"registry add https://github.com/HolyLab/HolyLabRegistry.git"
+Pkg.add("RegisterOptimize")
+```
+
+## Ecosystem
+
+RegisterOptimize sits at the optimization end of a multi-package pipeline:
+
+| Package | Role |
+|---------|------|
+| [RegisterMismatch.jl](https://github.com/HolyLab/RegisterMismatch.jl) | Compute block-wise mismatch between image pairs |
+| [RegisterFit.jl](https://github.com/HolyLab/RegisterFit.jl) | Fit quadratic models to mismatch data |
+| [RegisterPenalty.jl](https://github.com/HolyLab/RegisterPenalty.jl) | Regularization penalties for deformations |
+| [RegisterDeformation.jl](https://github.com/HolyLab/RegisterDeformation.jl) | Deformation types and coordinate transformations |
+| **RegisterOptimize.jl** | **Optimize the registration** |
+
+## Concepts
+
+### Mismatch data
+
+The main optimization routines work on *mismatch data* rather than raw images. For each block of the image grid, mismatch data records how much the two images differ as a function of local shift. This representation is precomputed by RegisterMismatch and interpolated before being passed to RegisterOptimize, enabling fast repeated evaluations during optimization.
+
+### Regularization
+
+Deformable registration is ill-posed: without constraints, the optimizer can find deformations that align images locally but are physically unrealistic. A regularization penalty (controlled by strength `λ`) penalizes deformations that deviate from an affine map. Larger `λ` gives smoother but less locally accurate registration; smaller `λ` gives more flexibility but risks overfitting.
+
+`fixed_λ` optimizes at a single user-supplied `λ`. `auto_λ` searches over a range and selects the best value automatically by fitting a sigmoid to the data-penalty curve.
+
+### Rigid vs. deformable
+
+- **Rigid** (`RegisterOptimize.optimize_rigid`): finds the rotation + translation that best aligns two raw images. Useful as an initial pass when images differ primarily by global motion.
+- **Deformable** (`fixed_λ`, `auto_λ`): finds a smooth spatially-varying warp field that aligns images locally.
+
+## Usage
+
+### Deformable registration at a fixed λ
+
+```julia
+using RegisterMismatch, RegisterFit, RegisterPenalty, RegisterOptimize
+
+# 1. Compute block-wise mismatch between images
+mms = mismatch_apertures(fixed, moving, aperture_centers, aperture_width, maxshift;
+                         normalization=:pixels)
+
+# 2. Fit quadratic models; get interpolated mismatch arrays
+thresh = 0.5^ndims(fixed) * length(fixed) / prod(gridsize)
+cs, Qs, mmis = mms2fit(mms, thresh)
+
+# 3. Set up a deformation grid and regularization penalty
+nodes = map(d -> range(1, stop=size(fixed, d), length=gridsize[d]), 1:ndims(fixed))
+ap = AffinePenalty(nodes, λ)   # λ controls regularization strength
+
+# 4. Optimize: returns the deformation and its total penalty
+ϕ, penalty = fixed_λ(cs, Qs, nodes, ap, mmis)
+```
+
+`ϕ` is a `GridDeformation` from RegisterDeformation.jl. Apply it with `transform(moving, ϕ)`.
+
+### Automatic λ selection
+
+When you don't know a suitable `λ` in advance, let `auto_λ` find it:
+
+```julia
+ϕ, penalty, λ_best, λ_all, datapenalty, quality =
+    auto_λ(fixed, moving, gridsize, maxshift, (λmin, λmax))
+```
+
+`auto_λ` tests a logarithmic sweep of `λ` values from `λmin` to `λmax`, fits a sigmoid to the data-penalty curve, and returns the deformation at the automatically chosen `λ_best`. Plot `datapenalty` vs `λ_all` (log x-axis) to verify the sigmoidal shape; if it is not sigmoidal, widen the range. A good first guess is `(λmin, λmax) = (1e-6, 100)`.
+
+If mismatch data are already computed, pass them directly to avoid recomputing:
+
+```julia
+ϕ, penalty, λ_best, λ_all, datapenalty, quality =
+    auto_λ(cs, Qs, nodes, mmis, (λmin, λmax))
+```
+
+### Rigid registration
+
+```julia
+using RegisterOptimize, CoordinateTransformations, LinearAlgebra
+
+tform0 = AffineMap(Matrix(1.0I, 2, 2), zeros(2))   # identity initial guess
+tform, fval = RegisterOptimize.optimize_rigid(fixed, moving, tform0, maxshift)
+```
+
+The returned `tform` is an `AffineMap` (rotation + translation) and `fval` is the normalized sum-of-squared-differences at the solution.
+
+### Temporal image sequences
+
+For a time series of images, add a temporal roughness penalty `λt` to penalize frame-to-frame variation in the deformation:
+
+```julia
+# Fixed spatial λ and temporal λt
+ϕs, penalty = fixed_λ(cs, Qs, nodes, ap, mmis; λt=0.1)
+
+# Or find a good λt automatically from the quadratic approximation
+λts, datapenalty = auto_λt(Es, cs, Qs, ap, (1e-6, 1.0))
+```
+
+`ϕs` is a `Vector{GridDeformation}`, one per frame. Plot `datapenalty` vs `λts` (log x-axis) to find the "kink" where `λt` begins to constrain the optimization.

From abe4de596c5fda365a99c0f11dc6bca91dd1931c Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 10:20:05 -0500
Subject: [PATCH 08/11] Remove freshen-package-status (all steps complete)

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 .claude/freshen-package-status | 11 -----------
 1 file changed, 11 deletions(-)
 delete mode 100644 .claude/freshen-package-status

diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status
deleted file mode 100644
index 9e5e151..0000000
--- a/.claude/freshen-package-status
+++ /dev/null
@@ -1,11 +0,0 @@
-DONE: design review
-DONE: API review
-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
-DONE: add or improve documentation

From 49c69a3c05e6b531a8d2d0ab6cad4ae70513238d Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 10:23:10 -0500
Subject: [PATCH 09/11] Remove CHANGELOG.md

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 CHANGELOG.md | 63 ----------------------------------------------------
 1 file changed, 63 deletions(-)
 delete mode 100644 CHANGELOG.md

diff --git a/CHANGELOG.md b/CHANGELOG.md
deleted file mode 100644
index fa039af..0000000
--- a/CHANGELOG.md
+++ /dev/null
@@ -1,63 +0,0 @@
-# Changelog
-
-## v1.0.0
-
-**Breaking changes**
-
-- `initial_deformation(ap, cs, Qs, ϕ_old, maxshift)` removed. The method body was
-  `error("This is broken, don't use it")` and is gone entirely.
-
-- `optimize(tform::AffineMap, mmis, nodes)` removed. Used a deprecated Optim API
-  (`interior`/`DifferentiableFunction`/`At_mul_Bt!`) that no longer exists.
-
-- `auto_λ`: the `stackidx::Integer` positional-first-argument overload is removed.
-  Pass `stackidx` as a keyword instead:
-  ```julia
-  # before
-  auto_λ(stackidx, cs, Qs, nodes, mmis, λrange)
-  # after
-  auto_λ(cs, Qs, nodes, mmis, λrange; stackidx)
-  ```
-
-- `optimize_rigid`: `SD` and `maxrot` are now keyword arguments (no defaults
-  change); `tol` is renamed `atol`; additional keyword arguments are forwarded
-  to Ipopt:
-  ```julia
-  # before
-  optimize_rigid(fixed, moving, tform0, maxshift, SD, maxrot; tol=1e-4)
-  # after
-  optimize_rigid(fixed, moving, tform0, maxshift; SD, maxrot, atol=1e-4)
-  ```
-
-- `rotation_gridsearch`: `SD` is now a keyword argument:
-  ```julia
-  # before
-  rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD)
-  # after
-  rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD)
-  ```
-
-- `initial_deformation`, `optimize!` (time-series), and `fixed_λ`: the temporal
-  penalty coefficient `λt` is now a keyword argument (default `nothing`):
-  ```julia
-  # before
-  initial_deformation(ap, λt, cs, Qs)
-  optimize!(ϕs, ϕs_old, dp, λt, mmis)
-  fixed_λ(cs, Qs, nodes, ap, λt, mmis)
-  # after
-  initial_deformation(ap, cs, Qs; λt)
-  optimize!(ϕs, ϕs_old, dp, mmis; λt)
-  fixed_λ(cs, Qs, nodes, ap, mmis; λt)
-  ```
-
-- `auto_λt` flat-array adapter overload removed. Pass SVector/SMatrix-typed
-  `cs`/`Qs` directly to the main `auto_λt(Es, cs, Qs, ap, λtrange)` overload.
-
-**Non-breaking improvements**
-
-- `auto_λ` flat-array overloads now accept `AbstractArray{<:Number}` and
-  `AbstractArray{Float64}` instead of requiring `Array{Tf}` / `Array{Float64}`,
-  allowing views and other array wrappers.
-
-- `auto_λ` flat-array overload: the element types of `cs`, `Qs`, and `mmis` are
-  now independent (previously all three had to share the same `Tf`).

From ce45ee65d884bd7c9bda94d04eba93ffa87e28a2 Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 10:31:12 -0500
Subject: [PATCH 10/11] Remove stale import of Linear from Interpolations

The _optimize! dispatch was widened to BSpline{<:Interpolations.Linear},
which uses the qualified name rather than the imported symbol.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 src/RegisterOptimize.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/RegisterOptimize.jl b/src/RegisterOptimize.jl
index 7d3c314..e33bff8 100644
--- a/src/RegisterOptimize.jl
+++ b/src/RegisterOptimize.jl
@@ -7,7 +7,7 @@ using CenterIndexedArrays: CenterIndexedArrays, CenterIndexedArray
 using CoordinateTransformations: CoordinateTransformations, AffineMap
 using ForwardDiff: ForwardDiff
 using Interpolations: Interpolations, AbstractExtrapolation, AbstractInterpolation,
-                      BSpline, InPlace, Linear, OnCell, Quadratic
+                      BSpline, InPlace, OnCell, Quadratic
 using Ipopt: Ipopt
 using IterativeSolvers: IterativeSolvers, cg
 using LinearAlgebra: LinearAlgebra, I, dot, mul!, tr

From a7be6c5eab133bcdc01261e2b36c44190972eca0 Mon Sep 17 00:00:00 2001
From: Dae Woo Kim <kim503@wustl.edu>
Date: Wed, 20 May 2026 10:47:50 -0500
Subject: [PATCH 11/11] =?UTF-8?q?Remove=20const=20fields:=20revert=20mutab?=
 =?UTF-8?q?le=20struct=20=E2=86=92=20struct?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 src/RegisterOptimize.jl | 56 ++++++++++++++++++++---------------------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/src/RegisterOptimize.jl b/src/RegisterOptimize.jl
index e33bff8..e4a656a 100644
--- a/src/RegisterOptimize.jl
+++ b/src/RegisterOptimize.jl
@@ -273,12 +273,12 @@ to_float(::Type{T}, A, B) where {T} = convert(Array{Float32}, A), convert(Array{
 ###
 ### Rigid registration from raw images, MathProg interface
 ###
-mutable struct RigidValue{N, A <: AbstractArray, I <: AbstractExtrapolation, SDT} <: MOI.AbstractNLPEvaluator
-    const fixed::A
-    const wfixed::A
-    const moving::I
-    const SD::SDT
-    const thresh
+struct RigidValue{N, A <: AbstractArray, I <: AbstractExtrapolation, SDT} <: MOI.AbstractNLPEvaluator
+    fixed::A
+    wfixed::A
+    moving::I
+    SD::SDT
+    thresh
 end
 
 function RigidValue(fixed::AbstractArray, moving::AbstractArray{T}, SD, thresh) where {T <: Real}
@@ -304,9 +304,9 @@ function (d::RigidValue)(x)
     return sum(abs2, f - m) / den
 end
 
-mutable struct RigidOpt{RV <: RigidValue, G} <: GradOnlyBoundsOnly
-    const rv::RV
-    const g::G
+struct RigidOpt{RV <: RigidValue, G} <: GradOnlyBoundsOnly
+    rv::RV
+    g::G
 end
 
 function RigidOpt(fixed, moving, SD, thresh)
@@ -456,10 +456,10 @@ function find_opt(P, b)
 end
 
 # A type for computing multiplication by the linear operator
-mutable struct AffineQHessian{AP <: AffinePenalty, M <: StaticMatrix, N, Φ}
-    const ap::AP
-    const Qs::Array{M, N}
-    const ϕ_old::Φ
+struct AffineQHessian{AP <: AffinePenalty, M <: StaticMatrix, N, Φ}
+    ap::AP
+    Qs::Array{M, N}
+    ϕ_old::Φ
 end
 
 function AffineQHessian(ap::AffinePenalty{T}, Qs::AbstractArray{TQ, N}, ϕ_old) where {T, TQ, N}
@@ -679,11 +679,11 @@ function u_as_vec(ϕs::Vector{D}, ::Type{T} = eltype(D)) where {D <: GridDeforma
 end
 
 
-mutable struct DeformOpt{D, Dold, DP, M} <: GradOnlyBoundsOnly
-    const ϕ::D
-    const ϕ_old::Dold
-    const dp::DP
-    const mmis::M
+struct DeformOpt{D, Dold, DP, M} <: GradOnlyBoundsOnly
+    ϕ::D
+    ϕ_old::Dold
+    dp::DP
+    mmis::M
 end
 # (d::DeformOpt)(x) = MOI.eval_objective(d, x)
 
@@ -735,12 +735,12 @@ function optimize!(ϕs::Vector{<:GridDeformation}, ϕs_old, dp::AffinePenalty{T,
     return optimize!(ϕs, ϕs_old, dp, mmisc; λt, kwargs...)
 end
 
-mutable struct DeformTseriesOpt{D, Dsold, DP, T, M} <: GradOnlyBoundsOnly
-    const ϕs::Vector{D}
-    const ϕs_old::Dsold
-    const dp::DP
-    const λt::T
-    const mmis::M
+struct DeformTseriesOpt{D, Dsold, DP, T, M} <: GradOnlyBoundsOnly
+    ϕs::Vector{D}
+    ϕs_old::Dsold
+    dp::DP
+    λt::T
+    mmis::M
 end
 
 # Using MOI is a legacy of the old Ipopt interface, but
@@ -1157,10 +1157,10 @@ function fit_sigmoid(data)
 end
 
 
-mutable struct SigmoidOpt{G, H} <: BoundsOnly
-    const data::Vector{Float64}
-    const g::G
-    const h::H
+struct SigmoidOpt{G, H} <: BoundsOnly
+    data::Vector{Float64}
+    g::G
+    h::H
 end
 
 SigmoidOpt(data::Vector{Float64}) = SigmoidOpt(data, y -> ForwardDiff.gradient(x -> sigpenalty(x, data), y), y -> ForwardDiff.hessian(x -> sigpenalty(x, data), y))