A summary of sub-pixel numerical extremum detection algorithm

Motivation

In the field of computer vision, it is often necessary to detect extreme position, such as sift key detection, template matching to obtain the maximum response position, statistical histogram peak position, edge detection and so on, sometimes only need pixel precision, sometimes sub-pixel accuracy is required. This paper attempts to summarize several common one-dimensional discrete data extreme detection methods, several algorithms mainly from the paper "A Comparison of algorithms for Subpixel Peak Detection", coupled with their own understanding and derivation.

Problem definition

Given the following discrete values, the extremum position is obtained. 125 is to observe the extremum.

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

If these discrete values are obtained from a distribution \ (f\) medium Spacing sample, their true extremum position should be between 120 and 125.

The following is a formal definition: given a set of discrete values, so \ (x\) is the observed extremum point position, its value is \ (f (x) \), the value of its left and right adjacent position is \ (f (x-1) \) and \ (f (x+1) \) , the true extreme point position is \ (x+\delta\), which makes \ (\hat{\delta}\) The estimated value of \ (\delta\) .

Algorithm

Assuming that the neighborhood of \ (x\) can be approximated by a model, such as Gaussian approximation, parabolic approximation, the extremum can be estimated from the model by using the neighborhood information of \ (x\) . There are different algorithms for using different models, as follows.

Gaussian approximation

A Gaussian function is as follows:

\[y = Y_{max} \cdot exp (-\frac{(X-\MU) ^2}{2\sigma^2}) \]
When \ (y_{max}=\frac{1}{\sqrt{2\sigma}\pi}\) is a standard Gaussian function, the shape of the

Suppose that the neighborhood of \ (x\) can be Gaussian approximated, with \ ((x, f (x)) \),\ ((X-1, F (x-1)) \),\ ((X+1, F (x+1)) \) The Gaussian function is fitted at three points to obtain the model parameter \ (\mu\) as the peak position,\ (\hat{\delta}=\mu-x\). The three points are taken into the above Gaussian function to obtain the logarithm at the same time:

\[\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))}\]

As you can see, the Gaussian approximation is equivalent to the parabolic approximation after the logarithm is taken .

Parabolic approximation

With a parabola approximation \ (x\) Local, you can bring \ ((x, f (x)) \),\ ((X-1, F (x-1)), \ ((x+1, F (x+1)) \) three points into \ (y=a (x-b) ^2+c\) parameter \ (b\) is the estimated extremum position, also can use Taylor expansion (Newton method) to find the extremum. The Taylor formula is actually a method that uses the higher-order derivative to pass the polynomial approximation function, and the following illustration intuitively understands this approximation, as shown by the sine curve near the origin by the Taylor formula:

Near the Taylor approximation \ (x\) , a parabolic approximation is obtained if only the second order is taken. Assuming higher order, the extremum is \ (f (x+\delta) \), then according to the Taylor formula,

\[f (X+\delta) = f (x) + F ' (x) \delta + \frac{1}{2} f ' (x) \delta^2 + O (\delta^3) \]

The derivative of the extremum is 0, here \ (x\) is a constant \ (\delta\) as a variable, both sides simultaneously on \ (\delta\) derivative, ignoring high-order items can be

\[f ' (X+\hat{\delta}) = F ' (x) + F ' (x) \hat{\delta} = 0\]

Use the first-order differential and second differential approximation \ (f ' (x) \) and \ (f ' (x) \) to

\[\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)}\]
is consistent with the result of the parameter being taken into the parabola, and the logarithm is consistent with the Gaussian approximation.

centroid algorithm

In physics, the center of mass of a distribution of mass in space are the unique point where the weighted RelA tive position of the distributed mass sums to zero, or the-the point where if a-is-applied it moves in the Directio N of the force without rotating.--center of mass wiki

If \ (x\),\ (x-1\),\ (x+1\) as a particle, will \ (f (x) \),\ (f (x-1) \),\ (f (x+1) \) as the mass of a particle, the centroid can be estimated as the extremum. The centroid position can be obtained according to the mass weighting and zero of the relative centroid position of the particle. ( r\) is the centroid coordinate,\ (m\) and \ (r\) are mass and coordinates respectively, then the centroid of \ (n\) particle is satisfied

\[\sum_{i=1}^n m_i (r_i-r) = 0\]

(M = \sum_{i=1}^n m_i\), the centroid coordinate is

\[r = \frac{1}{m} \sum_{i=1}^n m_ir_i\]

Brought in.

\[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)}\]

The above considerations are the centroid of the 3 particle system, and can also consider 5 particles, 7 particles, and even consider all points.

Linear interpolation

This model assumes linear growth and linear descent on both sides of the extremum, and the same rate of ascent and descent , i.e. \ (y=kx+b\), ascending side \ (k>0\), descent side \ (k<0\) , the absolute value is the same and can be used to solve the extremum position.

If \ (f (x+1) >f (x-1) \) The extremum is between \ ((x, x+1) \) , the equation can be listed

\[\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)}\]

Solution to

\[\hat{\delta}=\frac{1}{2}\frac{f (X+1)-F (x-1)}{f (x)-F (x-1)}\]

Similarly, if \ (f (x-1) >f (x+1) \) obtained

\[\hat{\delta}=\frac{1}{2}\frac{f (X+1)-F (x-1)}{f (x)-F (x+1)}\]

Numerical differential filtering

The method is to use the value of the derivative of the property of 0, in the differential filtering results are worthy of the derivative of the position of 0 , because the known extremum point in the vicinity of \ (x\) , so only in the vicinity of \ (x\) differential and interpolation can be. The value of the extremum point is estimated by the 0 points of the positive and negative lines of the extreme points, as shown in the interpolation.

The paper real-time numerical peak detector with 4-order and 8-order linear Filters ([1, 1, 0,-1, -1]\) and \ ([1,1,1,1,0,-1,-1,-1,-1]\), The corresponding function form is

\[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) \]

The 2-order form is \ (g_2 (x) = f (x-1)-f (x+1) \), and the performance of these filters is similar to that of a numerical differential filter.

When \ (f (x+1) >f (x-1) \) , the extremum point is between \ ((x, x+1) \) , \ (g (x) <0\),\ (g (x+1) >0\), Extreme point position is \ (g (x) \) and \ (g (x+1) \) line over 0 points, the slope is calculated

\[\hat{\delta} = \frac{g (x)}{g (x+1)-g (x)}\]

If \ (f (x-1) >f (x+1) \), then

\[\hat{\delta} = \frac{g (x-1)}{g (x-1)-G (x)}-1\]

Summarize

These numerical extremum detection methods are to obtain the observed Extremum (x\) and its neighborhood information, and then synthesize the neighborhood information to estimate the extremum position by interpolation in the respective model hypothesis. If you can know the real distribution of the value from, then directly fit the real distribution and then the extremum can be, but often we do not know what the real distribution is, even if we know the real distribution, sometimes in order to quickly calculate, will also take interpolation method to estimate the extremum, after all, the deviation acceptable effect is good enough. When applied, in order to anti-noise to the data first smooth and then to find the extremum, the specific method can be measured between accuracy and speed-the more closely used models and the real distribution of natural more accurate, if you do not know how to choose, the practice of comparison (because I do not know), after all, the great leader taught us- Practice is the only standard for testing truth !

Reference

• A Comparison of algorithms for Subpixel Peak Detection
• Real-time Numerical Peak detector

Personal Blog Address: Sub-pixel numerical extremum detection algorithm Summary

A summary of sub-pixel numerical extremum detection algorithm