动机
在计算机视觉领域，经常需要检测极值位置，比如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]
如果这些离散值是从某个分布ff中等间距采样获得，其真正的极值位置应位于120和125之间。
下面给出形式化的定义：给定一组离散值，令xx为观测到的极值点位置，其值为f(x)f(x)，其左右相邻位置的值为f(x−1)f(x−1)和f(x+1)f(x+1)，真正的极值点位置为x+δx+δ，令δ^δ^为δδ的估计值。
算法
假设xx的邻域可通过某个模型进行近似，如高斯近似、抛物线近似，则可以利用xx的邻域信息根据模型估计出极值。使用的模型不同就有不同的算法，具体如下。
高斯近似
一维高斯函数如下：
y=ymax⋅exp(−(x−μ)22σ2)y=ymax⋅exp(−(x−μ)22σ2)
当ymax=12σ√πymax=12σπ时为标准高斯函数，形如
假设xx的邻域可用高斯近似，用(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δ^=μ−x。将三点带入上面的高斯函数两边同时取对数求得：
δ^=12ln(f(x−1))−ln(f(x+1))ln(f(x−1))−2ln(f(x))+ln(f(x+1))δ^=12ln(f(x−1))−ln(f(x+1))ln(f(x−1))−2ln(f(x))+ln(f(x+1))
下面可以看到，高斯近似相当于取对数后的抛物线近似。
抛物线近似
使用抛物线近似xx的局部，可以将(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))三点带入y=a(x−b)2+cy=a(x−b)2+c求参数bb即为估计的极值位置，也可采用泰勒展开（牛顿法）来求极值。泰勒公式实际上是一种利用高阶导数通过多项式近似函数的方法，下面的图示可直观理解这种近似，图示为通过泰勒公式近似原点附近的正弦曲线：
正在上传…重新上传取消
泰勒近似xx附近，如只取到二阶则为抛物线近似。假设高阶可导，极值为f(x+δ)f(x+δ)，则根据泰勒公式，
f(x+δ)=f(x)+f′(x)δ+12f′′(x)δ2+O(δ3)f(x+δ)=f(x)+f′(x)δ+12f″(x)δ2+O(δ3)
极值处导数为0，这里xx为常数δδ为变量，两边同时对δδ求导，忽略高阶项可得
f′(x+δ^)=f′(x)+f′′(x)δ^=0f′(x+δ^)=f′(x)+f″(x)δ^=0
使用一阶微分和二阶微分近似f′(x)f′(x)和f′′(x)f″(x)得
δ^=−f′(x)f′′(x)=−(f(x+1)−f(x−1))/2(f(x+1)−f(x))−(f(x)−f(x−1))=12f(x−1)−f(x+1)f(x+1)−2f(x)+f(x−1)δ^=−f′(x)f″(x)=−(f(x+1)−f(x−1))/2(f(x+1)−f(x))−(f(x)−f(x−1))=12f(x−1)−f(x+1)f(x+1)−2f(x)+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
若将xx、x−1x−1、x+1x+1看成质点，将f(x)f(x)、f(x−1)f(x−1)、f(x+1)f(x+1)看成质点的质量，则可以把质心作为极值的估计。根据质点相对质心位置的质量加权和为零，可求得质心位置。令RR为质心坐标，mm和rr分别为质点质量和坐标，则nn个质点的质心满足
∑i=1nmi(ri−R)=0∑i=1nmi(ri−R)=0
令M=∑ni=1miM=∑i=1nmi，质心坐标为
R=1M∑i=1nmiriR=1M∑i=1nmiri
带入得
x+δ^=(x−1)f(x−1)+xf(x)+(x+1)f(x+1)f(x−1)+f(x)+f(x+1)x+δ^=(x−1)f(x−1)+xf(x)+(x+1)f(x+1)f(x−1)+f(x)+f(x+1)
δ^=f(x+1)−f(x−1)f(x−1)+f(x)+f(x+1)δ^=f(x+1)−f(x−1)f(x−1)+f(x)+f(x+1)
以上考虑的是3质点系统的质心，还可考虑5质点、7质点等，甚至考虑所有点。
线性插值
这个模型假设在极值两侧是线性增长和线性下降的，且上升和下降的速度相同，即y=kx+by=kx+b，上升侧k>0k>0，下降侧k<0k<0，两者绝对值相同，可以利用这个性质求解极值位置。
若f(x+1)>f(x−1)f(x+1)>f(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+δ)f(x)−f(x−1)x−(x−1)=f(x+δ)−f(x)x+δ−x=f(x+δ)−f(x+1)x+1−(x+δ)
解得
δ^=12f(x+1)−f(x−1)f(x)−f(x−1)δ^=12f(x+1)−f(x−1)f(x)−f(x−1)
同理，若f(x−1)>f(x+1)f(x−1)>f(x+1)求得
δ^=12f(x+1)−f(x−1)f(x)−f(x+1)δ^=12f(x+1)−f(x−1)f(x)−f(x+1)
数值微分滤波
这个方法是利用极值处导数为0的性质，在微分滤波结果上插值得到导数为0的位置，因已知极值点在xx附近，因此只需在xx附近做微分和插值即可。插值时取极值点两侧正负值连线的过零点作为极值点的估计，如下图所示
正在上传…重新上传取消
论文Real-time numerical peak detector中定义了4阶和8阶线性滤波器[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]，对应的函数形式为
g4(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)
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)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阶形式为g2(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)时，极值点位于(x,x+1)(x,x+1)之间，g(x)<0g(x)<0，g(x+1)>0g(x+1)>0，极值点位置为g(x)g(x)与g(x+1)g(x+1)连线的过零点，通过斜率求得
δ^=g(x)g(x+1)−g(x)δ^=g(x)g(x+1)−g(x)
若f(x−1)>f(x+1)f(x−1)>f(x+1)，则
δ^=g(x−1)g(x−1)−g(x)−1δ^=g(x−1)g(x−1)−g(x)−1
总结
这些数值极值检测方法均是先获取观测极值xx及其邻域信息，然后综合邻域信息在各自的模型假设下通过插值估计出极值位置。若能知道数值来自的真实分布，则直接拟合真实分布然后求极值即可，但往往我们并不知道真实的分布是什么，即使知道真实分布，有时为了快速计算，也会采取插值的方式来估计极值，毕竟偏差可接受效果足够好就可以了。应用时，为了抗噪可对数据先平滑然后求极值，具体采用何种方法可在准确和速度间权衡——所用模型与真实分布越相近自然越准确，如果实在不知道怎么选，就实践对比吧（因为我也不知道），毕竟伟大领袖教导过我们——实践是检验真理的唯一标准！