diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 228f9f9..f1a0057 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -43,4 +43,36 @@ jobs:
       - uses: codecov/codecov-action@v5
         with:
           file: lcov.info
-          token: ${{ secrets.CODECOV_TOKEN }}
\ No newline at end of file
+          token: ${{ secrets.CODECOV_TOKEN }}
+  docs:
+    name: Documentation
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+      statuses: write
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: '1'
+      - uses: actions/cache@v4
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-docs-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-docs-${{ env.cache-name }}-
+            ${{ runner.os }}-docs-
+            ${{ runner.os }}-
+      - name: registry_add
+        run: julia -e 'using Pkg; pkg"registry add General https://github.com/HolyLab/HolyLabRegistry.git"'
+      - name: Install dependencies
+        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy docs
+        uses: julia-actions/julia-docdeploy@v1
+        with:
+          install-package: false
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 9f002e4..c2bc4b3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "RegisterFit"
 uuid = "36121b08-3789-5198-aff2-59a3443d9b59"
+version = "1.0.0"
 authors = ["Tim Holy <tim.holy@gmail.com>"]
-version = "0.2.2"
 
 [deps]
 CenterIndexedArrays = "46a7138f-0d70-54e1-8ada-fb8296f91f24"
@@ -17,8 +17,11 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
+Aqua = "0.8"
 CenterIndexedArrays = "0.0, 0.1, 0.2, 1"
 CoordinateTransformations = "0.5, 0.6"
+Documenter = "1"
+ExplicitImports = "1"
 ImageBase = "0.1"
 ImageTransformations = "0.10"
 Interpolations = "0.12, 0.13, 0.14, 0.15, 0.16"
@@ -34,9 +37,12 @@ Test = "1"
 julia = "1.10"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e"
 ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ImageBase", "ImageTransformations", "Test"]
+test = ["Aqua", "Documenter", "ExplicitImports", "ImageBase", "ImageTransformations", "Test"]
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..2023059
--- /dev/null
+++ b/README.md
@@ -0,0 +1,69 @@
+# RegisterFit.jl
+
+[![CI](https://github.com/HolyLab/RegisterFit.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/RegisterFit.jl/actions/workflows/CI.yml)
+[![Coverage](https://codecov.io/gh/HolyLab/RegisterFit.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HolyLab/RegisterFit.jl)
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://HolyLab.github.io/RegisterFit.jl/stable)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://HolyLab.github.io/RegisterFit.jl/dev)
+
+RegisterFit computes affine transformations that minimize image registration
+mismatch. Given per-aperture mismatch data from
+[RegisterMismatch](https://github.com/HolyLab/RegisterMismatch.jl), it finds
+the best-fit affine transform by fitting each aperture's mismatch to a quadratic
+form and solving a global least-squares problem. It is part of the
+[HolyLab image registration pipeline](https://github.com/HolyLab).
+
+## Installation
+
+RegisterFit is distributed through the
+[HolyLabRegistry](https://github.com/HolyLab/HolyLabRegistry). Add the registry
+once, then install the package:
+
+```julia
+using Pkg
+pkg"registry add https://github.com/HolyLab/HolyLabRegistry.git"
+Pkg.add("RegisterFit")
+```
+
+## Quick start
+
+### Fitting a single aperture's mismatch to a quadratic
+
+```julia
+using RegisterFit, RegisterCore
+
+# Synthetic mismatch: parabola centred at shift (1, -2)
+num = [(i - 1)^2 + (j + 2)^2 for i in -5:5, j in -5:5]
+mm  = MismatchArray(num, ones(11, 11))
+
+E0, u0, Q = qfit(mm, 0.5)
+# E0 ≈ 0.0      — mismatch value at the minimum
+# u0 ≈ [1, -2]  — shift at the minimum
+# Q  ≈ I        — curvature matrix
+```
+
+### Rigid alignment via the Principal Axes Transform
+
+```julia
+using RegisterFit
+
+# Horizontal bar in fixed, vertical bar in moving (≈ 90° rotation)
+fixed  = zeros(5, 7); fixed[3, 2:6]  .= 1.0
+moving = zeros(7, 5); moving[2:6, 3] .= 1.0
+
+tfms = pat_rotation(fixed, moving)   # 2 candidate AffineMap transforms in 2D
+# Evaluate each candidate against mismatch data and pick the best
+```
+
+## Overview
+
+The typical workflow in the registration pipeline is:
+
+1. **Compute mismatch** — per-aperture `MismatchArray`s from `RegisterMismatch`
+2. **`mms2fit!`** — fit each aperture to a quadratic; returns per-aperture shifts
+   and curvature matrices
+3. **`mismatch2affine`** — global least-squares solve over all apertures to
+   produce a single `AffineMap`
+4. **Refine** — further optimisation in `RegisterDeformation`
+
+See the [documentation](https://HolyLab.github.io/RegisterFit.jl/stable) for
+the full API reference.
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..ba08347
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+RegisterCore = "67712758-55e7-5c3c-8e85-dda1d7758434"
+RegisterFit = "36121b08-3789-5198-aff2-59a3443d9b59"
+
+[compat]
+Documenter = "1"
+julia = "1.10"
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..8c4eda5
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,23 @@
+using Documenter
+using RegisterFit
+using RegisterCore
+
+DocMeta.setdocmeta!(RegisterFit, :DocTestSetup, :(using RegisterFit, RegisterCore); recursive=true)
+
+makedocs(;
+    modules=[RegisterFit],
+    sitename="RegisterFit.jl",
+    format=Documenter.HTML(;
+        canonical="https://HolyLab.github.io/RegisterFit.jl",
+    ),
+    checkdocs=:exports,
+    pages=[
+        "Home" => "index.md",
+        "API Reference" => "api.md",
+    ],
+)
+
+deploydocs(;
+    repo="github.com/HolyLab/RegisterFit.jl",
+    devbranch="master",
+)
diff --git a/docs/src/api.md b/docs/src/api.md
new file mode 100644
index 0000000..2990c39
--- /dev/null
+++ b/docs/src/api.md
@@ -0,0 +1,25 @@
+# API Reference
+
+## Global optimization
+
+```@docs
+mismatch2affine
+pat_rotation
+optimize_per_aperture
+```
+
+## Quadratic fitting
+
+```@docs
+qfit
+mms2fit!
+qbuild
+```
+
+## Displacement utilities
+
+```@docs
+uisvalid
+uclamp!
+principalaxes
+```
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..1b3f1a5
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,105 @@
+# RegisterFit.jl
+
+```@docs
+RegisterFit.RegisterFit
+```
+
+RegisterFit computes affine transformations that minimize image registration
+mismatch. It is part of the [HolyLab](https://github.com/HolyLab) image
+registration pipeline.
+
+## Installation
+
+RegisterFit is distributed through the
+[HolyLabRegistry](https://github.com/HolyLab/HolyLabRegistry). Add the registry
+once, then install:
+
+```julia
+using Pkg
+pkg"registry add https://github.com/HolyLab/HolyLabRegistry.git"
+Pkg.add("RegisterFit")
+```
+
+## Concepts
+
+### Apertures and mismatch data
+
+The registration pipeline divides an image into overlapping sub-regions called
+**apertures**. For each aperture, `RegisterMismatch` computes a `MismatchArray`
+that records the normalized squared difference between the fixed and moving images
+as a function of shift. Each `MismatchArray` element stores a numerator/denominator
+pair; only locations where `denom > thresh` are treated as valid.
+
+### Quadratic fitting
+
+Within a single aperture, the mismatch surface is often well approximated by a
+quadratic:
+
+```
+mm(u) ≈ E0 + (u − u0)ᵀ Q (u − u0)
+```
+
+[`qfit`](@ref) finds `E0` (minimum mismatch value), `u0` (shift at the minimum),
+and `Q` (curvature matrix) for one aperture. [`mms2fit!`](@ref) does this for an
+entire grid of apertures at once.
+
+### Global affine solve
+
+With one quadratic per aperture, [`mismatch2affine`](@ref) assembles a global
+linear system whose solution is the affine transformation (rotation + shear +
+translation) that best explains all aperture displacements simultaneously. The
+result is an `AffineMap` from
+[CoordinateTransformations.jl](https://github.com/JuliaGeometry/CoordinateTransformations.jl).
+
+### Principal Axes Transform
+
+[`pat_rotation`](@ref) provides a complementary rigid-alignment approach that
+matches the intensity-weighted covariance ellipsoids of two images. It is useful
+as an initial guess before running the quadratic optimization. Because an ellipsoid
+is symmetric, there are 2 candidate rotations in 2D and 4 in 3D; you must evaluate
+each against the mismatch data to pick the best one.
+
+## Typical workflow
+
+```julia
+using RegisterFit, RegisterCore
+
+# 1. Obtain per-aperture mismatch arrays from RegisterMismatch (not shown)
+#    mms :: Matrix{MismatchArray}  (gridsize matches spatial dims)
+
+# 2. Fit quadratics and prepare for interpolation
+cs, Qs, mmis = mms2fit!(mms, thresh)
+
+# 3. Solve for the global affine transform
+tform = mismatch2affine(mms, thresh, knots)
+
+# 4. (Optional) rigid pre-alignment with PAT
+candidates = pat_rotation(fixed, moving)
+# … pick best candidate, then refine with mismatch2affine / RegisterDeformation
+```
+
+## Quick examples
+
+### Fitting one aperture
+
+```julia
+using RegisterFit, RegisterCore
+
+num = [(i - 1)^2 + (j + 2)^2 for i in -5:5, j in -5:5]
+mm  = MismatchArray(num, ones(11, 11))
+
+E0, u0, Q = qfit(mm, 0.5)
+# E0 ≈ 0.0, u0 ≈ [1.0, -2.0], Q ≈ I
+```
+
+### Rigid alignment via Principal Axes Transform
+
+```julia
+using RegisterFit
+
+fixed  = zeros(5, 7); fixed[3, 2:6]  .= 1.0   # horizontal bar
+moving = zeros(7, 5); moving[2:6, 3] .= 1.0   # vertical bar
+
+tfms = pat_rotation(fixed, moving)   # 2 candidate AffineMap transforms
+# Select the candidate with the smallest mismatch
+```
diff --git a/src/RegisterFit.jl b/src/RegisterFit.jl
index 30b3b0d..fc483e6 100644
--- a/src/RegisterFit.jl
+++ b/src/RegisterFit.jl
@@ -1,10 +1,16 @@
 module RegisterFit
 
-using Interpolations, StaticArrays, Optim, CoordinateTransformations, NLsolve
-using Statistics, LinearAlgebra
-using RegisterPenalty, RegisterCore, CenterIndexedArrays
-
-using Base: @nloops, @nexprs, @nref, @nif
+using CenterIndexedArrays: CenterIndexedArrays, CenterIndexedArray
+using CoordinateTransformations: CoordinateTransformations, AffineMap
+using Interpolations: Interpolations
+using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, I, cholesky, det, diag, dot, eigen, mul!, svd
+using NLsolve: NLsolve, nlsolve
+using Optim: Optim
+using RegisterCore: RegisterCore, MismatchArray, NumDenom, argmin_mismatch, maxshift
+using RegisterPenalty: RegisterPenalty, interpolate_mm!
+using StaticArrays: StaticArrays, SArray, SVector, Size, StaticVector, similar_type
+using Statistics: Statistics, mean
+import Base.Cartesian: @nloops, @nexprs, @nref, @nif
 
 export
     mismatch2affine,
@@ -18,25 +24,21 @@ export
     uclamp!
 
 """
-# RegisterFit
-
-This module contains a number of functions that calculate affine
-transformations that minimize mismatch.  The functions are organized
-into categories:
-
-### Global optimization
+RegisterFit provides functions that compute affine transformations minimizing
+image registration mismatch, given per-aperture mismatch data from `RegisterMismatch`.
 
-- `mismatch2affine`: a transformation computed from mismatch data by least squares
-- `pat_rotation`: find the optimal rigid transform via a Principal Axes Transformation
-- `optimize_per_aperture`: naive registration using independent apertures
+## Global optimization
 
-### Utilities
+- [`mismatch2affine`](@ref): affine transform from mismatch data by least squares
+- [`pat_rotation`](@ref): rigid alignment via a Principal Axes Transformation
+- [`optimize_per_aperture`](@ref): naive per-aperture displacement search
 
-- `qfit`: fit a single aperture's mismatch data to a quadratic form
-- `mms2fit!`: prepare an array-of-mismatcharrays for optimization
-- `qbuild`: reconstruct mismatch data from a quadratic form
-- `uclamp!` and `uisvalid!`: check/enforce bounds on optimization
+## Utilities
 
+- [`qfit`](@ref): fit a single aperture's mismatch to a quadratic form
+- [`mms2fit!`](@ref): prepare an array of mismatch arrays for optimization
+- [`qbuild`](@ref): reconstruct mismatch data from a quadratic form
+- [`uclamp!`](@ref) and [`uisvalid`](@ref): enforce/check bounds on displacements
 """
 RegisterFit
 
@@ -45,20 +47,26 @@ const register_half = 0.5001
 const register_half_safe = 0.51
 
 """
-`tform = mismatch2affine(mms, thresh, knots)` returns an affine
-transformation that is a "best initial guess" for the transformation
-that would minimize the mismatch.  The mismatch is encoded in `mms`
-(of the format returned by RegisterMismatch), and `thresh` is the
-denominator-threshold that determines regions that have sufficient
-pixel/intensity overlap to be considered valid.  `knots` represents
-the aperture centers (see RegisterDeformation).
-
-The algorithm is based on fitting each aperture to a quadratic, and
-then performing a least-squares minimization of the
-sum-over-apertures.  The strength of this procedure is that it finds a
-global solution; however, it is based on a coarse approximation, the
-quadratic fit.  If you want to polish `tform` without relying on the
-quadratic fit, see `optimize`.
+    tform = mismatch2affine(mms, thresh, knots)
+
+Return an `AffineMap` that is a least-squares "best initial guess" for the
+transformation minimizing the mismatch.
+
+`mms` is an array of `MismatchArray`s (one per aperture, in the format returned
+by `RegisterMismatch`). `thresh` is the denominator threshold that determines
+which aperture regions have sufficient pixel/intensity overlap to be valid.
+`knots` specifies the aperture centers (see `RegisterDeformation`).
+
+The algorithm fits each aperture to a quadratic, then solves a global least-squares
+problem over all apertures — guaranteeing a global solution at the cost of the
+quadratic approximation. If `thresh` is too restrictive, it is halved up to three
+times before raising an error.
+
+# Returns
+- `tform::AffineMap` — affine transformation (linear map + translation)
+
+To refine `tform` beyond the quadratic approximation, see `optimize` in the
+registration pipeline.
 """
 function mismatch2affine(mms, thresh, knots)
     gridsize = size(mms)
@@ -137,22 +145,45 @@ end
 
 
 """
-`u = optimize_per_aperture(mms, thresh)` computes the "naive"
-displacement in each aperture to minimize the mismatch. Each aperture
-is examined independently of all others. `thresh` establishes a
-threshold for the mismatch data.
+    u = optimize_per_aperture(mms, thresh)
+
+Compute the naive per-aperture displacement that minimizes the mismatch, treating
+each aperture independently. `mms` is a `Vector` of `MismatchArray`s (one per
+aperture) and `thresh` is the denominator threshold.
+
+For each aperture, the first shift-dimension component of the `argmin` is
+recorded. The returned array has size `(1, length(mms))`.
 
-See also `indmin_mismatch`.
+See also `RegisterCore.argmin_mismatch`.
+
+# Returns
+- `u::Matrix{Float64}` of size `(1, n)` — first shift component at the minimum
+  for each of the `n` apertures
+
+# Examples
+```jldoctest
+julia> using RegisterCore
+
+julia> num1 = [(i - 1)^2 + j^2 for i in -5:5, j in -5:5];
+
+julia> num2 = [i^2 + j^2 for i in -5:5, j in -5:5];
+
+julia> mms = [MismatchArray(num1, ones(11, 11)), MismatchArray(num2, ones(11, 11))];
+
+julia> optimize_per_aperture(mms, 0.5)
+1×2 Matrix{Float64}:
+ 1.0  0.0
+```
 """
 function optimize_per_aperture(mms, thresh)
     gridsize = size(mms)
     nd = length(gridsize)
     u = zeros(nd, gridsize...)
     utmp = zeros(nd)
-    for (iblock, mm) in enumerate(mms)
-        I = indmin_mismatch(mm, thresh)
-        for idim in 1:nd
-            u[idim, iblock] = I[idim]
+    for (iblock,mm) in enumerate(mms)
+        I = argmin_mismatch(mm, thresh)
+        for idim = 1:nd
+            u[idim,iblock] = I[idim]
         end
     end
     return u
@@ -160,9 +191,39 @@ end
 
 
 """
-`r = qbuild(E0, u0, Q, maxshift)` builds an estimate of the mismatch
-ratio given the quadratic form parameters `E0, u0, Q` obtained from
-`qfit`.  Often useful for debugging or visualization.
+    r = qbuild(E0, umin, Q, maxshift)
+
+Build a `CenterIndexedArray` representing the quadratic mismatch approximation
+over the full shift domain `[-maxshift[d], maxshift[d]]`. The quadratic model is:
+
+```
+    r[u] = E0 + (u - umin)' * Q * (u - umin)
+```
+
+`E0`, `umin`, and `Q` are the outputs of [`qfit`](@ref). Useful for debugging
+and visualizing the quadratic fit.
+
+# Returns
+- `r::CenterIndexedArray` — evaluated mismatch, indexed from `-maxshift` to `+maxshift`
+
+# Examples
+```jldoctest
+julia> using RegisterCore
+
+julia> num = [(i - 1)^2 + (j + 2)^2 for i in -5:5, j in -5:5];
+
+julia> mm = MismatchArray(num, ones(11, 11));
+
+julia> E0, umin, Q = qfit(mm, 0.5);
+
+julia> r = qbuild(E0, umin, Q, (5, 5));
+
+julia> r[1, -2]
+0.0
+
+julia> r[0, 0]
+5.0
+```
 """
 function qbuild(E0::Real, umin::Vector, Q::Matrix, maxshift)
     d = length(maxshift)
@@ -184,12 +245,25 @@ function qbuild(E0::Real, umin::Vector, Q::Matrix, maxshift)
 end
 
 """
-`tf = uisvalid(u, maxshift)` returns `true` if all entries of `u` are
-within the allowed domain.
+    tf = uisvalid(u, maxshift)
+
+Return `true` if every entry of the displacement array `u` is within the allowed
+domain: `|u[idim, j]| < maxshift[idim] - 0.5001` for all aperture positions `j`
+and displacement dimensions `idim`. Returns `false` as soon as any entry violates
+this condition.
+
+# Examples
+```jldoctest
+julia> uisvalid([1.5, 0.5], (3, 3))
+true
+
+julia> uisvalid([2.5, 0.5], (3, 3))
+false
+```
 """
 function uisvalid(u::AbstractArray{T}, maxshift) where {T <: Number}
     nd = size(u, 1)
-    sztail = Base.tail(size(u))
+    sztail = size(u)[2:end]
     for j in CartesianIndices(sztail), idim in 1:nd
         if abs(u[idim, j]) >= maxshift[idim] - register_half
             return false
@@ -199,11 +273,28 @@ function uisvalid(u::AbstractArray{T}, maxshift) where {T <: Number}
 end
 
 """
-`u = uclamp!(u, maxshift)` clamps the values of `u` to the allowed domain.
+    uclamp!(u, maxshift)
+
+Clamp the entries of the displacement array `u` in-place so that each satisfies
+`|u[idim, j]| ≤ maxshift[idim] - 0.51`. Returns `u`.
+
+Accepts both numeric arrays of shape `(nd, apertures...)` and arrays whose elements
+are `StaticVector`s (e.g., `Array{SVector{N,T}}`); both representations are mutated
+in-place.
+
+# Examples
+```jldoctest
+julia> u = [4.0, -5.0];
+
+julia> uclamp!(u, (3, 3))
+2-element Vector{Float64}:
+  2.49
+ -2.49
+```
 """
 function uclamp!(u::AbstractArray{T}, maxshift) where {T <: Number}
     nd = size(u, 1)
-    sztail = Base.tail(size(u))
+    sztail = size(u)[2:end]
     for j in CartesianIndices(sztail), idim in 1:nd
         u[idim, j] = max(-maxshift[idim] + register_half_safe, min(u[idim, j], maxshift[idim] - register_half_safe))
     end
@@ -216,9 +307,31 @@ function uclamp!(u::AbstractArray{T}, maxshift) where {T <: StaticVector}
 end
 
 """
-`center, cov = principalaxes(img)` computes the principal axes of an
-image `img`.  `center` is the centroid of intensity, and `cov` the
-covariance matrix of the intensity.
+    center, cov = principalaxes(img)
+
+Compute the intensity-weighted centroid and covariance of image `img`.
+Coordinates are 1-based array indices. `NaN` pixels are ignored.
+
+# Returns
+- `center::Vector{T}` — intensity-weighted centroid, length `ndims(img)`
+- `cov::Matrix{T}` — `N×N` intensity-weighted covariance matrix
+
+# Examples
+```jldoctest
+julia> img = zeros(5, 5); img[3, 3] = 1.0;
+
+julia> center, cov = principalaxes(img);
+
+julia> center
+2-element Vector{Float64}:
+ 3.0
+ 3.0
+
+julia> cov
+2×2 Matrix{Float64}:
+ 0.0  0.0
+ 0.0  0.0
+```
 """
 function principalaxes(img::AbstractArray{T, N}) where {T, N}
     Ts = typeof(zero(T) / 1)
@@ -266,23 +379,43 @@ end
 end
 
 """
-`tfms = pat_rotation(fixed, moving, [SD=eye])` computes the Principal
-Axes Transform aligning the low-order moments of two images. The
-reference image is `fixed`, and `moving` is the raw moving image that
-you wish to align to `fixed`.  `SD` is a "spacedimensions" matrix, in
-some cases needed to ensure that rotations in *physical* space
-correspond to orthogonal matrices in array-index units.  For example,
-if your axes are not uniformly sampled, `SD = diagm(voxelspacing)`.
-
-If you're aligning many images to `fixed`, you may alternatively call
-this as `tfms = pat_rotation(fixedpa, moving, [SD=eye])`.  `fixedpa`
-is a `(center,cov)` tuple obtained from `principalaxes(fixed)`.
-
-`tfms` is a list of potential AffineTransform candidates.  PA data,
-being based on ellipsoids, are ambiguous up to rotations by 180
-degrees (i.e., sign-flips of even numbers of coordinates).
-Consequently, you may need to check all of the candidates for the one
-that produces the best alignment.
+    tfms = pat_rotation(fixed, moving)
+    tfms = pat_rotation(fixed, moving, SD)
+    tfms = pat_rotation(fixedpa, moving)
+    tfms = pat_rotation(fixedpa, moving, SD)
+
+Compute the Principal Axes Transform (PAT) aligning the low-order intensity
+moments of two images. `fixed` is the reference image and `moving` is the image
+to align. `fixedpa` is a `(center, cov)` tuple from [`principalaxes`](@ref),
+useful when aligning many images to the same reference to avoid recomputing its
+principal axes.
+
+`SD` is an optional spatial-dimensions matrix that accounts for non-isotropic
+sampling (e.g., `SD = Diagonal(voxelspacing)`). Defaults to the identity.
+
+Because intensity ellipsoids are ambiguous up to 180° rotations (sign-flips of
+an even number of coordinate axes), the function returns multiple candidate
+transforms. Evaluate alignment quality for each candidate and select the best.
+
+# Returns
+- `tfms::Vector{AffineMap}` — 2 candidates in 2D, 4 candidates in 3D
+
+# Examples
+```jldoctest
+julia> fixed = zeros(5, 7); fixed[3, 2:6] .= 1.0;   # horizontal bar
+
+julia> moving = zeros(7, 5); moving[2:6, 3] .= 1.0;  # vertical bar
+
+julia> tfms = pat_rotation(fixed, moving);
+
+julia> length(tfms)
+2
+
+julia> tfms[1].linear   # ≈ 90° rotation
+2×2 Matrix{Float64}:
+  0.0  1.0
+ -1.0  0.0
+```
 """
 function pat_rotation(
         fixedmoments::Tuple{Vector, Matrix}, moving::AbstractArray,
@@ -330,7 +463,7 @@ function pat_rotation(
     return tfms
 end
 
-pat_rotation(fixed::AbstractArray, moving::AbstractArray, SD = eye(ndims(fixed))) =
+pat_rotation(fixed::AbstractArray, moving::AbstractArray, SD = Matrix{Float64}(I, ndims(fixed), ndims(fixed))) =
     pat_rotation(principalaxes(fixed), moving, SD)
 
 function pat_at(S, SD, c, fmean, mmean)
@@ -378,21 +511,49 @@ end
 end
 
 """
-`E0, u0, Q = qfit(mm, thresh; [maxsep=size(mm), opt=true])` performs a
-quadratic fit of the mismatch data in `mm`.  On output, `u0` and `E0`
-hold the position and value, respectively, of the shift with smallest
-mismatch, and `Q` is a matrix representing the best fit to a model
+    E0, u0, Q = qfit(mm, thresh; maxsep=size(mm), opt=true)
+
+Perform a quadratic fit of the mismatch data in `mm`. Returns the mismatch value
+`E0` and shift `u0` at the minimum, plus the curvature matrix `Q` of the
+best-fit model:
 
 ```
-    mm ≈ E0 + (u-u0)'*Q*(u-u0)
+    mm ≈ E0 + (u - u0)' * Q * (u - u0)
 ```
-Only those shift-locations with `mm[i].denom > thresh` are used in
-performing the fit.
 
-`maxsep` allows you to restrict the fit to a region where each
-coordinate satisfies `|u[d]-u0[d]| <= maxsep[d]`. If `opt` is false,
-`Q` is a heuristic estimate of the best-fit `Q`. This can boost
-performance at the cost of accuracy.
+Only shift-locations where `mm[i].denom > thresh` are used. If no valid locations
+exist, returns `(zero(T), zeros(T, d), zeros(T, d, d))`.
+
+`maxsep` restricts the fit to shifts satisfying `|u[d] - u0[d]| ≤ maxsep[d]`.
+Setting `opt=false` uses a fast heuristic for `Q` instead of a full nonlinear
+solve, trading accuracy for speed.
+
+# Returns
+- `E0::T` — mismatch value at the fitted minimum
+- `u0::Vector{T}` — shift coordinates of the fitted minimum (length `ndims(mm)`)
+- `Q::Matrix{T}` — symmetric positive-semidefinite curvature matrix of size `(d, d)`
+
+# Examples
+```jldoctest
+julia> using RegisterCore
+
+julia> num = [(i - 1)^2 + (j + 2)^2 for i in -5:5, j in -5:5];
+
+julia> mm = MismatchArray(num, ones(11, 11));
+
+julia> E0, u0, Q = qfit(mm, 0.5);
+
+julia> E0
+0.0
+
+julia> u0
+2-element Vector{Float64}:
+  1.0
+ -2.0
+
+julia> Q ≈ [1.0 0.0; 0.0 1.0]
+true
+```
 """
 function qfit(mm::MismatchArray, thresh::Real; maxsep = size(mm), opt::Bool = true)
     return qfit(mm, thresh, maxsep, opt)
@@ -478,12 +639,40 @@ end
 end
 
 """
-`cs, Qs, mmis = mms2fit!(mms, thresh)` computes the shift and
-quadratic-fit values for the array-of-mismatcharrays `mms`, using a
-threshold of `thresh`. It also prepares `mms` for interpolation,
-modifying the data in-place (after computing `cs` and `Qs`).
+    cs, Qs, mmis = mms2fit!(mms, thresh)
+
+Compute per-aperture shifts and quadratic-fit matrices for the N-dimensional
+array-of-`MismatchArray`s `mms`, using `thresh` as the denominator threshold.
+Also prepares `mms` for interpolation, modifying it in-place after extracting
+`cs` and `Qs`.
 
-The return values are suited for input the `fixed_λ` and `auto_λ`.
+The dimension `N` of the container `mms` must equal the number of spatial
+dimensions of each `MismatchArray` element. For example, a 2×3 matrix of 2D
+mismatch arrays is valid; a `Vector` of 2D mismatch arrays is not.
+
+# Returns
+- `cs`: `Array{SVector{N,T},N}` — per-aperture shift positions
+- `Qs`: `Array{SMatrix{N,N,T},N}` — per-aperture quadratic curvature matrices
+- `mmis`: array of interpolated `MismatchArray`s, suitable for input to
+  `RegisterPenalty.fixed_λ` and `RegisterPenalty.auto_λ`
+
+# Examples
+```jldoctest
+julia> using RegisterCore
+
+julia> num1 = [(i - 1)^2 + (j + 2)^2 for i in -5:5, j in -5:5];
+
+julia> num2 = [(i + 1)^2 + (j - 1)^2 for i in -5:5, j in -5:5];
+
+julia> denom = ones(11, 11);
+
+julia> mms = reshape([MismatchArray(num1, denom), MismatchArray(num2, denom)], 1, 2);
+
+julia> cs, Qs, mmis = mms2fit!(mms, 0.5);
+
+julia> cs[1, 1] ≈ [1.0, -2.0]
+true
+```
 """
 function mms2fit!(mms::AbstractArray{A, N}, thresh) where {A <: MismatchArray, N}
     T = eltype(eltype(A))
diff --git a/test/runtests.jl b/test/runtests.jl
index b2b9c44..a429b47 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,97 +1,169 @@
-import RegisterFit
-using Test, CoordinateTransformations, Interpolations, ImageBase, ImageTransformations, LinearAlgebra
-using RegisterCore
-
-using RegisterUtilities
-
-@testset "qfit" begin
-    denom = ones(11, 11)
-    Q = rand(Float64, 2, 2); Q = Q' * Q
-    num = quadratic(11, 11, [1, -2], Q)
-    E0, cntr, Qf = @inferred(RegisterFit.qfit(MismatchArray(num, denom), 1.0e-3))
-    @test abs(E0) < eps()
-    @test cntr ≈ [1, -2]
-    @test Qf ≈ Q
-
-    num = num .+ 5
-    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), 1.0e-3)
-    @test E0 ≈ 5
-    @test cntr ≈ [1, -2]
-    @test Qf ≈ Q
-
-    num = quadratic(11, 13, [2, -4], Q)
-    thresh = 1.0e-3
-    scale = rand(Float64, size(num)) .+ thresh
-    denom = ones(size(num)) .* scale
-    @test all(denom .> thresh)
-    num = num .* scale
-    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
-    @test abs(E0) < eps()
-    @test cntr ≈ [2, -4]
-    @test Qf ≈ Q
-
-    # Degenerate solutions
-    Q = [1 0; 0 0]
-    denom = ones(13, 11)
-    num = quadratic(13, 11, [2, -4], Q)
-    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
-    @test abs(E0) < eps()
-    @test cntr[1] ≈ 2
-    @test Qf ≈ Q
-    a = rand(2) .+ 0.1
-    Q = a * a'
-    num = quadratic(13, 11, [2, -4], Q)
-    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
-    @test abs(E0) < eps()
-    @test abs(dot(cntr - [2, -4], a)) < eps()
-    @test ≈(Qf, Q, atol = 1.0e-12)
-
-    # Settings with very few above-threshold data points
-    # Just make sure there are no errors
-    denom0 = ones(5, 5)
-    Q = rand(Float64, 2, 2); Q = Q' * Q
-    num0 = quadratic(5, 5, [0, 0], Q)
-    denom = copy(denom0); denom[1:2, 1] *= 100; denom[5, 5] *= 100
-    num = copy(num0); num[1:2, 1] *= 100; num[5, 5] *= 100
-    thresh = 2
-    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
-
-    ### qbuild
-    A = RegisterFit.qbuild(2, [-1, 1], [0.3 0; 0 0.5], (5, 5))
-    v1 = 0.3 * ((-5:5) .+ 1) .^ 2
-    v2 = 0.5 * ((-5:5) .- 1) .^ 2
-    @test A.data ≈ v1 .+ v2' .+ 2
-end
-
-@testset "PAT" begin
-    # Principal Axes Transformation
-    fixed = zeros(10, 11)
-    fixed[2, 3:7] .= 1
-    fixed[3, 2:8] .= 1
-    moving = zeros(10, 11)
-    moving[3:7, 8] .= 1
-    moving[2:8, 7] .= 1
-    fmean, fvar = RegisterFit.principalaxes(fixed)
-    tfm = RegisterFit.pat_rotation((fmean, fvar), moving)
-    for i in 1:2
-        S = tfm[i].linear
-        @test abs(S[1, 2]) ≈ 1
-        @test abs(S[2, 1]) ≈ 1
-        @test abs(S[1, 1]) < 1.0e-8
-        @test abs(S[2, 2]) < 1.0e-8
-    end
-
-    F = meanfinite(abs.(fixed); dims = (1, 2))[1]
-
-    df = zeros(2)
-    movinge = extrapolate(interpolate(moving, BSpline(Linear())), NaN)
-    origin_dest = center(movinge)
-    origin_src = center(fixed)
-    for i in 1:2
-        # mov = TransformedArray(movinge, tfm[i])
-        translation = tfm[i].translation - tfm[i].linear * origin_dest + origin_src
-        tform = AffineMap(tfm[i].linear, translation)
-        df[i] = meanfinite(abs.(fixed - [movinge(tform([idx[1], idx[2]])...) for idx in CartesianIndices(fixed)]); dims = (1, 2))[1]
-    end
-    @test minimum(df) < 1.0e-4 * F
-end
+import RegisterFit
+using Test, Aqua, ExplicitImports, Documenter, CoordinateTransformations, Interpolations, ImageBase, ImageTransformations, LinearAlgebra
+using RegisterCore
+
+using RegisterUtilities
+
+@testset "Doctests" begin
+    DocMeta.setdocmeta!(RegisterFit, :DocTestSetup, :(using RegisterFit, RegisterCore); recursive=true)
+    doctest(RegisterFit; manual=false)
+end
+
+@testset "Aqua" begin
+    Aqua.test_all(RegisterFit)
+end
+
+@testset "ExplicitImports" begin
+    test_explicit_imports(RegisterFit)
+end
+
+@testset "qfit" begin
+    denom = ones(11,11)
+    Q = rand(Float64,2,2); Q = Q'*Q
+    num = quadratic(11, 11, [1,-2], Q)
+    E0, cntr, Qf = @inferred(RegisterFit.qfit(MismatchArray(num, denom), 1e-3))
+    @test abs(E0) < eps()
+    @test cntr ≈ [1,-2]
+    @test Qf ≈ Q
+
+    num = num.+5
+    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), 1e-3)
+    @test E0 ≈ 5
+    @test cntr ≈ [1,-2]
+    @test Qf ≈ Q
+
+    num = quadratic(11, 13, [2,-4], Q)
+    thresh = 1e-3
+    scale = rand(Float64,size(num)) .+ thresh
+    denom = ones(size(num)).*scale
+    @test all(denom .> thresh)
+    num = num.*scale
+    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
+    @test abs(E0) < eps()
+    @test cntr ≈ [2,-4]
+    @test Qf ≈ Q
+
+    # Degenerate solutions
+    Q = [1 0; 0 0]
+    denom = ones(13, 11)
+    num = quadratic(13, 11, [2,-4], Q)
+    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
+    @test abs(E0) < eps()
+    @test cntr[1] ≈ 2
+    @test Qf ≈ Q
+    a = rand(2).+0.1
+    Q = a*a'
+    num = quadratic(13, 11, [2,-4], Q)
+    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
+    @test abs(E0) < eps()
+    @test abs(dot(cntr-[2,-4], a)) < eps()
+    @test ≈(Qf, Q, atol=1e-12)
+
+    # Settings with very few above-threshold data points
+    # Just make sure there are no errors
+    denom0 = ones(5,5)
+    Q = rand(Float64,2,2); Q = Q'*Q
+    num0 = quadratic(5, 5, [0,0], Q)
+    denom = copy(denom0); denom[1:2,1] *= 100; denom[5,5] *= 100
+    num = copy(num0); num[1:2,1] *= 100; num[5,5] *= 100
+    thresh = 2
+    E0, cntr, Qf = RegisterFit.qfit(MismatchArray(num, denom), thresh)
+
+    ### qbuild
+    A = RegisterFit.qbuild(2, [-1,1], [0.3 0; 0 0.5], (5,5))
+    v1 = 0.3*((-5:5).+1).^2
+    v2 = 0.5*((-5:5).-1).^2
+    @test A.data ≈ v1.+v2'.+2
+end
+
+@testset "uisvalid and uclamp!" begin
+    # maxshift must exceed register_half (0.5001) for the bounds to be non-trivial
+    maxshift = (3, 3, 4)
+    u_valid = reshape([1.0, -2.0, 0.5], 3, 1)
+    @test RegisterFit.uisvalid(u_valid, maxshift)
+    u_invalid = reshape([2.6, -2.0, 0.5], 3, 1)
+    @test !RegisterFit.uisvalid(u_invalid, maxshift)
+
+    u_clamp = reshape([5.0, -6.0, 0.1], 3, 1)
+    RegisterFit.uclamp!(u_clamp, maxshift)
+    @test abs(u_clamp[1]) < maxshift[1]
+    @test abs(u_clamp[2]) < maxshift[2]
+    @test abs(u_clamp[3]) < maxshift[3]
+end
+
+@testset "qfit edge cases" begin
+    # All below threshold → zero return
+    num = rand(5, 5)
+    denom = zeros(5, 5)
+    E0, c, Q = RegisterFit.qfit(MismatchArray(num, denom), 1.0)
+    @test E0 == 0
+    @test all(c .== 0)
+    @test all(Q .== 0)
+end
+
+@testset "optimize_per_aperture" begin
+    # 1D grid of 2D mismatch arrays; only first shift component is stored per aperture
+    # argmin_mismatch trims edges, so use shifts within -1:1 for a 5×5 array
+    Q = [1.0 0; 0 1.0]
+    mm1 = MismatchArray(quadratic(5, 5, [1, 0], Q), ones(5, 5))
+    mm2 = MismatchArray(quadratic(5, 5, [-1, 0], Q), ones(5, 5))
+    mms = [mm1, mm2]
+    u = RegisterFit.optimize_per_aperture(mms, 0.5)
+    @test size(u) == (1, 2)
+    @test u[1, 1] ≈ 1
+    @test u[1, 2] ≈ -1
+end
+
+@testset "mms2fit!" begin
+    # Grid dimensionality must match the shift dimensionality of each mismatch array
+    Q = [1.0 0; 0 1.0]
+    mm1 = MismatchArray(quadratic(5, 5, [1, -1], Q), ones(5, 5))
+    mm2 = MismatchArray(quadratic(5, 5, [0,  1], Q), ones(5, 5))
+    mm3 = MismatchArray(quadratic(5, 5, [-1, 0], Q), ones(5, 5))
+    mm4 = MismatchArray(quadratic(5, 5, [1,  0], Q), ones(5, 5))
+    mms = reshape([mm1, mm2, mm3, mm4], 2, 2)
+    cs, Qs, mmis = RegisterFit.mms2fit!(mms, 0.5)
+    @test size(cs) == (2, 2)
+    @test cs[1, 1] ≈ [1.0, -1.0] atol=1e-10
+    @test cs[2, 1] ≈ [0.0,  1.0] atol=1e-10
+end
+
+@testset "PAT" begin
+    # Principal Axes Transformation
+    fixed = zeros(10,11)
+    fixed[2,3:7] .= 1
+    fixed[3,2:8] .= 1
+    moving = zeros(10,11)
+    moving[3:7,8] .= 1
+    moving[2:8,7] .= 1
+    fmean, fvar = RegisterFit.principalaxes(fixed)
+    tfm = RegisterFit.pat_rotation((fmean, fvar), moving)
+    for i = 1:2
+        S = tfm[i].linear
+        @test abs(S[1,2]) ≈ 1
+        @test abs(S[2,1]) ≈ 1
+        @test abs(S[1,1]) < 1e-8
+        @test abs(S[2,2]) < 1e-8
+    end
+
+    # Test the array-input convenience wrapper
+    tfms2 = RegisterFit.pat_rotation(fixed, moving)
+    @test length(tfms2) == length(tfm)
+    for i in eachindex(tfm)
+        @test tfms2[i].linear ≈ tfm[i].linear
+    end
+
+    F = meanfinite(abs.(fixed); dims = (1,2))[1]
+
+    df = zeros(2)
+    movinge = extrapolate(interpolate(moving, BSpline(Linear())), NaN)
+    origin_dest = center(movinge)
+    origin_src = center(fixed)
+    for i = 1:2
+        # mov = TransformedArray(movinge, tfm[i])
+        translation = tfm[i].translation - tfm[i].linear*origin_dest + origin_src
+        tform = AffineMap(tfm[i].linear,translation)
+        df[i] = meanfinite(abs.(fixed-[movinge(tform([idx[1], idx[2]])...) for idx in CartesianIndices(fixed)]); dims = (1,2))[1]
+    end
+    @test minimum(df) < 1e-4*F
+end