亞像素數值極值檢測算法總結

動機

在計算機視覺領域，經常需要檢測極值位置，比如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}{2\sigma^2}) \]

當\(y_{max}=\frac{1}{\sqrt{2\sigma}\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)) - 2\ln(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\)及其鄰域信息，然后綜合鄰域信息在各自的模型假設下通過插值估計出極值位置。若能知道數值來自的真實分布，則直接擬合真實分布然后求極值即可，但往往我們並不知道真實的分布是什么，即使知道真實分布，有時為了快速計算，也會采取插值的方式來估計極值，畢竟偏差可接受效果足夠好就可以了。應用時，為了抗噪可對數據先平滑然后求極值，具體采用何種方法可在准確和速度間權衡——所用模型與真實分布越相近自然越准確，如果實在不知道怎么選，就實踐對比吧（因為我也不知道），畢竟偉大領袖教導過我們——實踐是檢驗真理的唯一標准！

參考

• A Comparison of Algorithms for Subpixel Peak Detection
• Real-time numerical peak detector

個人博客地址：亞像素數值極值檢測算法總結