亚像素数值极值检测算法总结

悠悠 2022-05-19 12:37 217阅读 0赞

# 动机 #

在计算机视觉领域，经常需要检测极值位置，比如SIFT关键点检测、模板匹配获得最大响应位置、统计直方图峰值位置、边缘检测等等，有时只需要像素精度就可以，有时则需要亚像素精度。本文尝试总结几种常用的一维离散数据极值检测方法，几个算法主要来自论文《A Comparison of Algorithms for Subpixel Peak Detection》，加上自己的理解和推导。

# 问题定义 #

给定如下离散值，求其极值位置。可知125为观察极值。

\[ 60 , 80 , 100 , 120 , 125 , 105 , 70 , 55 \] \[60, 80, 100, 120, 125, 105, 70, 55\] \[60,80,100,120,125,105,70,55\]

如果这些离散值是从某个分布 f f f中等间距采样获得，其真正的极值位置应位于120和125之间。

下面给出形式化的定义：给定一组离散值，令 x x x为观测到的极值点位置，其值为 f ( x ) f(x) f(x)，其左右相邻位置的值为 f ( x − 1 ) f(x-1) f(x−1)和 f ( x + 1 ) f(x+1) f(x\+1)，真正的极值点位置为 x + δ x+\\delta x\+δ，令 δ ^ \\hat\{\\delta\} δ^为 δ \\delta δ的估计值。

# 算法 #

假设\*\* x x x的邻域\*\*可通过某个模型进行近似，如高斯近似、抛物线近似，则可以利用 x x x的邻域信息根据模型估计出极值。使用的模型不同就有不同的算法，具体如下。

## 高斯近似 ##

一维[高斯函数][Link 1]如下：

y = y m a x ⋅ e x p ( − ( x − μ ) 2 2 σ 2 ) y = y\_\{max\} \\cdot exp(-\\frac\{(x-\\mu)^2\}\{2\\sigma^2\}) y=ymax⋅exp(−2σ2(x−μ)2)  
当 y m a x = 1 2 σ π y\_\{max\}=\\frac\{1\}\{\\sqrt\{2\\sigma\}\\pi\} ymax=2σπ1时为标准高斯函数，形如

![标准高斯函数][izOixs.png]

假设 x x x的邻域可用高斯近似，用 ( x , f ( x ) ) (x, f(x)) (x,f(x))、 ( x − 1 , f ( x − 1 ) ) (x-1, f(x-1)) (x−1,f(x−1))、 ( x + 1 , f ( x + 1 ) ) (x+1, f(x+1)) (x\+1,f(x\+1))三点对高斯函数进行拟合，获得模型参数 μ \\mu μ即为峰值位置， δ ^ = μ − x \\hat\{\\delta\}=\\mu - x δ^=μ−x。将三点带入上面的高斯函数两边同时取对数求得：

δ ^ = 1 2 ln ⁡ ( f ( x − 1 ) ) − ln ⁡ ( f ( x + 1 ) ) ln ⁡ ( f ( x − 1 ) ) − 2 ln ⁡ ( f ( x ) ) + ln ⁡ ( f ( x + 1 ) ) \\hat\{\\delta\} = \\frac\{1\}\{2\} \\frac\{\\ln(f(x-1)) - \\ln(f(x+1))\}\{\\ln(f(x-1)) - 2\\ln(f(x)) + \\ln(f(x+1))\} δ^=21ln(f(x−1))−2ln(f(x))\+ln(f(x\+1))ln(f(x−1))−ln(f(x\+1))

下面可以看到，**高斯近似相当于取对数后的抛物线近似**。

## 抛物线近似 ##

使用抛物线近似 x x x的局部，可以将 ( x , f ( x ) ) (x, f(x)) (x,f(x))、 ( x − 1 , f ( x − 1 ) ) (x-1, f(x-1)) (x−1,f(x−1))、 ( x + 1 , f ( x + 1 ) ) (x+1, f(x+1)) (x\+1,f(x\+1))三点带入 y = a ( x − b ) 2 + c y=a(x-b)^2+c y=a(x−b)2\+c求参数 b b b即为估计的极值位置，也可采用[泰勒展开][Link 2]（[牛顿法][Link 3]）来求极值。泰勒公式实际上是一种利用高阶导数通过多项式近似函数的方法，下面的图示可直观理解这种近似，图示为通过泰勒公式近似原点附近的正弦曲线：

![泰勒近似正弦曲线][izzoIP.png]

泰勒近似 x x x附近，如只取到二阶则为抛物线近似。假设高阶可导，极值为 f ( x + δ ) f(x+\\delta) f(x\+δ)，则根据泰勒公式，

f ( x + δ ) = f ( x ) + f ′ ( x ) δ + 1 2 f ′ ′ ( x ) δ 2 + O ( δ 3 ) f(x+\\delta) = f(x) + f&\#x27;(x)\\delta + \\frac\{1\}\{2\} f&\#x27;&\#x27;(x)\\delta^2 + O(\\delta^3) f(x\+δ)=f(x)\+f′(x)δ\+21f′′(x)δ2\+O(δ3)

极值处导数为0，这里 x x x为常数 δ \\delta δ为变量，两边同时对 δ \\delta δ求导，忽略高阶项可得

f ′ ( x + δ ^ ) = f ′ ( x ) + f ′ ′ ( x ) δ ^ = 0 f&\#x27;(x+\\hat\{\\delta\}) = f&\#x27;(x) + f&\#x27;&\#x27;(x)\\hat\{\\delta\} = 0 f′(x\+δ^)=f′(x)\+f′′(x)δ^=0

使用一阶微分和二阶微分近似 f ′ ( x ) f&\#x27;(x) f′(x)和 f ′ ′ ( x ) f&\#x27;&\#x27;(x) f′′(x)得

δ ^ = − f ′ ( x ) f ′ ′ ( x ) = − ( f ( x + 1 ) − f ( x − 1 ) ) / 2 ( f ( x + 1 ) − f ( x ) ) − ( f ( x ) − f ( x − 1 ) ) = 1 2 f ( x − 1 ) − f ( x + 1 ) f ( x + 1 ) − 2 f ( x ) + f ( x − 1 ) \\hat\{\\delta\} = - \\frac\{f&\#x27;(x)\}\{f&\#x27;&\#x27;(x)\} = - \\frac\{(f(x+1)-f(x-1))/2\}\{(f(x+1)-f(x))-(f(x) - f(x-1))\}= \\frac\{1\}\{2\}\\frac\{f(x-1)-f(x+1)\}\{f(x+1)-2f(x)+ f(x-1)\} δ^=−f′′(x)f′(x)=−(f(x\+1)−f(x))−(f(x)−f(x−1))(f(x\+1)−f(x−1))/2=21f(x\+1)−2f(x)\+f(x−1)f(x−1)−f(x\+1)  
与带入抛物线求参数的结果是一致的，加上对数则与高斯近似一致。

## 质心算法 ##

> In physics, the **center of mass** of a distribution of mass in space is **the unique point where the weighted relative position of the distributed mass sums to zero**, or the point where if a force is applied it moves in the direction of the force without rotating.——[Center of mass wiki][]

![质心][izz7Pf.png]

若将 x x x、 x − 1 x-1 x−1、 x + 1 x+1 x\+1看成质点，将 f ( x ) f(x) f(x)、 f ( x − 1 ) f(x-1) f(x−1)、 f ( x + 1 ) f(x+1) f(x\+1)看成质点的质量，则可以把质心作为极值的估计。根据质点相对质心位置的质量加权和为零，可求得质心位置。令 R R R为质心坐标， m m m和 r r r分别为质点质量和坐标，则 n n n个质点的质心满足

∑ i = 1 n m i ( r i − R ) = 0 \\sum\_\{i=1\}^n m\_i(r\_i - R) = 0 i=1∑nmi(ri−R)=0

令 M = ∑ i = 1 n m i M = \\sum\_\{i=1\}^n m\_i M=∑i=1nmi，质心坐标为

R = 1 M ∑ i = 1 n m i r i R = \\frac\{1\}\{M\} \\sum\_\{i=1\}^n m\_ir\_i R=M1i=1∑nmiri

带入得

x + δ ^ = ( x − 1 ) f ( x − 1 ) + x f ( x ) + ( x + 1 ) f ( x + 1 ) f ( x − 1 ) + f ( x ) + f ( x + 1 ) x + \\hat\{\\delta\} = \\frac\{(x-1)f(x-1)+xf(x)+(x+1)f(x+1)\}\{f(x-1)+f(x)+f(x+1)\} x\+δ^=f(x−1)\+f(x)\+f(x\+1)(x−1)f(x−1)\+xf(x)\+(x\+1)f(x\+1)

δ ^ = f ( x + 1 ) − f ( x − 1 ) f ( x − 1 ) + f ( x ) + f ( x + 1 ) \\hat\{\\delta\} = \\frac\{f(x+1)-f(x-1)\}\{f(x-1)+f(x)+f(x+1)\} δ^=f(x−1)\+f(x)\+f(x\+1)f(x\+1)−f(x−1)

以上考虑的是3质点系统的质心，还可考虑5质点、7质点等，甚至考虑所有点。

## 线性插值 ##

这个模型假设在极值两侧是线性增长和线性下降的，且**上升和下降的速度相同**，即 y = k x + b y=kx+b y=kx\+b，上升侧 k > 0 k>0 k>0，下降侧 k < 0 k<0 k<0，两者绝对值相同，可以利用这个性质求解极值位置。

若 f ( x + 1 ) > f ( x − 1 ) f(x+1)>f(x-1) f(x\+1)>f(x−1)则极值位于 ( x , x + 1 ) (x, x+1) (x,x\+1)之间，可列等式

f ( x ) − f ( x − 1 ) x − ( x − 1 ) = f ( x + δ ) − f ( x ) x + δ − x = f ( x + δ ) − f ( x + 1 ) x + 1 − ( x + δ ) \\frac\{f(x) - f(x-1)\}\{x-(x-1)\} = \\frac\{f(x+\\delta)-f(x)\}\{x+\\delta - x\} = \\frac\{f(x+\\delta)-f(x+1)\}\{x+1-(x+\\delta)\} x−(x−1)f(x)−f(x−1)=x\+δ−xf(x\+δ)−f(x)=x\+1−(x\+δ)f(x\+δ)−f(x\+1)

解得

δ ^ = 1 2 f ( x + 1 ) − f ( x − 1 ) f ( x ) − f ( x − 1 ) \\hat\{\\delta\}=\\frac\{1\}\{2\}\\frac\{f(x+1)-f(x-1)\}\{f(x)-f(x-1)\} δ^=21f(x)−f(x−1)f(x\+1)−f(x−1)

同理，若 f ( x − 1 ) > f ( x + 1 ) f(x-1)>f(x+1) f(x−1)>f(x\+1)求得

δ ^ = 1 2 f ( x + 1 ) − f ( x − 1 ) f ( x ) − f ( x + 1 ) \\hat\{\\delta\}=\\frac\{1\}\{2\}\\frac\{f(x+1)-f(x-1)\}\{f(x)-f(x+1)\} δ^=21f(x)−f(x\+1)f(x\+1)−f(x−1)

## 数值微分滤波 ##

这个方法是利用极值处导数为0的性质，**在微分滤波结果上插值得到导数为0的位置**，因已知极值点在 x x x附近，因此只需在 x x x附近做微分和插值即可。插值时取极值点两侧正负值连线的过零点作为极值点的估计，如下图所示

![Linear  interpolation of the peak position][Linear interpolation of the peak position]

论文[Real-time numerical peak detector][]中定义了4阶和8阶线性滤波器 \[ 1 , 1 , 0 , − 1 , − 1 \] \[1, 1, 0, -1, -1\] \[1,1,0,−1,−1\]和 \[ 1 , 1 , 1 , 1 , 0 , − 1 , − 1 , − 1 , − 1 \] \[1,1,1,1,0,-1,-1,-1,-1\] \[1,1,1,1,0,−1,−1,−1,−1\]，对应的函数形式为

g 4 ( x ) = f ( x − 2 ) + f ( x − 1 ) − f ( x + 1 ) − f ( x + 2 ) g\_4(x)=f(x-2)+f(x-1)-f(x+1)-f(x+2) g4(x)=f(x−2)\+f(x−1)−f(x\+1)−f(x\+2)

g 8 ( x ) = f ( x − 4 ) + f ( x − 3 ) + f ( x − 2 ) + f ( x − 1 ) − f ( x + 1 ) − f ( x + 2 ) − f ( x + 3 ) − f ( x + 4 ) g\_8(x)=f(x-4)+f(x-3)+f(x-2)+f(x-1) \\\\ -f(x+1)-f(x+2)-f(x+3)-f(x+4) g8(x)=f(x−4)\+f(x−3)\+f(x−2)\+f(x−1)−f(x\+1)−f(x\+2)−f(x\+3)−f(x\+4)

2阶形式为 g 2 ( x ) = f ( x − 1 ) − f ( x + 1 ) g\_2(x) = f(x-1) -f(x+1) g2(x)=f(x−1)−f(x\+1)，这些滤波器的表现与数值微分滤波器相似。

当 f ( x + 1 ) > f ( x − 1 ) f(x+1)>f(x-1) f(x\+1)>f(x−1)时，极值点位于 ( x , x + 1 ) (x, x+1) (x,x\+1)之间， g ( x ) < 0 g(x)<0 g(x)<0， g ( x + 1 ) > 0 g(x+1)>0 g(x\+1)>0，极值点位置为 g ( x ) g(x) g(x)与 g ( x + 1 ) g(x+1) g(x\+1)连线的过零点，通过斜率求得

δ ^ = g ( x ) g ( x + 1 ) − g ( x ) \\hat\{\\delta\} = \\frac\{g(x)\}\{g(x+1)-g(x)\} δ^=g(x\+1)−g(x)g(x)

若 f ( x − 1 ) > f ( x + 1 ) f(x-1)>f(x+1) f(x−1)>f(x\+1)，则

δ ^ = g ( x − 1 ) g ( x − 1 ) − g ( x ) − 1 \\hat\{\\delta\} = \\frac\{g(x-1)\}\{g(x-1)-g(x)\} - 1 δ^=g(x−1)−g(x)g(x−1)−1

# 总结 #

\*\*这些数值极值检测方法均是先获取观测极值 x x x及其邻域信息，然后综合邻域信息在各自的模型假设下通过插值估计出极值位置。\*\*若能知道数值来自的真实分布，则直接拟合真实分布然后求极值即可，但往往我们并不知道真实的分布是什么，即使知道真实分布，有时为了快速计算，也会采取插值的方式来估计极值，毕竟偏差可接受效果足够好就可以了。应用时，为了抗噪可对数据先平滑然后求极值，具体采用何种方法可在准确和速度间权衡——所用模型与真实分布越相近自然越准确，如果实在不知道怎么选，就实践对比吧（因为我也不知道），毕竟伟大领袖教导过我们——**实践是检验真理的唯一标准**！

# 参考 #

*  [A Comparison of Algorithms for Subpixel Peak Detection][]
 *  [Real-time numerical peak detector][]

个人博客地址：[亚像素数值极值检测算法总结][Link 4]

[Link 1]: https://wiki2.org/en/Gaussian_function
[izOixs.png]: /images/20220519/1be5c59f6e68480db7577bd84fe6dac0.png
[Link 2]: https://wiki2.org/en/Taylor_series
[Link 3]: https://wiki2.org/en/Newton%27s_method
[izzoIP.png]: /images/20220519/2f9c8adb1eef47119236d6bbc691b93a.png
[Center of mass wiki]: https://wiki2.org/en/Center_of_mass
[izz7Pf.png]: /images/20220519/aee265d91a3c4510a8cfe28eda95856a.png
[Linear interpolation of the peak position]: /images/20220519/e67581687f6c4dcb8ddc7c57dba486c5.png
[Real-time numerical peak detector]: https://dl.acm.org/citation.cfm?id=10782
[A Comparison of Algorithms for Subpixel Peak Detection]: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.380.5886&rep=rep1&type=pdf
[Link 4]: https://blog.shinelee.me/2018/08-04-%E4%BA%9A%E5%83%8F%E7%B4%A0%E6%95%B0%E5%80%BC%E6%9E%81%E5%80%BC%E6%A3%80%E6%B5%8B%E7%AE%97%E6%B3%95%E6%80%BB%E7%BB%93.html