diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 1b512f8..ff62964 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -13,7 +13,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.10'
+          - 'min'
           - '1'
 #          - 'nightly'
         os:
@@ -23,7 +23,7 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
@@ -47,6 +47,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v5
         with:
           file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
\ No newline at end of file
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
index 9880dd0..aabc045 100644
--- a/.github/workflows/TagBot.yml
+++ b/.github/workflows/TagBot.yml
@@ -28,8 +28,5 @@ jobs:
     steps:
       - uses: JuliaRegistries/TagBot@v1
         with:
-          token: ${{ secrets.GITHUB_TOKEN }}
-          # Edit the following line to reflect the actual name of the GitHub Secret containing your private key
-          ssh: ${{ secrets.DOCUMENTER_KEY }}
-          # ssh: ${{ secrets.NAME_OF_MY_SSH_PRIVATE_KEY_SECRET }}
+          token: ${{ secrets.TAGBOT_TOKEN }}
           registry: HolyLab/HolyLabRegistry
\ No newline at end of file
diff --git a/API_REVIEW_PLAN.md b/API_REVIEW_PLAN.md
new file mode 100644
index 0000000..f844057
--- /dev/null
+++ b/API_REVIEW_PLAN.md
@@ -0,0 +1,190 @@
+# API Review Plan
+<!-- Auto-generated by /review-api. Edit freely, but preserve chunk IDs and status values. -->
+
+## Metadata
+- **Kind**: `api`
+- **Package**: RegisterOptimize
+- **Source review date**: 2026-05-18
+- **Current version**: 0.3.6
+
+## Stated values
+Aggressive modernization. Package is pre-1.0 (v0.3.6) so breaking changes are cheap. All Tier 1 breaking changes, all Tier 2 non-breaking improvements, and all Tier 3 internal consistency fixes are approved. No deprecation shims — clean break throughout.
+
+## Release strategy
+- **Pre-breaking-release**: `no`
+- **Inter-cluster releases**: `no` — all breaking changes batch into a single v1.0.0
+
+## Baseline
+- Tests pass on the starting commit: `yes` (200/200 passing after prep work; see CHUNK-001 notes)
+- `Test.detect_ambiguities(RegisterOptimize)` count: `0`
+- Working tree clean: `yes` (at commit `7f64c5a`)
+- Current `Project.toml` version: `0.3.6`
+
+## Decisions
+<!-- Answers to `decide` chunks land here, with the chunk ID. -->
+
+## Chunks
+
+### CHUNK-001: preflight
+- **Kind**: `preflight`
+- **Originating finding**: n/a
+- **Cluster**: none
+- **Breaking**: no
+- **Description**: Establish baseline — confirm tests pass, record `Test.detect_ambiguities` count, confirm working tree is clean.
+- **Depends on**: none
+- **Verification**: full test suite green; `Test.detect_ambiguities(RegisterOptimize)` count recorded
+- **Status**: `complete`
+- **Notes**: Working tree was dirty at session start. Three prep commits landed before baseline could be recorded:
+  1. `fea6bbd` — CI/TagBot modernization (julia-actions v2, codecov v5, TagBot to HolyLabRegistry).
+  2. `ed25eee` — `Project.toml` compat bumps for `CachedInterpolations`, `CenterIndexedArrays`, and the `Register*` family to `1`.
+  3. `7f64c5a` — Source/test fixes required by the v1 deps: `TimeHessian` `*`/`mul!` now call `penalty!(y, ϕs, P.λt)` (arg order swapped per RegisterPenalty v1); the temporal-penalty test now calls `quadratic(denom, cs[:,t], Qs[:,:,t])` (denom-first per RegisterUtilities v1). None of these are part of the API review scope but were necessary to get a green baseline.
+
+  Baseline tests: 200 pass / 0 fail / 0 error (Julia 1.12.6, 4 test groups: Register Optimize 148, Minimization to mismatch data 36, Optimization with a temporal penalty 12, Optimization with linear interpolation of mismatch data 4). Ambiguity count: 0. Test run emits a deprecation warning from Optim about `x_tol` — pre-existing, not addressed here.
+
+---
+
+### CHUNK-002: remove-dead-code
+- **Kind**: `implement`
+- **Originating finding**: Tier 1 / T1-D (`initial_deformation` Method 3, line 462) and T1-E (`optimize` function, line 1115)
+- **Cluster**: dead-code
+- **Breaking**: yes
+- **Description**: Remove two exported/semi-public methods that are already broken at runtime:
+  1. `initial_deformation(ap, cs, Qs, ϕ_old, maxshift)` (Method 3) — body immediately calls `error("This is broken, don't use it")`.
+  2. `optimize(tform::AffineMap, mmis, nodes)` — docstring says "not updated yet, probably broken"; uses old Optim API.
+  Remove both method definitions and any associated tests that exist solely to exercise the broken path. Update docstrings for `initial_deformation` to no longer mention Method 3.
+- **Depends on**: CHUNK-001
+- **Verification**: tests pass; `initial_deformation` still works for remaining methods; no references to the removed signatures remain
+- **Status**: `not-started`
+- **Notes**:
+
+---
+
+### CHUNK-003: stackidx-to-keyword
+- **Kind**: `implement`
+- **Originating finding**: Tier 1 / T1-A (`auto_λ` Method 3, line 906)
+- **Cluster**: none
+- **Breaking**: yes
+- **Description**: `auto_λ(stackidx::Integer, cs, Qs, nodes, mmis, λrange; kwargs...)` puts a configuration integer (`stackidx`) before data. Modern Julia convention is data-first.
+  Proposed: remove this overload as a separate method and instead add `stackidx` as a keyword to the relevant `auto_λ` overload(s):
+  `auto_λ(cs, Qs, nodes, mmis, λrange; stackidx=nothing, kwargs...)`.
+  When `stackidx` is not `nothing`, slice `cs`, `Qs`, `mmis` along the stack dimension before proceeding. Migrate or remove any tests for the old positional form.
+- **Depends on**: CHUNK-001
+- **Verification**: existing `auto_λ` calls without `stackidx` unaffected; new keyword form works; tests pass
+- **Status**: `not-started`
+- **Notes**:
+
+---
+
+### CHUNK-004: optimize-rigid-modernize
+- **Kind**: `implement`
+- **Originating finding**: Tier 1 / T1-B (line 81), Tier 2 / T2-A, T2-B
+- **Cluster**: rigid-search
+- **Breaking**: yes
+- **Description**: Three related changes to `optimize_rigid` (line 81):
+  1. **T1-B — positional-to-keyword**: `SD` and `maxrot` are currently positional arguments with defaults. Move both to keyword arguments:
+     `optimize_rigid(fixed, moving, A::AffineMap, maxshift; SD=I_matrix, maxrot=π, atol=1e-4, print_level=0, max_iter=3000, thresh=0, kwargs...)`
+     where `I_matrix` is constructed from `ndims(fixed)` at call time.
+  2. **T2-A — `tol` → `atol`**: Rename `tol` to `atol` (consistent with `Base.isapprox`, `IterativeSolvers`, etc.). No deprecation shim — clean rename.
+  3. **T2-B — `kwargs...` passthrough**: Replace the fixed keyword list (minus `atol`, `print_level`, `max_iter`, `thresh` which stay as explicit named keywords for discoverability) with a forwarded `kwargs...` so callers can pass any Ipopt option.
+  Update callers of `optimize_rigid` within the package (e.g., `rotation_gridsearch` calls it).
+- **Depends on**: CHUNK-001
+- **Verification**: `optimize_rigid` works with all-keyword call; `rotation_gridsearch` (which calls it) still works; tests pass
+- **Status**: `not-started`
+- **Notes**: `I_matrix` default for `SD` should be lazy — compute inside the function body as `Matrix{Float64}(I, ndims(fixed), ndims(fixed))`.
+
+---
+
+### CHUNK-005: rotation-gridsearch-SD-keyword
+- **Kind**: `implement`
+- **Originating finding**: Tier 1 / T1-C (`rotation_gridsearch`, line 191)
+- **Cluster**: rigid-search
+- **Breaking**: yes
+- **Description**: `rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD=...)` — `SD` is a positional argument with a default. Move `SD` to a keyword argument:
+  `rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD=I_matrix)`
+  where `I_matrix` is constructed lazily from `ndims(fixed)` inside the function body (as for CHUNK-004).
+- **Depends on**: CHUNK-004
+- **Verification**: existing calls without `SD` unaffected; explicit-SD calls updated to keyword form; tests pass
+- **Status**: `not-started`
+- **Notes**: Depends on CHUNK-004 because `rotation_gridsearch` calls `optimize_rigid`, and the `SD` default construction pattern should match.
+
+---
+
+### CHUNK-006: lambda-t-unify
+- **Kind**: `implement`
+- **Originating finding**: Tier 3 / T3-A — `λt` placement inconsistency across `initial_deformation` Method 4 (line 1052), `optimize!` Methods 3–4 (lines 730/743), `fixed_λ` Method 2 (line 817)
+- **Cluster**: lambda-t
+- **Breaking**: yes
+- **Description**: `λt` (the temporal penalty coefficient) appears at different positional slots in different functions:
+  - `initial_deformation(ap, λt, cs, Qs)` — position 2
+  - `optimize!(ϕs, ϕs_old, dp, λt, mmis)` — position 4
+  - `fixed_λ(cs, Qs, nodes, ap, λt, mmis)` — position 5
+
+  Unify by making `λt` a keyword argument across all three functions with a default of `nothing` (or `0`) meaning "no temporal penalty." This collapses the temporal and non-temporal overloads:
+  - `initial_deformation(ap, cs, Qs; λt=nothing)` — dispatch to `TimeHessian` internally when `λt !== nothing`
+  - `optimize!(ϕs, ϕs_old, dp, mmis; λt=nothing, kwargs...)` — use time-series path when `λt !== nothing`
+  - `fixed_λ(cs, Qs, nodes, ap, mmis; λt=nothing, ϕs_old=identity, mu_init=0.1, kwargs...)` — unified signature
+
+  Remove the now-redundant separate overloads. Update all internal callers. Update tests.
+- **Depends on**: CHUNK-001
+- **Verification**: single-timepoint calls (`λt` omitted) work; time-series calls (`λt=value`) work; ambiguity count unchanged or reduced; tests pass
+- **Status**: `not-started`
+- **Notes**: This is the most complex chunk — `TimeHessian` construction and dispatch currently happens via overload selection; it needs to move inside the function body as a conditional. Implement carefully with a test for each path.
+
+---
+
+### CHUNK-007: auto-lambda-t-interface
+- **Kind**: `decide`
+- **Originating finding**: Tier 3 / T3-B — `auto_λt` Method 1 (line 1019) accepts `λtrange`, Method 2 (line 1043) accepts scalar `λt`
+- **Cluster**: lambda-t
+- **Breaking**: no (decide chunk — no changes made until decision is recorded)
+- **Description**: `auto_λt` has two methods with different last-argument semantics:
+  - Method 1: `auto_λt(Es, cs, Qs, ap, λtrange)` — iterates over a range
+  - Method 2: `auto_λt(Es, cs::Array{Tf}, Qs::Array{Tf}, ap, λt)` — takes a scalar and repacks raw flat arrays
+
+  Method 2 appears to be a raw-format repacking helper (like the `Array{Tf}` overloads elsewhere). If it delegates to Method 1 with `λtrange = λt:λt` (a trivial range), the scalar `λt` here represents a fixed value, not a range to search — which makes calling `auto_λt` with a scalar semantically different from calling it with a range, and the function name `auto_λt` ("auto-select λt") doesn't make sense for a scalar input.
+
+  **Decision needed**: Should Method 2 (a) be renamed to reflect that it's a raw-format adapter, (b) be removed and callers updated to use Method 1 directly, or (c) keep as-is with improved documentation?
+- **Depends on**: CHUNK-006
+- **Verification**: depends on decision
+- **Status**: `not-started`
+- **Notes**: Record the decision in the `## Decisions` section above with this chunk ID, then the implementer will act on it.
+
+---
+
+### CHUNK-008: widen-array-types
+- **Kind**: `implement`
+- **Originating finding**: Tier 2 / T2-C (`auto_λ` Method 5, line 918) and Tier 3 / T3-C (`auto_λ` Method 4, line 913)
+- **Cluster**: none
+- **Breaking**: no
+- **Description**: Two related type-widening changes in legacy-format `auto_λ` overloads:
+  1. **T2-C** (line 918): `cs::Array{Float64}`, `Qs::Array{Float64}`, `mmis::Array{Float64}` — widen to `AbstractArray{Float64}` (or `AbstractArray{<:AbstractFloat}` if the downstream code tolerates it).
+  2. **T3-C** (line 913): `cs::Array{Tf}`, `Qs::Array{Tf}`, `mmis::Array{Tf}` where `Tf<:Number` — all three share the same `Tf`. Relax to independent type parameters or use `promote_type` inside the function body.
+  Neither change breaks existing callers (widening is non-breaking).
+- **Depends on**: CHUNK-001
+- **Verification**: existing tests pass; confirm `Float32` arrays are now accepted where they previously weren't
+- **Status**: `not-started`
+- **Notes**: Be careful about ambiguity with other overloads when widening — run `Test.detect_ambiguities` after.
+
+---
+
+### CHUNK-009: version-bump
+- **Kind**: `version-bump`
+- **Originating finding**: n/a
+- **Cluster**: none
+- **Breaking**: yes
+- **Description**: Bump version from 0.3.6 → 1.0.0 (first stable release, incorporating all breaking changes). Update CHANGELOG (or create one if absent). Tag the release.
+- **Depends on**: CHUNK-002, CHUNK-003, CHUNK-004, CHUNK-005, CHUNK-006, CHUNK-007, CHUNK-008
+- **Verification**: full test suite green; `Project.toml` version updated; no half-finished cluster; git tag created
+- **Status**: `not-started`
+- **Notes**: Only cut after all other chunks are complete and CHUNK-007 decision is resolved.
+
+## Dropped findings
+<!-- Items the user chose not to act on, with one-line reasons. Preserved for institutional memory. -->
+*(none — all findings approved)*
+
+## Session Log
+- 2026-05-18: Initial review and plan generated by /review-api. 9 chunks total (1 preflight, 6 implement, 1 decide, 1 version-bump).
+- 2026-05-18: Implemented CHUNK-001 (preflight). Baseline: 200/200 tests pass, 0 ambiguities at commit `7f64c5a`. Pre-existing uncommitted work (CI modernization, Register* compat bump to v1, and adapting two call sites to the v1 `penalty!` / `quadratic` signatures) was committed first to obtain a green baseline. Next up: CHUNK-002 (remove-dead-code).
+
+## Open Questions
+- CHUNK-007: How should the `auto_λt` scalar vs. range interface inconsistency be resolved? See chunk for options.
diff --git a/API_REVIEW_SESSION.md b/API_REVIEW_SESSION.md
new file mode 100644
index 0000000..516c94b
--- /dev/null
+++ b/API_REVIEW_SESSION.md
@@ -0,0 +1,52 @@
+# Session Handoff — 2026-05-18
+
+## Plan
+API_REVIEW_PLAN.md — RegisterOptimize v0.3.6
+
+## What was just completed
+CHUNK-001: preflight
+Established the baseline. Working tree had unrelated in-progress work; landed
+three prep commits (CI/TagBot modernization, Register* compat bump to v1, and
+two-line adaptations of `TimeHessian` and a test for the v1 `penalty!` and
+`quadratic` signatures) before the test suite went green. Final baseline:
+200/200 tests pass, 0 ambiguities, at commit `7f64c5a`.
+
+## Key decisions / shim choices
+- User chose to commit the dirty working tree's changes rather than stash or
+  revert, then to fix forward against RegisterPenalty/RegisterUtilities v1
+  rather than pin to older deps.
+- The Optim `x_tol` deprecation warning is pre-existing and intentionally not
+  addressed in CHUNK-001.
+
+## State of the codebase
+- Files modified (since session start, across three commits):
+  `.github/workflows/CI.yml`, `.github/workflows/TagBot.yml`, `Project.toml`,
+  `src/RegisterOptimize.jl`, `test/register_optimize.jl`.
+- Test suite: 200/200 pass (Julia 1.12.6).
+- Ambiguity count: 0 (baseline).
+- Staged but uncommitted: no. Untracked: `.claude/`, `API_REVIEW_PLAN.md`,
+  `API_REVIEW_SESSION.md`.
+
+## Cluster status
+- (no cluster — preflight)
+- Remaining clusters: `dead-code` (1 chunk), `rigid-search` (2 chunks),
+  `lambda-t` (2 chunks), plus standalone CHUNK-003, CHUNK-008, CHUNK-009.
+
+## Next chunk
+CHUNK-002: remove-dead-code — remove `initial_deformation` Method 3 (line
+462, `error("This is broken…")`) and the broken `optimize(::AffineMap, …)`
+overload (line ~1115). Breaking change.
+
+## Watch out for
+- Baseline was reached only after fixing-forward to RegisterPenalty v1 and
+  RegisterUtilities v1. Other call sites of `penalty!` (lines 442, 449, 549,
+  556, 600, 612, 708, 714, 762, 767, 965, 969, 987, 991, 1174) all use the
+  longer signatures and were unaffected — but any future change in the v1
+  signatures could break them.
+- CHUNK-006 (lambda-t-unify) touches the same `TimeHessian` machinery whose
+  arg-order was just fixed. Double-check that `penalty!(y, ϕs, λt)` continues
+  to apply when `λt` becomes optional.
+- Julia 1.12 is in use locally; CI matrix is `min` (= 1.10 per Project.toml
+  `julia` compat) and `1` (latest). The local Manifest currently resolves
+  Stdlib versions to 1.11 / 1.12 — this is fine for testing but worth
+  remembering if results ever diverge.
diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..fa039af
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,63 @@
+# 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`).
diff --git a/Project.toml b/Project.toml
index d040f59..51cda05 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 = "0.3.6"
+version = "1.0.0"
 
 [deps]
 CachedInterpolations = "b9709bfb-d23d-5560-8276-8c35c4b76823"
@@ -27,8 +27,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-CachedInterpolations = "0.1"
-CenterIndexedArrays = "0.2"
+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"
@@ -39,10 +39,10 @@ MathOptInterface = "1"
 NLsolve = "3, 4"
 Optim = "0.18, 0.19, 0.20, 1"
 ProgressMeter = "0.7, 0.8, 0.9, 1"
-RegisterCore = "0.2"
-RegisterDeformation = "0.3, 0.4"
-RegisterFit = "0.1, 0.2"
-RegisterPenalty = "0.1, 0.2, 0.3"
+RegisterCore = "1"
+RegisterDeformation = "1"
+RegisterFit = "1"
+RegisterPenalty = "1"
 StaticArrays = "0.10, 0.11, 0.12, 1"
 Statistics = "1"
 julia = "1.10"
diff --git a/src/RegisterOptimize.jl b/src/RegisterOptimize.jl
index 7b50fb9..e1e2c4a 100644
--- a/src/RegisterOptimize.jl
+++ b/src/RegisterOptimize.jl
@@ -65,8 +65,7 @@ MOI.eval_constraint_jacobian(::BoundsOnly, J, x) = nothing
 ### Rigid registration from raw images
 ###
 """
-`tform = optimize_rigid(fixed, moving, tform0, maxshift, [SD = eye];
-[thresh=0, tol=1e-4, print_level=0])` optimizes a rigid transformation
+`tform = optimize_rigid(fixed, moving, tform0, maxshift; kwargs...)` optimizes a rigid transformation
 (rotation + shift) to minimize the mismatch between `fixed` and
 `moving`.
 
@@ -77,10 +76,19 @@ enforces a certain amount of sum-of-squared-intensity overlap between
 the two images; with non-zero `thresh`, it is not permissible to
 "align" the images by shifting one entirely out of the way of the
 other.
+
+Keyword arguments:
+- `SD`: sample-spacing matrix (default: identity of size `ndims(fixed)`)
+- `maxrot`: maximum rotation in radians (default: `π`)
+- `atol`: Ipopt convergence tolerance (default: `1e-4`)
+- `print_level`: Ipopt verbosity (default: `0`)
+- `max_iter`: maximum Ipopt iterations (default: `3000`)
+- `thresh`: minimum intensity-overlap threshold (default: `0`)
+- Any additional keyword arguments are forwarded to Ipopt as options.
 """
-function optimize_rigid(fixed, moving, A::AffineMap, maxshift,
-                SD = Matrix{Float64}(I,size(A.linear,1),size(A.linear,1)),
-                maxrot=pi; thresh=0, tol=1e-4, print_level=0, max_iter=3000)
+function optimize_rigid(fixed, moving, A::AffineMap, maxshift;
+                SD=nothing, maxrot=π, atol=1e-4, print_level=0, max_iter=3000, thresh=0, kwargs...)
+    SD === nothing && (SD = Matrix{Float64}(I, ndims(fixed), ndims(fixed)))
     objective = RigidOpt(to_float(fixed, moving)..., SD, thresh)
     # Convert initial guess into parameter vector
     R = SD*A.linear/SD
@@ -92,10 +100,11 @@ function optimize_rigid(fixed, moving, A::AffineMap, maxshift,
     # Set up and run the solver
     model = Model(optimizer_with_attributes(Ipopt.Optimizer,
                                             "hessian_approximation" => "limited-memory",
-                                             "print_level" => print_level,
-                                             "tol" => tol,
-                                             "max_iter" => max_iter,
-                                             "sb" => "yes"))
+                                            "print_level" => print_level,
+                                            "tol" => atol,
+                                            "max_iter" => max_iter,
+                                            "sb" => "yes",
+                                            (string(k) => v for (k, v) in kwargs)...))
 
     # ub = T[fill(maxrot, length(p0)-length(maxshift)); [maxshift...]]
     # MOI.loadproblem!(model, length(p0), 0, -ub, ub, T[], T[], :Min, objective)
@@ -180,15 +189,16 @@ function grid_rotations(maxradians, rgridsz, SD)
 end
 
 """
-`best_tform, best_mm = rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD =Matrix{Float64}(I,ndims(fixed),ndims(fixed))))`
+`best_tform, best_mm = rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD=I_matrix)`
 Tries a grid of rotations to align `moving` to `fixed`.  Also calculates the best translation (`maxshift` pixels
 or less) to align the images after performing the rotation. Returns an AffineMap that captures both the
 best rotation and shift out of the values searched, along with the mismatch value after applying that transform (`best_mm`).
 
-For more on how the arguments `maxradians`, `rgridsz`, and `SD` influence the search, see the documentation for
-`grid_rotations`.
+`SD` defaults to the identity matrix of size `ndims(fixed)`. For more on how `maxradians`, `rgridsz`,
+and `SD` influence the search, see the documentation for `grid_rotations`.
 """
-function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD = Matrix{Float64}(I,ndims(fixed),ndims(fixed)))
+function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD=nothing)
+    SD === nothing && (SD = Matrix{Float64}(I, ndims(fixed), ndims(fixed)))
     rgridsz = [rgridsz...]
     nd = ndims(moving)
     @assert nd == ndims(fixed)
@@ -284,21 +294,21 @@ MOI.eval_objective_gradient(d::RigidOpt, grad_f, x) =
 ### quadratic-fit mismatch data
 ###
 """
-`u0, converged = initial_deformation(ap::AffinePenalty, cs, Qs;
-[ϕ_old=identity])` 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)`
+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
 
 ```
     ap(ϕ(u0)) + ∑_i (u0[i]-cs[i])' * Qs[i] * (u0[i]-cs[i])
 ```
 where `ϕ(u0)` is the deformation associated with `u0`.
 
-If `ϕ_old` is not the identity, it must be interpolating.
+Pass `λt` to include a temporal roughness penalty (requires
+`cs::AbstractArray{SVector{…}}` and `Qs::AbstractArray{SMatrix{…}}`).
 """
-function initial_deformation(ap::AffinePenalty, cs, Qs)
+function initial_deformation(ap::AffinePenalty, cs, Qs; λt=nothing)
     _initial_deformation(ap, cs, Qs)
 end
 
@@ -333,11 +343,20 @@ end
 
 cs2u(::Type{V}, cs) where {V} = V[V((c...,)) for c in cs]
 
-function initial_deformation(ap::AffinePenalty{T,N}, cs::AbstractArray{V}, Qs::AbstractArray{M}) where {T,N,V<:SVector,M<:SMatrix}
+function initial_deformation(ap::AffinePenalty{T,N}, cs::AbstractArray{V}, Qs::AbstractArray{M}; λt=nothing) where {T,N,V<:SVector,M<:SMatrix}
     Tv = eltype(V)
     eltype(M) == Tv || error("element types of cs ($(eltype(V))) and Qs ($(eltype(M))) must match")
     size(M,1) == size(M,2) == length(V) || throw(DimensionMismatch("size $(size(M)) of Qs matrices is inconsistent with cs vectors of size $(size(V))"))
-    _initial_deformation(convert(AffinePenalty{Tv,N}, ap), cs, Qs)
+    if λt === nothing
+        _initial_deformation(convert(AffinePenalty{Tv,N}, ap), cs, Qs)
+    else
+        length(V) == N || throw(DimensionMismatch("Dimensionality $N of ap does not match $(length(V))"))
+        apc = convert(AffinePenalty{Tv,N}, ap)
+        b = prep_b(Tv, cs, Qs)
+        P = TimeHessian(AffineQHessian(apc, Qs, identity), convert(Tv, λt))
+        x, isconverged = find_opt(P, b)
+        convert_to_fixed(SVector{N,Tv}, x, size(cs)), isconverged
+    end
 end
 
 function to_full(ap::AffinePenalty{T,N}, Qs) where {T,N}
@@ -459,105 +478,6 @@ function _affine_part!(g, ap::AffinePenalty{T,N}, u) where {T,N}
     s
 end
 
-function initial_deformation(ap::AffinePenalty{T,N}, cs, Qs, ϕ_old, maxshift) where {T,N}
-    error("This is broken, don't use it")
-    b = prep_b(T, cs, Qs)
-    # In case the grid is really big, solve iteratively
-    # (The matrix is not sparse, but matrix-vector products can be
-    # computed efficiently.)
-    P0 = AffineQHessian(ap, Qs, identity)
-    x0 = find_opt(P0, b)
-    P = AffineQHessian(ap, Qs, ϕ_old)
-    x = find_opt(P, b, maxshift, x0)
-    u0 = convert_to_fixed(x, (N,size(cs)...)) #reinterpret(SVector{N,eltype(x)}, x, size(cs))
-end
-
-# type for minimization with composition (which turns the problem into
-# a nonlinear problem)
-mutable struct InitialDefOpt{AQH,B} <: GradOnlyBoundsOnly
-    P::AQH
-    b::B
-end
-
-function find_opt(P::AffineQHessian{AP,M,N,Φ}, b, maxshift, x0) where {AP,M,N,Φ<:GridDeformation}
-    objective = InitialDefOpt(P, b)
-#=
-    solver = IpoptSolver(hessian_approximation="limited-memory",
-                         print_level=0,
-                         sb="yes")
-    m = MOI.NonlinearModel(solver)
-    T = eltype(b)
-    n = length(b)
-    ub1 = T[maxshift...] - T(RegisterFit.register_half)
-    ub = repeat(ub1, outer=[div(n, length(maxshift))])
-    MOI.loadproblem!(m, n, 0, -ub, ub, T[], T[], :Min, objective)
-    MOI.setwarmstart!(m, x0)
-    MOI.optimize!(m)
-    stat = MOI.status(m)
-    stat == :Optimal || @warn("Solution was not optimal")
-    MOI.getsolution(m)
-=#
-    T = eltype(b)
-    n = length(b)
-
-    model = Model(optimizer_with_attributes(Ipopt.Optimizer,
-                                            "hessian_approximation"=>"limited-memory",
-                                            "print_level" => 0,
-                                            "sb" => "yes"))
-    ub1 = T[maxshift...] .- T(RegisterFit.register_half)
-    ub = repeat(ub1, outer=[div(n, length(maxshift))])
-    @variable(model, -ub[i] <= x[i in 1:n] <= ub[i], start = x0[i])
-    @operator(model, op_objective, n, (x...) -> MOI.eval_objective(objective, collect(x)),
-                            (g, x...) -> MOI.eval_objective_gradient(objective, g, collect(x)))
-    @objective(model, Min, op_objective(x...))
-    fval0 =  MOI.eval_objective(objective, x0)
-    isfinite(fval0) || error("Initial value must be finite")
-    JuMP.optimize!(model)
-
-    stat = termination_status(model)
-    stat == LOCALLY_SOLVED || @warn("Solution was not optimal")
-    JuMP.value.(x)
-end
-
-# We omit the constant term ∑_i cs[i]'*Qs[i]*cs[i], since it won't
-# affect the solution
-MOI.eval_objective(d::InitialDefOpt, x::AbstractVector) =
-    _eval_f(d.P, d.b, x)
-
-function _eval_f(P::AffineQHessian{AffinePenalty{T,N}}, b, x::AbstractVector) where {T,N}
-    gridsize = size(P.Qs)
-    n = prod(gridsize)
-    u  = convert_to_fixed(x, (N,gridsize...))# reinterpret(SVector{N,T}, x, gridsize)
-    bf = convert_to_fixed(b, (N,gridsize...))# reinterpret(SVector{N,T}, b, gridsize)
-    λ = P.ap.λ
-    P.ap.λ = λ*n/2
-    val = affine_part!(nothing, P, u)
-    P.ap.λ = λ
-    for i = 1:n
-        val += ((u[i]' * P.Qs[i] * u[i])/2 - bf[i]'*u[i])[1]
-    end
-    val
-end
-
-function MOI.eval_objective_gradient(d::InitialDefOpt, grad_f, x)
-    P, b = d.P, d.b
-    copy!(grad_f, P*x-b)
-end
-
-function affine_part!(g, P::AffineQHessian{AP,M,N,Φ}, u) where {AP,M,N,Φ<:GridDeformation}
-    ϕ_c, g_c = compose(P.ϕ_old, GridDeformation(u, P.ϕ_old.nodes))
-    penalty!(g, P.ap, ϕ_c, g_c)
-end
-
-function affine_part!(::Nothing, P::AffineQHessian{AP,M,N,Φ}, u) where {AP,M,N,Φ<:GridDeformation}
-    # Sadly, with GradientNumbers this gives an error I haven't traced
-    # down (might be a Julia bug)
-    # ϕ_c = P.ϕ_old(GridDeformation(u, P.ϕ_old.nodes))
-    # penalty!(nothing, P.ap, ϕ_c)
-    u_c = RegisterDeformation._compose(P.ϕ_old.u, u, P.ϕ_old.nodes)
-    penalty!(nothing, P.ap, u_c)
-end
-
 ###
 ### Optimize (via descent) a deformation to mismatch data
 ###
@@ -722,14 +642,15 @@ end
 # end
 
 """
-`ϕs, penalty = optimize!(ϕs, ϕs_old, dp, λt, mmis; kwargs...)`
+`ϕs, penalty = optimize!(ϕs, ϕs_old, dp, mmis; λt=nothing, kwargs...)`
 optimizes a sequence `ϕs` of deformations for an image sequence with
 mismatch data `mmis`, using a spatial deformation penalty `dp` and
-temporal penalty coefficient `λt`.
+temporal penalty coefficient `λt` (default `0`, meaning no temporal coupling).
 """
-function optimize!(ϕs, ϕs_old, dp::AffinePenalty, λt, mmis; kwargs...)
+function optimize!(ϕs::Vector{<:GridDeformation}, ϕs_old, dp::AffinePenalty, mmis; λt=nothing, kwargs...)
     T = eltype(eltype(first(mmis)))
-    objective = DeformTseriesOpt(ϕs, ϕs_old, dp, λt, mmis)
+    λt_val = λt === nothing ? zero(T) : T(λt)
+    objective = DeformTseriesOpt(ϕs, ϕs_old, dp, λt_val, mmis)
     uvec = u_as_vec(ϕs)
     df = OnceDifferentiable(x->MOI.eval_objective(objective, x),
                                 (g,x)->MOI.eval_objective_gradient(objective, g, x),uvec)
@@ -740,11 +661,11 @@ function optimize!(ϕs, ϕs_old, dp::AffinePenalty, λt, mmis; kwargs...)
     _copy!(ϕs, Optim.minimizer(results)), Optim.minimum(results)
 end
 
-function optimize!(ϕs, ϕs_old, dp::AffinePenalty{T,N}, λt, mmis::Array{Tf}) where {Tf<:Number, T, N}
+function optimize!(ϕs::Vector{<:GridDeformation}, ϕs_old, dp::AffinePenalty{T,N}, mmis::Array{Tf}; λt=nothing, kwargs...) where {Tf<:Number, T, N}
     ND = NumDenom{Tf}
     mmisr = reshape(reinterpret(ND, vec(mmis)), tail(size(mmis)))
     mmisc = cachedinterpolators(mmisr, N, ntuple(d->(size(mmisr,d)+1)>>1, N))
-    optimize!(ϕs, ϕs_old, dp, λt, mmisc)
+    optimize!(ϕs, ϕs_old, dp, mmisc; λt, kwargs...)
 end
 
 mutable struct DeformTseriesOpt{D,Dsold,DP,T,M} <: GradOnlyBoundsOnly
@@ -779,61 +700,58 @@ function _copy!(ϕs::Vector{D}, x::Array{T}) where {D<:GridDeformation,T<:Number
 end
 
 """
-`ϕ, penalty = fixed_λ(cs, Qs, nodes, affinepenalty, mmis)` computes an
-optimal deformation `ϕ` 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)` 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!`).
+
+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}`.
 
 See also: `auto_λ`.
 """
-function fixed_λ(cs, Qs, nodes::NTuple{N}, ap::AffinePenalty{T,N}, mmis; ϕs_old = identity, mu_init=0.1, kwargs...) where {T,N}
+function fixed_λ(cs, Qs, nodes::NTuple{N}, ap::AffinePenalty{T,N}, mmis; λt=nothing, ϕs_old=identity, mu_init=0.1, kwargs...) where {T,N}
     maxshift = map(x->(x-1)>>1, size(first(mmis)))
-    u0, isconverged = initial_deformation(ap, cs, Qs)
-    if !isconverged
-        Base.warn_once("initial_deformation failed to converge with λ = ", ap.λ)
-        if any(x->!isfinite(x), convert_from_fixed(u0))
-            u0 = cs2u(SVector{N,T}, cs)
+    if λt === nothing
+        u0, isconverged = initial_deformation(ap, cs, Qs)
+        if !isconverged
+            Base.warn_once("initial_deformation failed to converge with λ = ", ap.λ)
+            if any(x->!isfinite(x), convert_from_fixed(u0))
+                u0 = cs2u(SVector{N,T}, cs)
+            end
         end
+        uclamp!(u0, maxshift)
+        ϕ = GridDeformation(u0, nodes)
+        local mismatch
+        while mu_init > 1e-16
+            ϕ, mismatch, mismatch0 = optimize!(ϕ, ϕs_old, ap, mmis; mu_strategy="monotone", mu_init, kwargs...)
+            mismatch <= mismatch0 && break
+            mu_init /= 10
+            @show mu_init
+        end
+        ϕ, mismatch
+    else
+        λtT = T(λt)
+        apT = convert(AffinePenalty{T,N}, ap)
+        print("Calculating initial guess (this may take a while)...")
+        u0, isconverged = initial_deformation(apT, cs, Qs; λt=λtT)
+        println("done")
+        if !isconverged
+            Base.warn_once("initial_deformation failed to converge with λ = ", ap.λ, ", λt = ", λt)
+        end
+        uclamp!(u0, maxshift)
+        colons = ntuple(ColonFun, Val(N))
+        ϕs = [GridDeformation(u0[colons..., i], nodes) for i = 1:size(u0)[end]]
+        println("Starting optimization.")
+        optimize!(ϕs, ϕs_old, apT, mmis; λt=λtT, kwargs...)
     end
-    uclamp!(u0, maxshift)
-    ϕ = GridDeformation(u0, nodes)
-    local mismatch
-    while mu_init > 1e-16
-        ϕ, mismatch, mismatch0 = optimize!(ϕ, ϕs_old, ap, mmis; mu_strategy="monotone", mu_init=mu_init, kwargs...)
-        mismatch <= mismatch0 && break
-        mu_init /= 10
-        @show mu_init
-    end
-    ϕ, mismatch
-end
-
-"""
-`ϕs, penalty = fixed_λ(cs, Qs, nodes, affinepenalty, λt, mmis)`
-computes an optimal vector-of-deformations `ϕs` for an image sequence,
-using an temporal penalty coefficient `λt`.
-"""
-function fixed_λ(cs::AbstractArray{SVector{N,T}}, Qs::AbstractArray{SMatrix{N,N,T,L}}, nodes::NTuple{N}, ap::AffinePenalty{TP,N}, λt, mmis; ϕs_old = identity, mu_init=0.1, kwargs...) where {T,N,TP,L}
-    λtT = T(λt)
-    apT = convert(AffinePenalty{T,N}, ap)
-    maxshift = map(x->(x-1)>>1, size(first(mmis)))
-    print("Calculating initial guess (this may take a while)...")
-    u0, isconverged = initial_deformation(apT, λtT, cs, Qs)
-    println("done")
-    if !isconverged
-        Base.warn_once("initial_deformation failed to converge with λ = ", ap.λ, ", λt = ", λt)
-    end
-    uclamp!(u0, maxshift)
-    colons = ntuple(ColonFun, Val(N))
-    ϕs = [GridDeformation(u0[colons..., i], nodes) for i = 1:size(u0)[end]]
-    local mismatch
-    println("Starting optimization.")
-    optimize!(ϕs, ϕs_old, apT, λtT, mmis; kwargs...)
 end
 
 # This version re-packs variables as read from the .jld file
-function fixed_λ(cs::Array{Tf}, Qs::Array{Tf}, nodes::NTuple{N}, ap::AffinePenalty{T,N}, λt, mmis::Array{Tf}; kwargs...) where {Tf<:Number,T,N}
+function fixed_λ(cs::Array{Tf}, Qs::Array{Tf}, nodes::NTuple{N}, ap::AffinePenalty{T,N}, mmis::Array{Tf}; λt=nothing, kwargs...) where {Tf<:Number,T,N}
     csr = unsafe_wrap(Array, convert(Ptr{SVector{N,Tf}}, pointer(cs)), (tail(size(cs))...,))
     Qsr = unsafe_wrap(Array, convert(Ptr{similar_type(SArray,Tf,Size(N,N))}, pointer(Qs)), (tail(tail(size(Qs)))...,))
     if length(mmis) > 10^7
@@ -843,7 +761,7 @@ function fixed_λ(cs::Array{Tf}, Qs::Array{Tf}, nodes::NTuple{N}, ap::AffinePena
     ND = NumDenom{Tf}
     mmisr = unsafe_wrap(Array, convert(Ptr{ND}, pointer(mmis)), (tail(size(mmis))...,))
     mmisc = cachedinterpolators(mmisr, N, ntuple(d->(size(mmisr,d)+1)>>1, N))
-    fixed_λ(csr, Qsr, nodes, ap, λt, mmisc; kwargs...)
+    fixed_λ(csr, Qsr, nodes, ap, mmisc; λt, kwargs...)
 end
 
 ###
@@ -877,9 +795,9 @@ 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, `auto_λ(stackindex, cs, Qs,
-nodes, mmis, (λmin, λmax))` will perform the analysis on the chosen
-`stackindex`.
+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_λ`.
@@ -898,27 +816,29 @@ function auto_λ(fixed::AbstractArray{R}, moving::AbstractArray{S}, gridsize::NT
     auto_λ(cs, Qs, nodes, mmis, λrange; kwargs...)
 end
 
-function auto_λ(cs, Qs, nodes::NTuple{N}, mmis, λrange; kwargs...) where N
+function auto_λ(cs, Qs, nodes::NTuple{N}, mmis, λrange; stackidx=nothing, kwargs...) where N
+    if stackidx !== nothing
+        colons = ntuple(d->Colon(), ndims(cs)-1)
+        cs   = cs[  colons..., stackidx]
+        Qs   = Qs[  colons..., stackidx]
+        mmis = mmis[colons..., stackidx]
+    end
     ap = AffinePenalty{Float64,N}(nodes, λrange[1])  # default to affine-residual penalty, Ipopt needs Float64
     auto_λ(cs, Qs, nodes, ap, mmis, λrange; kwargs...)
 end
 
-function auto_λ(stackidx::Integer, cs, Qs, nodes::NTuple{N}, mmis, λrange; kwargs...) where N
-    cs1 = cs[ntuple(d->Colon(),ndims(cs)-1)..., stackidx];
-    Qs1 = Qs[ntuple(d->Colon(),ndims(Qs)-1)..., stackidx];
-    mmis1 = mmis[ntuple(d->Colon(),ndims(mmis)-1)..., stackidx];
-    auto_λ(cs1, Qs1, nodes, mmis1, λrange; kwargs...)
-end
-
-function auto_λ(cs::Array{Tf}, Qs::Array{Tf}, nodes::NTuple{N}, ap::AffinePenalty{T,N}, mmis::Array{Tf}, λrange; kwargs...) where {Tf<:Number,T,N}
+function auto_λ(cs::AbstractArray{<:Number}, Qs::AbstractArray{<:Number}, nodes::NTuple{N}, ap::AffinePenalty{T,N}, mmis::AbstractArray{<:Number}, λrange; kwargs...) where {T,N}
     # Ipopt requires Float64
     auto_λ(convert(Array{Float64}, cs), convert(Array{Float64}, Qs), nodes, ap, convert(Array{Float64}, mmis), λrange; kwargs...)
 end
 
-function auto_λ(cs::Array{Float64}, Qs::Array{Float64}, nodes::NTuple{N}, ap::AffinePenalty{T,N}, mmis::Array{Float64}, λrange; kwargs...) where {T,N}
-    csr = reshape(reinterpret(SVector{N,Float64}, vec(cs)), tail(size(cs)))
-    Qsr = reshape(reinterpret(similar_type(SArray,Float64,Size(N,N)), vec(Qs)), tail(tail(size(Qs))))
-    mmisr = reshape(reinterpret(NumDenom{Float64}, vec(mmis)), tail(size(mmis)))
+function auto_λ(cs::AbstractArray{Float64}, Qs::AbstractArray{Float64}, nodes::NTuple{N}, ap::AffinePenalty{T,N}, mmis::AbstractArray{Float64}, λrange; kwargs...) where {T,N}
+    cs64 = cs isa Array{Float64} ? cs : Array(cs)
+    Qs64 = Qs isa Array{Float64} ? Qs : Array(Qs)
+    mmis64 = mmis isa Array{Float64} ? mmis : Array(mmis)
+    csr = reshape(reinterpret(SVector{N,Float64}, vec(cs64)), tail(size(cs64)))
+    Qsr = reshape(reinterpret(similar_type(SArray,Float64,Size(N,N)), vec(Qs64)), tail(tail(size(Qs64))))
+    mmisr = reshape(reinterpret(NumDenom{Float64}, vec(mmis64)), tail(size(mmis64)))
     mmisc = cachedinterpolators(mmisr, N, ntuple(d->(size(mmisr,d)+1)>>1, N))
     ap64 = convert(AffinePenalty{Float64,N}, ap)
     auto_λ(csr, Qsr, nodes, ap64, mmisc, λrange; kwargs...)
@@ -1025,7 +945,7 @@ function auto_λt(Es, cs, Qs, ap, λtrange)
     λts = Vector{typeof(λt)}(undef, n)
     @showprogress 1 "Calculating quadratic penalty as a function of λt: " for λindex = 1:n
         λts[λindex] = λt
-        u0, isconverged = initial_deformation(ap, λt, cs, Qs)
+        u0, isconverged = initial_deformation(ap, cs, Qs; λt)
         if !isconverged
             println("initial_deformation failed to converge with λ = ", ap.λ, ", λt = ", λt)
         end
@@ -1040,26 +960,10 @@ function auto_λt(Es, cs, Qs, ap, λtrange)
     λts, datapenalty
 end
 
-function auto_λt(Es, cs::Array{Tf}, Qs::Array{Tf}, ap::AffinePenalty{T,N}, λt) where {Tf<:Number,T,N}
-    csr = reshape(reinterpret(SVector{N,Tf}, vec(cs)), tail(size(cs)))
-    Qsr = reshape(reinterpret(similar_type(SArray,Tf,Size(N,N)), vec(Qs)), tail(tail(size(Qs))))
-    auto_λt(Es, csr, Qsr, ap, λt)
-end
 
 ###
 ### Whole-experiment optimization with a temporal roughness penalty
 ###
-function initial_deformation(ap::AffinePenalty{T,N}, λt, cs::AbstractArray{V}, Qs::AbstractArray{M}) where {T,N,V<:SVector,M<:SMatrix}
-    Tv = eltype(V)
-    eltype(M) == Tv || error("element types of cs ($(eltype(V))) and Qs ($(eltype(M))) must match")
-    length(V) == N || throw(DimensionMismatch("Dimensionality $N of ap does not match $(length(V))"))
-    size(M,1) == size(M,2) == N || throw(DimensionMismatch("size $(size(M)) of Qs matrices is inconsistent with cs vectors of size $(size(V))"))
-    apc = convert(AffinePenalty{Tv,N}, ap)
-    b = prep_b(Tv, cs, Qs)
-    P = TimeHessian(AffineQHessian(apc, Qs, identity), convert(Tv, λt))
-    x, isconverged = find_opt(P, b)
-    convert_to_fixed(SVector{N,Tv}, x, size(cs)), isconverged
-end
 
 struct TimeHessian{AQH<:AffineQHessian,T}
     aqh::AQH
@@ -1074,7 +978,7 @@ function (*)(P::TimeHessian{AQH}, x::AbstractVector) where AQH
     yv = P.aqh*x
     ϕs = vec2vecϕ(P.aqh.Qs, x)
     y = convert_to_fixed(SVector{size(P.aqh.Qs[1],1),eltype(yv)}, yv, size(P.aqh.Qs))
-    penalty!(y, P.λt, ϕs)
+    penalty!(y, ϕs, P.λt)
     yv
 end
 
@@ -1083,7 +987,7 @@ function LinearAlgebra.mul!(y::AbstractVector,
                               x::AbstractVector) where AQH
     mul!(y, P.aqh, x)
     ϕs = vec2vecϕ(P.aqh.Qs, x)
-    penalty!(y, P.λt, ϕs)
+    penalty!(y, ϕs, P.λt)
     y
 end
 
@@ -1098,95 +1002,6 @@ end
 end
 
 
-###
-### Mismatch-based optimization of affine transformation
-###
-### NOTE: not updated yet, probably broken
-"""
-`tform = optimize(tform0, mms, nodes)` performs descent-based
-minimization of the total mismatch penalty as a function of the
-parameters of an affine transformation, starting from an initial guess
-`tform0`.  While this is unlikely to yield very accurate results for
-large rotations or skews (the mismatch data are themselves suspect in
-such cases), it can be helpful for polishing small deformations.
-
-For a good initial guess, see `mismatch2affine`.
-"""
-function optimize(tform::AffineMap, mmis, nodes)
-    gridsize = size(mmis)
-    N = length(gridsize)
-    ndims(tform) == N || error("Dimensionality of tform is $(ndims(tform)), which does not match $N for nums/denoms")
-    mm = first(mmis)
-    mxs = maxshift(mm)
-    T = eltype(eltype(mm))
-    # Compute the bounds
-    asz = arraysize(nodes)
-    center = T[(asz[i]+1)/2 for i = 1:N]
-    X = zeros(T, N+1, prod(gridsize))
-    for (i, node) in enumerate(eachnode(nodes))
-        X[1:N,i] = node - center
-        X[N+1,i] = 1
-    end
-    bound = convert(Vector{T}, [mxs .- register_half; Inf])
-    lower = repeat(-bound, outer=[1,size(X,2)])
-    upper = repeat( bound, outer=[1,size(X,2)])
-    # Extract the parameters from the initial guess
-    Si = tform.linear
-    displacement = tform.translation
-    A = convert(Matrix{T}, [Si-Matrix{Float64}(I,N,N) displacement; zeros(1,N) 1])
-    # Determine the blocks that start in-bounds
-    AX = A*X
-    keep = trues(gridsize)
-    for j = 1:length(keep)
-        for idim = 1:N
-            xi = AX[idim,j]
-            if xi < -mxs[idim]+register_half_safe || xi > mxs[idim]-register_half_safe
-                keep[j] = false
-                break
-            end
-        end
-    end
-    if !any(keep)
-        @show tform
-        @warn("No valid blocks were found")
-        return tform
-    end
-    ignore = !keep[:]
-    lower[:,ignore] = -Inf
-    upper[:,ignore] =  Inf
-    # Assemble the objective and constraints
-
-    constraints = Optim.ConstraintsL(X', lower', upper')
-    gtmp = Array{SVector{N,T}}(undef, gridsize)
-    objective = (x,g) -> affinepenalty!(g, x, mmis, keep, X', gridsize, gtmp)
-    @assert typeof(objective(A', T[])) == T
-    result = interior(DifferentiableFunction(x->objective(x,T[]), Optim.dummy_g!, objective), A', constraints, method=:cg)
-    @assert Optim.converged(result)
-    Aopt = result.minimum'
-    Siopt = Aopt[1:N,1:N] + Matrix{Float64}(I,N,N)
-    displacementopt = Aopt[1:N,end]
-    AffineMap(convert(Matrix{T}, Siopt), convert(Vector{T}, displacementopt)), result.f_minimum
-end
-
-function affinepenalty!(g, At, mmis, keep, Xt, gridsize::NTuple{N}, gtmp) where N
-    u = _calculate_u(At, Xt, gridsize)
-    @assert eltype(u) == eltype(At)
-    val = penalty!(gtmp, u, mmis, keep)
-    @assert isa(val, eltype(At))
-    if !isempty(g)
-        T = eltype(eltype(gtmp))
-        nblocks = size(Xt,1)
-        At_mul_Bt!(g, Xt, [reshape(reinterpret(T,vec(gtmp)),(N,nblocks)); zeros(1,nblocks)])
-    end
-    val
-end
-
-function _calculate_u(At, Xt, gridsize::NTuple{N}) where N
-    Ut = Xt*At
-    u = Ut[:,1:size(Ut,2)-1]'                   # discard the dummy dimension
-    reshape(reinterpret(SVector{N, eltype(u)}, vec(u)), gridsize) # put u in the shape of the grid
-end
-
 ###
 ### Fitting to a sigmoid
 ###
diff --git a/test/register_optimize.jl b/test/register_optimize.jl
index e5f4c3b..125fe34 100644
--- a/test/register_optimize.jl
+++ b/test/register_optimize.jl
@@ -153,7 +153,7 @@ end
     v = zeros(size(P,1)); v[end] = 1
     @test ≈(P*v, vec(Ac[end,:]), atol=0.0001)
     ux = initial_guess_direct(A, 1.0, csr, Qsr)
-    u, isconverged = RegisterOptimize.initial_deformation(ap, 1.0, cs, Qs)
+    u, isconverged = RegisterOptimize.initial_deformation(ap, cs, Qs; λt=1.0)
     @test isconverged
     @test size(u) == size(ux)
     @test eltype(u) == SVector{2,Float64}
@@ -273,14 +273,14 @@ end
     cs = cat([5,-3], [0,0], [3,-1], dims=2)
     gridsize = (2,2)
     denom = ones(15,15)
-    mms = tighten([quadratic(cs[:,t], Qs[:,:,t], denom) for i = 1:gridsize[1], j = 1:gridsize[2], t = 1:3])
+    mms = tighten([quadratic(denom, cs[:,t], Qs[:,:,t]) for i = 1:gridsize[1], j = 1:gridsize[2], t = 1:3])
     mmis = RegisterPenalty.interpolate_mm!(mms)
     nodes = (range(1, stop=100, length=gridsize[1]), range(1, stop=99, length=gridsize[2]))
     ap = RegisterPenalty.AffinePenalty(nodes, 1.0)
     u = randn(2, gridsize..., 3)
     ϕs = tighten([GridDeformation(u[:,:,:,t], nodes) for t = 1:3])
     g = similar(u)
-    ϕs, fval = RegisterOptimize.optimize!(ϕs, identity, ap, 1.0, mmis)
+    ϕs, fval = RegisterOptimize.optimize!(ϕs, identity, ap, mmis; λt=1.0)
     c = 1/prod(gridsize)  # not sure about this
     A = [2c+1 -1 0; -1 2 -1; 0 -1 2c+1]
     target = (A\(2c*cs'))'