動機
在計算機視覺領域，經常需要檢測極值位置，比如SIFT關鍵點檢測、模板匹配獲得最大響應位置、統計直方圖峰值位置、邊緣檢測等等，有時只需要畫素精度就可以，有時則需要亞畫素精度。本文嘗試總結幾種常用的一維離散資料極值檢測方法，幾個演算法主要來自論文《A Comparison of Algorithms for Subpixel Peak Detection》，加上自己的理解和推導。
問題定義
給定如下離散值，求其極值位置。可知125為觀察極值。
[[60, 80, 100, 120, 125, 105, 70, 55]]
如果這些離散值是從某個分佈(f)中等間距取樣獲得，其真正的極值位置應位於120和125之間。
下面給出形式化的定義：給定一組離散值，令(x)為觀測到的極值點位置，其值為(f(x))，其左右相鄰位置的值為(f(x-1))和(f(x+1))，真正的極值點位置為(x+delta)，令(hat{delta})為(delta)的估計值。
演算法
假設(x)的鄰域可通過某個模型進行近似，如高斯近似、拋物線近似，則可以利用(x)的鄰域資訊根據模型估計出極值。使用的模型不同就有不同的演算法，具體如下。
高斯近似
一維高斯函式如下：
[y = y_{max} cdot exp(-frac{(x-mu)^2}{2sigma^2})]
當(y_{max}=frac{1}{sqrt{2sigma}pi})時為標準高斯函式，形如
假設(x)的鄰域可用高斯近似，用((x, f(x)))、((x-1, f(x-1)))、((x+1, f(x+1)))三點對高斯函式進行擬合，獲得模型引數(mu)即為峰值位置，(hat{delta}=mu – x)。將三點帶入上面的高斯函式兩邊同時取對數求得：
[hat{delta} = frac{1}{2} frac{ln(f(x-1)) – ln(f(x+1))}{ln(f(x-1)) – 2ln(f(x)) + ln(f(x+1))}]
下面可以看到，高斯近似相當於取對數後的拋物線近似。
拋物線近似
使用拋物線近似(x)的區域性，可以將((x, f(x)))、((x-1, f(x-1)))、((x+1, f(x+1)))三點帶入(y=a(x-b)^2+c)求引數(b)即為估計的極值位置，也可採用泰勒展開（牛頓法）來求極值。泰勒公式實際上是一種利用高階導數通過多項式近似函式的方法，下面的圖示可直觀理解這種近似，圖示為通過泰勒公式近似原點附近的正弦曲線：
泰勒近似(x)附近，如只取到二階則為拋物線近似。假設高階可導，極值為(f(x+delta))，則根據泰勒公式，
[f(x+delta) = f(x) + f`(x)delta + frac{1}{2} f“(x)delta^2 + O(delta^3)]
極值處導數為0，這裡(x)為常數(delta)為變數，兩邊同時對(delta)求導，忽略高階項可得
[f`(x+hat{delta}) = f`(x) + f“(x)hat{delta} = 0]
使用一階微分和二階微分近似(f`(x))和(f“(x))得
[hat{delta} = – frac{f`(x)}{f“(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)}]
與帶入拋物線求引數的結果是一致的，加上對數則與高斯近似一致。
質心演算法
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
若將(x)、(x-1)、(x+1)看成質點，將(f(x))、(f(x-1))、(f(x+1))看成質點的質量，則可以把質心作為極值的估計。根據質點相對質心位置的質量加權和為零，可求得質心位置。令(R)為質心座標，(m)和(r)分別為質點質量和座標，則(n)個質點的質心滿足
[sum_{i=1}^n m_i(r_i – R) = 0]
令(M = sum_{i=1}^n m_i)，質心座標為
[R = frac{1}{M} sum_{i=1}^n m_ir_i]
帶入得
[x + hat{delta} = frac{(x-1)f(x-1)+xf(x)+(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)}]
以上考慮的是3質點系統的質心，還可考慮5質點、7質點等，甚至考慮所有點。
線性插值
這個模型假設在極值兩側是線性增長和線性下降的，且上升和下降的速度相同，即(y=kx+b)，上升側(k>0)，下降側(k<0)，兩者絕對值相同，可以利用這個性質求解極值位置。
若(f(x+1)>f(x-1))則極值位於((x, x+1))之間，可列等式
[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)}]
解得
[hat{delta}=frac{1}{2}frac{f(x+1)-f(x-1)}{f(x)-f(x-1)}]
同理，若(f(x-1)>f(x+1))求得
[hat{delta}=frac{1}{2}frac{f(x+1)-f(x-1)}{f(x)-f(x+1)}]
數值微分濾波
這個方法是利用極值處導數為0的性質，在微分濾波結果上插值得到導數為0的位置，因已知極值點在(x)附近，因此只需在(x)附近做微分和插值即可。插值時取極值點兩側正負值連線的過零點作為極值點的估計，如下圖所示
論文Real-time numerical peak detector中定義了4階和8階線性濾波器([1, 1, 0, -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_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)]
2階形式為(g_2(x) = f(x-1) -f(x+1))，這些濾波器的表現與數值微分濾波器相似。
當(f(x+1)>f(x-1))時，極值點位於((x, x+1))之間，(g(x)<0)，(g(x+1)>0)，極值點位置為(g(x))與(g(x+1))連線的過零點，通過斜率求得
[hat{delta} = frac{g(x)}{g(x+1)-g(x)}]
若(f(x-1)>f(x+1))，則
[hat{delta} = frac{g(x-1)}{g(x-1)-g(x)} – 1]
總結
這些數值極值檢測方法均是先獲取觀測極值(x)及其鄰域資訊，然後綜合鄰域資訊在各自的模型假設下通過插值估計出極值位置。若能知道數值來自的真實分佈，則直接擬合真實分佈然後求極值即可，但往往我們並不知道真實的分佈是什麼，即使知道真實分佈，有時為了快速計算，也會採取插值的方式來估計極值，畢竟偏差可接受效果足夠好就可以了。應用時，為了抗噪可對資料先平滑然後求極值，具體採用何種方法可在準確和速度間權衡——所用模型與真實分佈越相近自然越準確，如果實在不知道怎麼選，就實踐對比吧（因為我也不知道），畢竟偉大領袖教導過我們——實踐是檢驗真理的唯一標準！
參考
個人部落格地址：亞畫素數值極值檢測演算法總結