Harris corner measure

The Harris corner detector algorithm relies on a central principle: at a corner, the image intensity will change largely in multiple directions. This can alternatively be formulated by examining the changes of intensity due to shifts in a local window. Around a corner point, the image intensity will change greatly when the window is shifted in an arbitrary direction. Following this intuition and through a clever decomposition, the Harris detector uses the second moment matrix as the basis of its corner decisions. (See corner detection for a more complete derivation). The matrix ${\displaystyle A}$, has also been called the autocorrelation matrix and has values closely related to the derivatives of image intensity.

${\displaystyle A(\mathbf {x} )=\sum _{p,q}w(p,q){\begin{bmatrix}I_{x}^{2}(p,q)&I_{x}I_{y}(p,q)\\I_{x}I_{y}(p,q)&I_{y}^{2}(p,q)\\\end{bmatrix}}}$

where ${\displaystyle I_{x}}$ and ${\displaystyle I_{y}}$ are the respective derivatives (of pixel intensity) in the ${\displaystyle x}$ and ${\displaystyle y}$ direction at point (${\displaystyle p}$,${\displaystyle q}$); ${\displaystyle p}$ and ${\displaystyle q}$ are the position parameters of the weighting function w. The off-diagonal entries are the product of ${\displaystyle I_{x}}$ and ${\displaystyle I_{y}}$, while the diagonal entries are squares of the respective derivatives. The weighting function ${\displaystyle w(x,y)}$ can be uniform, but is more typically an isotropic, circular Gaussian,

${\displaystyle w(x,y)=g(x,y,\sigma )={\frac {1}{2\pi \sigma ^{2}}}e^{\left(-{\frac {x^{2}+y^{2}}{2\sigma ^{2}}}\right)}}$

that acts to average in a local region while weighting those values near the center more heavily.

As it turns out, this ${\displaystyle A}$ matrix describes the shape of the autocorrelation measure as due to shifts in window location. Thus, if we let ${\displaystyle \lambda _{1}}$ and ${\displaystyle \lambda _{2}}$ be the eigenvalues of ${\displaystyle A}$, then these values will provide a quantitative description of how the autocorrelation measure changes in space: its principal curvatures. As Harris and Stephens (1988) point out, the ${\displaystyle A}$ matrix centered on corner points will have two large, positive eigenvalues.^[8] Rather than extracting these eigenvalues using methods like singular value decomposition, the Harris measure based on the trace and determinant is used:

${\displaystyle R=\det(A)-\alpha \operatorname {trace} ^{2}(A)=\lambda _{1}\lambda _{2}-\alpha (\lambda _{1}+\lambda _{2})^{2}}$

where ${\displaystyle \alpha }$ is a constant. Corner points have large, positive eigenvalues and would thus have a large Harris measure. Thus, corner points are identified as local maxima of the Harris measure that are above a specified threshold.

${\displaystyle {\begin{aligned}\{x_{c}\}={\big \{}x_{c}\mid R(x_{c})>R(x_{i}),\forall x_{i}\in W(x_{c}){\big \}},\\R(x_{c})>t_{\text{threshold}}\end{aligned}}}$

where ${\displaystyle \{x_{c}\}}$ are the set of all corner points, ${\displaystyle R(x)}$ is the Harris measure calculated at ${\displaystyle x}$, ${\displaystyle W(x_{c})}$ is an 8-neighbor set centered on ${\displaystyle x_{c}}$ and ${\displaystyle t_{\text{threshold}}}$ is a specified threshold.

8-Point neighborhood

Gaussian scale-space

A Gaussian scale space representation of an image is the set of images that result from convolving a Gaussian kernel of various sizes with the original image. In general, the representation can be formulated as:

${\displaystyle L(\mathbf {x} ,s)=G(s)\otimes I(\mathbf {x} )}$

where ${\displaystyle G(s)}$ is an isotropic, circular Gaussian kernel as defined above. The convolution with a Gaussian kernel smooths the image using a window the size of the kernel. A larger scale, ${\displaystyle s}$, corresponds to a smoother resultant image. Mikolajczyk and Schmid (2001) point out that derivatives and other measurements must be normalized across scales.^[9] A derivative of order ${\displaystyle m}$, ${\displaystyle D_{i_{1},...i_{m}}}$, must be normalized by a factor ${\displaystyle s^{m}}$ in the following manner:

${\displaystyle D_{i_{1},\dots ,i_{m}}(\mathbf {x} ,s)=s^{m}L_{i_{1},\dots ,i_{m}}(\mathbf {x} ,s)}$

These derivatives, or any arbitrary measure, can be adapted to a scale space representation by calculating this measure using a set of scales recursively where the ${\displaystyle n}$th scale is ${\displaystyle s_{n}=k^{n}s_{0}}$. See scale space for a more complete description.

Combining Harris detector across Gaussian scale-space

The Harris-Laplace detector combines the traditional 2D Harris corner detector with the idea of a Gaussian scale space representation in order to create a scale-invariant detector. Harris-corner points are good starting points because they have been shown to have good rotational and illumination invariance in addition to identifying the interesting points of the image.^[10] However, the points are not scale invariant and thus the second-moment matrix must be modified to reflect a scale-invariant property. Let us denote, ${\displaystyle M=\mu (\mathbf {x} ,\sigma _{\mathit {I}},\sigma _{\mathit {D}})}$ as the scale adapted second-moment matrix used in the Harris-Laplace detector.

${\displaystyle M=\mu (\mathbf {x} ,\sigma _{\mathit {I}},\sigma _{\mathit {D}})=\sigma _{D}^{2}g(\sigma _{I})\otimes {\begin{bmatrix}L_{x}^{2}(\mathbf {x} ,\sigma _{D})&L_{x}L_{y}(\mathbf {x} ,\sigma _{D})\\L_{x}L_{y}(\mathbf {x} ,\sigma _{D})&L_{y}^{2}(\mathbf {x} ,\sigma _{D})\end{bmatrix}}}$^[11]

where ${\displaystyle g(\sigma _{I})}$ is the Gaussian kernel of scale ${\displaystyle \sigma _{I}}$ and ${\displaystyle \mathbf {x} =(x,y)}$. Similar to the Gaussian-scale space, ${\displaystyle L(\mathbf {x} )}$ is the Gaussian-smoothed image. The ${\displaystyle \mathbf {\otimes } }$ operator denotes convolution. ${\displaystyle L_{x}(\mathbf {x} ,\sigma _{D})}$ and ${\displaystyle L_{y}(\mathbf {x} ,\sigma _{D})}$ are the derivatives in their respective direction applied to the smoothed image and calculated using a Gaussian kernel with scale ${\displaystyle \sigma _{D}}$. In terms of our Gaussian scale-space framework, the ${\displaystyle \sigma _{I}}$ parameter determines the current scale at which the Harris corner points are detected.

Building upon this scale-adapted second-moment matrix, the Harris-Laplace detector is a twofold process: applying the Harris corner detector at multiple scales and automatically choosing the characteristic scale.

Multi-scale Harris corner points

The algorithm searches over a fixed number of predefined scales. This set of scales is defined as:

${\displaystyle {\sigma _{1}\dots \sigma _{n}}={k^{1}\sigma _{0}\dots k^{n}\sigma _{0}}}$

Mikolajczyk and Schmid (2004) use ${\displaystyle k=1.4}$. For each integration scale, ${\displaystyle \sigma _{I}}$, chosen from this set, the appropriate differentiation scale is chosen to be a constant factor of the integration scale: ${\displaystyle \sigma _{D}=s\sigma _{I}}$. Mikolajczyk and Schmid (2004) used ${\displaystyle s=0.7}$.^[11] Using these scales, the interest points are detected using a Harris measure on the ${\displaystyle \mu (\mathbf {x} ,\sigma _{\mathit {I}},\sigma _{\mathit {D}})}$ matrix. The cornerness, like the typical Harris measure, is defined as:

${\displaystyle {\mathit {cornerness}}=\det(\mu (\mathbf {x} ,\sigma _{\mathit {I}},\sigma _{\mathit {D}}))-\alpha \operatorname {trace} ^{2}(\mu (\mathbf {x} ,\sigma _{\mathit {I}},\sigma _{\mathit {D}}))}$

Like the traditional Harris detector, corner points are those local (8 point neighborhood) maxima of the cornerness that are above a specified threshold.

Characteristic scale identification

An iterative algorithm based on Lindeberg (1998) both spatially localizes the corner points and selects the characteristic scale.^[6] The iterative search has three key steps, that are carried for each point ${\displaystyle \mathbf {x} }$ that were initially detected at scale ${\displaystyle \sigma _{I}}$ by the multi-scale Harris detector (${\displaystyle k}$ indicates the ${\displaystyle kth}$ iteration):

• Choose the scale ${\displaystyle \sigma _{I}^{(k+1)}}$ that maximizes the Laplacian-of-Gaussians (LoG) over a predefined range of neighboring scales. The neighboring scales are typically chosen from a range that is within a two scale-space neighborhood. That is, if the original points were detected using a scaling factor of ${\displaystyle 1.4}$ between successive scales, a two scale-space neighborhood is the range ${\displaystyle t\in [0.7,\dots ,1.4]}$. Thus the Gaussian scales examined are: ${\displaystyle \sigma _{I}^{(k+1)}=t\sigma _{I}^{k}}$. The LoG measurement is defined as:

${\displaystyle |\operatorname {LoG} (\mathbf {x} ,\sigma _{I})|=\sigma _{I}^{2}\left|L_{xx}(\mathbf {x} ,\sigma _{I})+L_{yy}(\mathbf {x} ,\sigma _{I})\right|}$

where ${\displaystyle L_{xx}}$ and ${\displaystyle L_{yy}}$ are the second derivatives in their respective directions.^[12] The ${\displaystyle \sigma _{I}^{2}}$ factor (as discussed above in Gaussian scale-space) is used to normalize the LoG across scales and make these measures comparable, thus making a maximum relevant. Mikolajczyk and Schmid (2001) demonstrate that the LoG measure attains the highest percentage of correctly detected corner points in comparison to other scale-selection measures.^[9] The scale which maximizes this LoG measure in the two scale-space neighborhood is deemed the characteristic scale, ${\displaystyle \sigma _{I}^{(k+1)}}$, and used in subsequent iterations. If no extrema, or maxima of the LoG is found, this point is discarded from future searches.

• Using the characteristic scale, the points are spatially localized. That is to say, the point ${\displaystyle \mathbf {x} ^{(k+1)}}$ is chosen such that it maximizes the Harris corner measure (cornerness as defined above) within an 8×8 local neighborhood.
• Stopping criterion: ${\displaystyle \sigma _{I}^{(k+1)}==\sigma _{I}^{(k)}}$ and ${\displaystyle \mathbf {x} ^{(k+1)}==\mathbf {x} ^{(k)}}$.

If the stopping criterion is not met, then the algorithm repeats from step 1 using the new ${\displaystyle k+1}$ points and scale. When the stopping criterion is met, the found points represent those that maximize the LoG across scales (scale selection) and maximize the Harris corner measure in a local neighborhood (spatial selection).

Mathematical theory

The Harris-Laplace detected points are scale invariant and work well for isotropic regions that are viewed from the same viewing angle. In order to be invariant to arbitrary affine transformations (and viewpoints), the mathematical framework must be revisited. The second-moment matrix ${\displaystyle \mathbf {\mu } }$ is defined more generally for anisotropic regions:

${\displaystyle \mu (\mathbf {x} ,\Sigma _{I},\Sigma _{D})=\det(\Sigma _{D})g(\Sigma _{I})*(\nabla L(\mathbf {x} ,\Sigma _{D})\nabla L(\mathbf {x} ,\Sigma _{D})^{T})}$

where ${\displaystyle \Sigma _{I}}$ and ${\displaystyle \Sigma _{D}}$ are covariance matrices defining the differentiation and the integration Gaussian kernel scales. Although this may look significantly different from the second-moment matrix in the Harris-Laplace detector; it is in fact, identical. The earlier ${\displaystyle \mu }$ matrix was the 2D-isotropic version in which the covariance matrices ${\displaystyle \Sigma _{I}}$ and ${\displaystyle \Sigma _{D}}$ were 2x2 identity matrices multiplied by factors ${\displaystyle \sigma _{I}}$ and ${\displaystyle \sigma _{D}}$, respectively. In the new formulation, one can think of Gaussian kernels as a multivariate Gaussian distributions as opposed to a uniform Gaussian kernel. A uniform Gaussian kernel can be thought of as an isotropic, circular region. Similarly, a more general Gaussian kernel defines an ellipsoid. In fact, the eigenvectors and eigenvalues of the covariance matrix define the rotation and size of the ellipsoid. Thus we can easily see that this representation allows us to completely define an arbitrary elliptical affine region over which we want to integrate or differentiate.

The goal of the affine invariant detector is to identify regions in images that are related through affine transformations. We thus consider a point ${\displaystyle \mathbf {x} _{L}}$ and the transformed point ${\displaystyle \mathbf {x} _{R}=A\mathbf {x} _{L}}$, where A is an affine transformation. In the case of images, both ${\displaystyle \mathbf {x} _{R}}$ and ${\displaystyle \mathbf {x} _{L}}$ live in ${\displaystyle R^{2}}$ space. The second-moment matrices are related in the following manner:^[3]

${\displaystyle {\begin{aligned}\mu (\mathbf {x} _{L},\Sigma _{I,L},\Sigma _{D,L})&{}=A^{T}\mu (\mathbf {x} _{R},\Sigma _{I,R},\Sigma _{D,R})A\\M_{L}&{}=\mu (\mathbf {x} _{L},\Sigma _{I,L},\Sigma _{D,L})\\M_{R}&{}=\mu (\mathbf {x} _{R},\Sigma _{I,R},\Sigma _{D,R})\\M_{L}&{}=A^{T}M_{R}A\\\Sigma _{I,R}&{}=A\Sigma _{I,L}A^{T}{\text{ and }}\Sigma _{D,R}=A\Sigma _{D,L}A^{T}\end{aligned}}}$

where ${\displaystyle \Sigma _{I,b}}$ and ${\displaystyle \Sigma _{D,b}}$ are the covariance matrices for the ${\displaystyle b}$ reference frame. If we continue with this formulation and enforce that

${\displaystyle {\begin{aligned}\Sigma _{I,L}=\sigma _{I}M_{L}^{-1}\\\Sigma _{D,L}=\sigma _{D}M_{L}^{-1}\end{aligned}}}$

where ${\displaystyle \sigma _{I}}$ and ${\displaystyle \sigma _{D}}$ are scalar factors, one can show that the covariance matrices for the related point are similarly related:

${\displaystyle {\begin{aligned}\Sigma _{I,R}=\sigma _{I}M_{R}^{-1}\\\Sigma _{D,R}=\sigma _{D}M_{R}^{-1}\end{aligned}}}$

By requiring the covariance matrices to satisfy these conditions, several nice properties arise. One of these properties is that the square root of the second-moment matrix, ${\displaystyle M^{\tfrac {1}{2}}}$ will transform the original anisotropic region into isotropic regions that are related simply through a pure rotation matrix ${\displaystyle R}$. These new isotropic regions can be thought of as a normalized reference frame. The following equations formulate the relation between the normalized points ${\displaystyle x_{R}^{'}}$ and ${\displaystyle x_{L}^{'}}$:

${\displaystyle {\begin{aligned}A=M_{R}^{-{\tfrac {1}{2}}}RM_{L}^{\tfrac {1}{2}}\\x_{R}^{'}=M_{R}^{\tfrac {1}{2}}x_{R}\\x_{L}^{'}=M_{L}^{\tfrac {1}{2}}x_{L}\\x_{L}^{'}=Rx_{R}^{'}\\\end{aligned}}}$

The rotation matrix can be recovered using gradient methods likes those in the SIFT descriptor. As discussed with the Harris detector, the eigenvalues and eigenvectors of the second-moment matrix, ${\displaystyle M=\mu (\mathbf {x} ,\Sigma _{I},\Sigma _{D})}$ characterize the curvature and shape of the pixel intensities. That is, the eigenvector associated with the largest eigenvalue indicates the direction of largest change and the eigenvector associated with the smallest eigenvalue defines the direction of least change. In the 2D case, the eigenvectors and eigenvalues define an ellipse. For an isotropic region, the region should be circular in shape and not elliptical. This is the case when the eigenvalues have the same magnitude. Thus a measure of the isotropy around a local region is defined as the following:

${\displaystyle {\mathcal {Q}}={\frac {\lambda _{\min }(M)}{\lambda _{\max }(M)}}}$

where ${\displaystyle \lambda }$ denote eigenvalues. This measure has the range ${\displaystyle [0\dots 1]}$. A value of ${\displaystyle 1}$ corresponds to perfect isotropy.

Iterative algorithm

Using this mathematical framework, the Harris affine detector algorithm iteratively discovers the second-moment matrix that transforms the anisotropic region into a normalized region in which the isotropic measure is sufficiently close to one. The algorithm uses this shape adaptation matrix, ${\displaystyle U}$, to transform the image into a normalized reference frame. In this normalized space, the interest points' parameters (spatial location, integration scale and differentiation scale) are refined using methods similar to the Harris-Laplace detector. The second-moment matrix is computed in this normalized reference frame and should have an isotropic measure close to one at the final iteration. At every ${\displaystyle k}$th iteration, each interest region is defined by several parameters that the algorithm must discover: the ${\displaystyle U^{(k)}}$ matrix, position ${\displaystyle \mathbf {x} ^{(k)}}$, integration scale ${\displaystyle \sigma _{I}^{(k)}}$ and differentiation scale ${\displaystyle \sigma _{D}^{(k)}}$. Because the detector computes the second-moment matrix in the transformed domain, it's convenient to denote this transformed position as ${\displaystyle \mathbf {x} _{w}^{(k)}}$ where ${\displaystyle U^{(k)}\mathbf {x} _{w}^{(k)}=\mathbf {x^{(k)}} }$.