牛顿拉普森亚像素搜索：从整像素到1/100像素的精度精化

简介：牛顿拉普森（NR）迭代法结合亚像素搜索的MATLAB实现，面向图像处理、机器视觉及数值计算等领域的研发人员与学习者，解决特征点精确定位、亚像素位移估计及非线性方程求解问题。该压缩包共包含2个m脚本文件，整体体积仅3KB，结构紧凑，便于快速阅读和直接嵌入到已有的图像配准或目标跟踪流程中。代码完整实现了NR迭代的核心步骤，涵盖误差函数构建、雅可比矩阵与海森矩阵计算、线性化迭代公式求解，并通过数字图像相关等示例清晰演示亚像素搜索的收敛过程，同时留有收敛阈值和迭代步长等参数调整入口。目前已吸引374人学习下载，特别适合需要快速掌握NR算法原理、验证亚像素精度方法或开展相关课题研究的开发者参考。借助这份源码，读者既能理解从梯度信息到位移增量的完整推导链路，也能将算法灵活迁移至光学测量、视频分析等具体应用场景。

1. NR 亚像素搜索：和 5G NR 无关，却能决定视觉系统的精度上限

如果你是被“NR”这个词带进来的，先确认一下：这里不是 5G New Radio，也不是 paging、PSS/SSS 那套物理层流程，而是 Newton-Raphson。在机器视觉里，Newton-Raphson 常和亚像素搜索绑定在一起，解决一个很实际的问题：整像素匹配只能给出像素网格上的答案，而产线对位、晶圆定位、PCB 对差需要 1/10 甚至 1/100 像素的重复精度。常见做法是先做一次粗搜索锁定整数像素位置，再以这个位置为起点，用 NR 对连续化的相似度曲面求极值点。适合做视觉检测、图像配准或底层算法库的工程师读。后文会给出可以直接落地的 Python 实现，以及参数怎么设、不收敛时看什么。

2. 牛顿拉普森亚像素搜索的数学模型：把相关峰变成可优化曲面

亚像素搜索的前提是“本来只能算整数值”的相似度函数可以被扩展到连续空间。模板匹配的分数只在每个整数像素采样处存在，NR 又依赖函数的一阶和二阶导数，所以目标函数必须先连续化。另一个容易被忽视的点是：亚像素问题本质上是一个极小值优化问题，而不是直接求相关峰坐标。用 NR 做优化，需要同时计算梯度向量和 Hessian 矩阵，这一步做对了，后面的代码才稳定。

2.1 为什么亚像素搜索必须让目标函数连续

以最常用的平方差和（SSD）作为匹配度量，目标函数可以写成：

E(p)=Σ_r [ I(x+p+r) - T(r) ]²

其中 p=(dx,dy) 是待求的亚像素平移，r 是模板内的像素坐标，T(r) 是模板灰度。整像素搜索时 p 只在整数集合里取值；要让 NR 迭代，p 必须变成实数值，I(x+p+r) 就需要从离散图像中插值得到。插值阶数直接决定曲面光滑程度和结果偏置。常用选择如下：

我一般用三次卷积或 B 样条做基准测试。双线性虽然快，但导数在整数像素上不连续，NR 的梯度会跳变，收敛后容易停在离正确值 0.1 像素左右的位置。目标函数的连续化选型，直接影响后续所有步骤，值得先用合成数据跑一遍再定。

2.2 牛顿拉普森迭代：从求根到求极值

牛顿法最初用于求 f(x)=0 的根，迭代式是 x←x-f(x)/f'(x)。求函数极值时，令梯度 g(p)=∇E(p)=0，也就是把“求梯度为零的点”变成一个求根问题，于是得到多维牛顿-拉普森更新：

p←p - H⁻¹g

其中 g 是梯度向量，H 是 Hessian 矩阵。如果 E 在局部是二次曲面，这一步就能直接到极值；如果不是，也能保证在 Hessian 正定且步长不过大的情况下快速收敛。相比梯度下降，NR 不需要手工调学习率，因为 H⁻¹ 相当于对梯度做了一次包含曲率信息的白化变换。代价是 Hessian 必须可逆且半正定；一旦遇到平坦区域或周期纹理，H 会病态，更新方向就会乱跳。这个问题留到第 4 章解决，先看目标函数本身怎么展开。

2.3 以 SSD 为例的梯度与 Hessian 推导

令残差 r_i=I(x_i+p)-T(x_i)，雅可比向量 J_i=∇I(x_i+p)，即图像在采样点处的灰度梯度。对 SSD 求导，梯度为：

g = Σ_i 2r_i J_i

精确 Hessian 是：

H = Σ_i 2(J_i J_iᵀ + r_i ∇²I)

实际代码里几乎都会忽略残差乘以图像二阶导那一项，因为收敛点附近残差均值接近零，且插值二阶导会带来噪声。于是 Hessian 近似为：

H ≈ Σ_i 2J_i J_iᵀ

这个近似叫 Gauss-Newton，是牛顿-拉普森在最小二乘场景下的标准变体。它不再需要算插值函数的二阶导，稳定性好很多。为了确认自己的梯度解析式没有算错，可以先用有限差分校验。下面这个函数就是通用的数值梯度计算：

import numpy as np def numeric_gradient(func, p, eps=1e-2): g = np.zeros_like(p) for i in range(p.size): dp = np.zeros_like(p) dp[i] = eps g[i] = (func(p + dp) - func(p - dp)) / (2 * eps) return g

中心差分比单侧差分误差小一阶，eps 通常取 1e-3 到 1e-2。注意，如果图像灰度是整数类型，eps 小于 1 时插值输出变化量可能淹没在噪声里，建议先把图像转成 float32。这套模型只用二维平移就能说明 NR 亚像素搜索的骨架，扩展到仿射参数只是把 p 和 H 的维度变大，方法不变。

3. Python 实现牛顿拉普森亚像素搜索：从插值采样到 NR 迭代闭环

直接给一个最小可运行实现。这个实现没有调第三方优化器，目的是让 NR 的每一步都可见。输入是一张目标图像 img、一个从参考图像裁出的模板 T，以及整像素粗匹配中心 (r0,c0)；输出是亚像素偏移 p=(dy,dx)。实际工程里你会在模板匹配之后调用它，把整数峰值精化到小数位置。

3.1 三次样条采样与坐标约定

首先建立坐标网格。模板 T 的每个像素坐标用 (row, col) 表示，row 方向对应偏移 dy，col 方向对应 dx。采样函数如下：

from scipy.ndimage import map_coordinates def sample_at(img, r0, c0, grid_r, grid_c, offset): r = np.clip(r0 + offset[0] + grid_r, 0, img.shape[0] - 1) c = np.clip(c0 + offset[1] + grid_c, 0, img.shape[1] - 1) vals = map_coordinates( img, [r.ravel(), c.ravel()], order=3, mode='nearest' ) return vals.reshape(grid_r.shape)

map_coordinates 的第一个数组对应行坐标，第二个对应列坐标，传入顺序错了会让结果沿对角线错位。order=3 表示使用三次 B 样条，导数连续，适合 NR。mode='nearest' 用来处理越界，但边界附近插值会变得平坦，所以提取模板时最好保留至少两个像素的边界余量。

3.2 数值梯度与 Hessian 计算

完整计算 Hessian 需要二阶差分。对二维 p 来说，H 是一个 2×2 矩阵，需要计算对角线和非对角线四项。实现如下：

def numeric_hessian(func, p, eps=1e-2): n = p.size H = np.zeros((n, n)) for i in range(n): for j in range(i, n): dp_i = np.zeros_like(p) dp_j = np.zeros_like(p) dp_i[i] = eps dp_j[j] = eps if i == j: H[i, j] = ( func(p + dp_i) - 2.0 * func(p) + func(p - dp_i) ) / (eps * eps) else: H[i, j] = ( func(p + dp_i + dp_j) - func(p + dp_i - dp_j) - func(p - dp_i + dp_j) + func(p - dp_i - dp_j) ) / (4.0 * eps * eps) H[j, i] = H[i, j] return H

非对角线项用四差分公式，中心差分的精度比反复套一阶差分更高，而且不需要额外调用过多函数。代价是每次 Hessian 计算要多次求 cost，模板尺寸小的场景可以接受；实时系统请换解析梯度或 Gauss-Newton。

3.3 NR 主循环与收敛判据

下面是主函数。初始 p 设为 0，表示从整像素位置开始修偏。每次迭代求解 HΔ=g，然后执行 p←p-Δ。如果成本没有下降，用回溯步长收住；如果 Hessian 奇异，用最小二乘兜底。

def refine_nr(img, T, r0, c0, max_iter=20, tol=1e-6, eps=1e-2): h, w = T.shape grid_r, grid_c = np.mgrid[0:h, 0:w] def cost(p): sampled = sample_at(img, r0, c0, grid_r, grid_c, p) return 0.5 * np.mean((sampled - T) ** 2) p = np.zeros(2, dtype=float) for it in range(max_iter): g = numeric_gradient(cost, p, eps) H = numeric_hessian(cost, p, eps) try: step = np.linalg.solve(H, g) except np.linalg.LinAlgError: step, _, _, _ = np.linalg.lstsq(H, g, rcond=1e-8) p_new = p - step if cost(p_new) < cost(p): p = p_new else: alpha = 0.5 while alpha > 1e-3 and cost(p - alpha * step) >= cost(p): alpha *= 0.5 p = p - alpha * step if np.max(np.abs(step)) < tol: break if np.max(np.abs(p)) > 1.5: p = np.clip(p, -1.5, 1.5) return p, it + 1

逻辑说明：step 是牛顿步长，正常情况一次性更新；成本不降时说明当前二次模型过度外推，用回溯减半，保证目标函数单调下降。tol 控制的是参数更新量而不是成本值，周期纹理下成本差距很小，参数仍需稳定。clamp 到 ±1.5 像素是防止跳到相邻相关峰，粗匹配误差超过 2 个像素时这个限制会失效，需要先做金字塔。

4. NR 亚像素搜索参数调优：Hessian 修正、步长限制与不收敛排查

NR 参数很少，但每个参数都可能让结果完全跑偏。最核心的是 Hessian 的数值稳定性，其次是 eps 和 clamp 的配合。许多失败案例不是算法不对，而是 Hessian 没有做正则化，导致步长直接越出相关峰。

4.1 关键参数表与调优顺序

调优顺序我先调 eps，再调 clamp，最后才碰 Hessian 正则化。eps 影响所有导数，其他参数只能缓解症状。把 eps 从 1e-2 往小调，观察输出偏移是否稳定，如果连续变化说明插值有足够精度；如果跳来跳去，就要换更高阶插值。

4.2 Levenberg-Marquardt 阻尼：Hessian 不正定时的兜底

最简单的 Hessian 修正是用 LM 阻尼。给 H 加上一个倍数单位矩阵或对角矩阵，让 Hessian 至少半正定。对角阻尼更适合 NR 亚像素搜索，因为 dx、dy 的尺度可能不同。实现如下：

def solve_lm(H, g, lam): diag = np.diag(H) * lam H_reg = H + np.diag(diag) try: return np.linalg.solve(H_reg, g) except np.linalg.LinAlgError: return np.linalg.lstsq(H_reg, g, rcond=1e-8)[0]

在 NR 主循环里，每次迭代先关闭试解，如果成本没有下降，就增大 lam 重新算；如果成本下降，则把 lam 减小。常见的启发式是：

if cost(p - step) < cost(p): lam = max(lam / 3, 1e-9) p = p - step else: lam = min(lam * 3, 1e6)

lam 很小时是纯 NR，lam 很大时接近梯度下降，中间过渡就是 LM 对 NR 的鲁棒化改造。这套机制在参数从 2 个扩到 6 个时尤其重要，因为仿射参数的尺度差异会让 Hessian 条件数变大。

4.3 不收敛诊断：看梯度量级、行列式和成本下降

不收敛通常有三种表现。第一种是 step 持续很大，p 在相邻峰之间跳来跳去，打印 Hessian 行列式会发现它接近零或反复变号，原因是模板落在平坦区域或周期性结构上。第二种是成本下降但 p 停在离真实位置 0.2 像素处，这多半是插值阶数不够或模板边缘边界被 clip。第三种是第一次迭代 cost 反而增大，说明初始点离极值太远，或者图像有异常噪声点。

诊断时不要盲改参数，先打印内部状态：

if it % 5 == 0: det = np.linalg.det(H) print(f"iter={it}, p={p}, |g|={np.linalg.norm(g):.3e}, " f"detH={det:.3e}, cost={cost(p):.6f}")

提示：如果行列式一直是零，说明模板内容太“空”，亚像素结果本质上不可靠。此时应该返回整像素结果，而不是强行让 NR 吐一个数字。

边界问题也很常见。模板太靠近图边时，sample_at 的 clip 会让部分采样点重复边缘像素，导致梯度虚假变小。我一般要求模板边缘离图像边界至少 3 个像素，否则先扩边再搜索。

5. 从平移扩展到仿射：牛顿拉普森亚像素搜索的参数空间扩张

很多实际场景不是纯平移。相机轻微旋转、镜头畸变、被检物体倾斜都会让模板发生仿射变化。NR 的好处是参数只是从 2 个变成 N 个，迭代骨架完全一样。代价是 Hessian 维度和数值差分成本变高，需要更规范地组织采样。

5.1 仿射参数模型与采样矩阵构造

用六个参数描述坐标映射：

r_out = a·r + b·c + c c_out = d·r + e·c + f

其中 (a,b,d,e) 控制缩放旋转，(c,f) 是行列方向的平移。模板坐标要减去中心，这样可以降低 a、b、d、e 和 c、f 之间的相关性。直接使用原始像素坐标时，大偏移和小旋转会互相耦合，Hessian 条件数会差一个量级。采样代码是：

def affine_transform(grid_r, grid_c, params): a, b, c, d, e, f = params r_out = a * grid_r + b * grid_c + c c_out = d * grid_r + e * grid_c + f return r_out, c_out

放到 NR 里，cost 函数变成：先用仿射变换生成采样坐标，再用 map_coordinates 取灰度，后续梯度、阻尼全是通用逻辑。唯一需要额外小心的是参数尺度不同：c、f 是像素单位，a、e 是无量纲缩放量，b、d 还和模板尺寸相乘。数值差分 eps 对所有参数取同一个值会造成病态，建议给平移参数用 1e-2，给缩放参数用 1e-4 或更小。

5.2 六参数牛顿更新的数值雅可比

用数值差分处理六参数 Hessian，单次迭代至少要调用 cost 数十次，适合离线标定或小模板。工程上更稳的做法是引入解析雅可比：对采样位置求偏导时，灰度梯度可以复用插值函数的导数值，再用链式法则把图像梯度和参数映射的偏导乘起来。这里给出一个简化版本：

def refine_affine(img, T, grid_r, grid_c, p0): h, w = T.shape def cost(p): r, c = affine_transform(grid_r, grid_c, p) vals = map_coordinates( img, [r.ravel(), c.ravel()], order=3, mode='nearest' ) return 0.5 * np.mean((vals.reshape(T.shape) - T) ** 2) p = np.array(p0, dtype=float) for _ in range(20): g = numeric_gradient(cost, p, eps=1e-3) H = numeric_hessian(cost, p, eps=1e-3) H_reg = H + np.diag(np.diag(H) * 1e-6) step = np.linalg.solve(H_reg, g) p = p - step if np.max(np.abs(step)) < 1e-7: break return p

参数说明：p0 一般由粗匹配或上一级金字塔给出；eps 这里取 1e-3，是因为缩放参数的取值范围通常很小。如果 p 中既有像素偏移又有缩放量，更好的做法是把工作空间归一化，让每个参数在数值上都是同一量级，再做差分和 H 正则化。否则 LM 的 diag(H)×lam 会偏向量级大的参数，修正失效。

5.3 多尺度金字塔：扩大搜索盆地

NR 是局部优化算法，初始点必须在最优解附近的收敛半径内。纯平移时收敛半径大约 1 到 2 个像素；进入仿射参数后，旋转角度的收敛半径可能只有几度。要处理粗匹配误差大的情况，常见做法是图像金字塔。上层图像缩小后，原本 10 像素的匹配误差变成 2.5 像素，粗定位更容易落在正确的峰附近；然后在每一层做 NR 精化，把结果放大后再传给下一层作为初始值。每层都要重新提取模板，不能只对原图做一次金字塔。这样 NR 亚像素搜索就不是一个孤立函数，而是整个配准流程的最后一环。

6. 验证技巧：用正弦条纹平移图给 NR 亚像素搜索做基准

最后给一个可直接复制的验证方法，用合成图像测精度和稳定性。之所以用正弦条纹而不是纯随机噪声，是因为条纹图案会出现周期性旁峰，能同时检验 NR 是否收敛到主峰；而随机纹理的匹配曲面太平滑，测不出边界问题。

import numpy as np from scipy.ndimage import gaussian_filter, map_coordinates # 生成带纹理的基准图 rng = np.random.default_rng(42) base = rng.normal(size=(256, 256)).astype(np.float32) base = gaussian_filter(base, 1.0) # 从基准图裁模板 r0, c0, half = 120, 150, 16 T = base[r0-half:r0+half+1, c0-half:c0+half+1] # 已知亚像素平移，构造待测图 true_d = np.array([0.37, 0.21]) rr, cc = np.mgrid[0:256, 0:256] rows = rr - true_d[0] cols = cc - true_d[1] I_shift = map_coordinates( base, [rows.ravel(), cols.ravel()], order=3, mode='nearest' ).reshape(base.shape) # 调用 NR 亚像素搜索 est, iters = refine_nr(I_shift, T, r0, c0) print(f"true={true_d}, est={est}, err={est - true_d}, iters={iters}")

验证脚本要循环多组偏移，比如从 0.0 到 0.9 每隔 0.1 取一组，得到误差曲线。判断标准有两个维度：一是误差均值，反映插值是否有系统偏置；二是误差标准差，反映 NR 对初始噪声的敏感性。如果误差在 0.5px 附近出现固定峰值，优先怀疑插值核本身带有相位偏差；如果误差随模板位置变化，则要考虑边界裁剪影响。把这个脚本写成一个 pytest 参数化用例，每次改插值阶数或参数正则化时跑一遍，比打日志直观得多。

这套验证不复杂，但能把插值、梯度和 Hessian 的问题分开暴露。真正上线前还应该把合成测试换成真实图像序列，用同一位置重复采集 100 张，统计输出偏移的重复性，再决定是否把 clamp 放宽到 ±2px。

本文还有配套的精品资源，点击获取