The SIFT algorithm performs keypoint detection in a multi-scale space, and the Gaussian kernel is the only linear kernel that can complete the scale transformation. We take the convolution of the Gaussian kernel function G(x, y, σ) and the image I(x, y) as the scale space L(x, y, σ) of the image, that is, the LOG (Laplacian of Gaussian) operator: (1)L(x,y,σ)=G(x,y,σ)∗I(x,y)

G(x, y, σ) is a scale-variable Gaussian function. (2)G(x,y,σ)=12πσ2e-(x2+y2)2σ2 Convolution is represented by the symbol ∗, the scale factor is σ, and each pixel position of an image is (x, y). Large scale is used if we need to know the general features of the image, whereas small scale is used if we need to know the detailed features of the image. This is because the larger the scale, the more the image is filtered.

In order to make the key points of detection more stable, David. Lowe et al. proposed to convolved the image with Gaussian difference function to obtain the extreme value of the scale space, i.e., (3)D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))∗I(x,y)=L(x,y,kσ)           -L(x,y,σ) Where k represents the multiple relationship in adjacent images, and k=2 is taken in this paper.

We usually use the response value of the DOG operator to approximate the value of the σ²∇²G operator, because the Gaussian difference operator can be equivalent to the Laplace operator to some extent.

Where, the relationship between D(x, y, σ) and σ²∇²G can be calculated as follows: (4)∂G∂σ=σ∇2G Calculate the difference of Equation (4) as follows: (5)σ∇2G=∂G∂σ≈G(x,y,kσ)-G(x,y,σ)kσ-σ That is: (6)G(x,y,kσ)-G(x,y,σ)≈(k-1)σ2∇2G

D(x, y, σ) can be approximated as σ²∇²G, because the constant (k−1) on the right side of the equation in Formula (3–6) generally does not interfere with the position of the key point. Therefore, this also proves that the DOG operator can be approximated as the LOG operator.

The advantages of the DOG operator are self-evident. (1) Since the information (size and scale) of each layer of images generated in the pyramid is contained in the DOG operator, we can directly use the spatial scale images obtained in Formula (1−3) and obtain the features in the images we need without recalculating the scale again; (2) The LOG operator takes less time to compute Gaussian convolution kernels. This is because compared with the LOG operator, DOG operator only uses one convolution kernels to compute. (3) Since DOG operator can be simplified and approximate to LOG operator, DOG has excellent properties of LOG. For example, LOG operator has invariability on noise, and it will be more stable than other detection methods when carrying out characteristic point detection, such as DoH of Hessian, Harris characteristic point detection, etc.

The strong edge response of DOG operator is caused by the weak response of some feature points, so the key points obtained are not necessarily stable. We can generally use a three-dimensional quadratic function to accurately determine the scale and location information of key points, so that those unstable and low contrast key points are also removed, so that the stability of registration is improved, and the anti-noise performance is also strengthened.

For the difference function D(x, y, σ), expand the second order of Taylor's formula: (7)D(x)=D+∂DT∂xx+12xT∂2D∂x2x Assume the extreme point x^  of x and assume that Equation (7) can be derived and the left-hand side of the equation should be 0. (8)x^=-∂2D-1∂x2∂D∂x By substituting the values obtained in Formula (8) back into Formula (7), feature points with low contrast can be removed: (9)D(x^)=D+12∂DT∂xx^ If the displacement of x^ in all directions is >0.5 relative to the interpolation center point, this point needs to be eliminated because the center point may have shifted to a nearby point. And when |D(x^)|<0.03, the response value point also needs to be deleted because it is too small so it is not stable when disturbed by noise.

Because the image edge will have a great influence on the stability of key points, if some key points are close to or even on the image boundary, then these points cannot be used in image registration because they are extremely unstable. Therefore, in order to prevent the detected key points from being affected by noise, it is not enough to delete only the points with low contrast. Since points on the boundary are easily unstable when affected by noise, and it is difficult for us to accurately determine the position of points on the boundary, it is necessary to delete key points on the boundary.

Strong edge response is the shortcoming that DOG operator cannot overcome. In order to determine whether some key points are at the image boundary, we can use the principal curvature to judge. Hessian matrix H is a method to calculate the magnitude of the principal curvature: (10)H=[DxxDxyDxyDyy] Among them, we only need to find the ratio between the eigenvalues and eigenvalues in the H matrix instead of figuring out the eigenvalues one by one. The value of D can be obtained by calculating the difference of the gray values of the surrounding points. If we set the minimum eigenvalue of H as β = λ_min and the maximum eigenvalue as α = λ_max, then the magnitude and determinant of the trace of H can be expressed by α and β : (11){Tr(H)=Dxx+Dyy=α+βDet(H)=DxxDyy-(Dyy)2=αβ Let α = γβ, α is represented by β and proportional coefficient γ : (12)Tr(H)2Det(H)=(α+β)2αβ=(γβ+β)2γβ2=(γ+1)2γ From Equation (12), we can see that the equation depends only on the ratio of α to β and not on their size. (γ+1)2γ decreases as the ratio of α to β decreases. To minimize (γ+1)2γ, all you need is alpha and beta to be equal. Therefore, we can determine whether the principal curvature is less than the threshold γ through Equation (13). (13)Tr(H)2Det(H)<(γ+1)2γ It works best when γ = 10. That is, we keep the key points where the ratio of alpha to beta is <10, and eliminate the ones that don't meet the requirement.

When describing the key points, it is necessary to assign a main direction to them one by one. The main purpose is to realize the invariance of image rotation, and the main idea is based on the gradient direction and magnitude nature of amplitude. Difference is used to solve the magnitude and direction of the gradient in a circle with a radius of 3 × 1.5σ at the center of the key point: (14){m(x,y)=(L(x+1,y)−L(x−1,y))2+(L(x,y+1)−L(x,y−1))2θ(x,y)=tan−1((L(x,y+1)−L(x,y−1))L(x+1,y)−L(x−1,y)) The sampling is centered on feature points, and histogram is used to describe the gradient and direction. L represents the scale of feature points. The CVD has 36 CVD directions in the histogram, so it ranges from 0 to 360°. It is necessary to use gaussian function to weight the magnitude of the modulus of the gradient when calculating the direction of the histogram.

The direction of the gradient around the key point is represented by the peak value of the histogram. In general, the main direction of the key point is the highest in the histogram, while the direction around the key point is distributed at the rest of the peak value. In order to enhance the robustness of matching, usually there will be multiple auxiliary directions of SIFT feature points, that is, when they are >80% of the peak value, in order to make them have better robustness in image alignment. Although only 15% of these key points may have multiple auxiliary directions, the stability of key point registration has been greatly improved.

Each key point contains three kinds of information: direction, scale, and location. The key points are described so that they are stable without interference from external factors, so that they are invariant to changes in Angle or light. The uniqueness of the descriptor is important in order to achieve higher accuracy in key point matching. Generally, in order to make the key point in a more appropriate scale, it is necessary to sample the surrounding pixels, and then use the normalized correlation algorithm to match the pixel gray. However, simple correlation is very sensitive to affine, non-rigid deformation and 3D perspective changes, which can lead to mismatching of samples.

Feature vectors can be obtained through the structural order of key points. And the SIFT descriptor is a vector that has a lot of dimensions because it computes the gradient histogram. Descriptors need to be computed at the image scale because it is directly related to the scale size. The neighborhood of the key point is divided into 4 × 4 areas with side length of 3σ, i.e., 16 small areas. σ represents the scale size. The key point's neighborhood needs to expand to 152σ because it needs to be interpolated. We get the critical neighborhood and finally 152σ because there's a certain amount of rotation.

In order to make the key point invariant to light changes, the feature vectors need to be uniformly processed, that is, divided into the interval [0, 1]. When unifying key points, the values below 0.2 remain unchanged, and those above 0.2 are fixed as 0.2. In this way, key points are unique. Therefore, SIFT operator has good robustness for noise interference and affine changes, with invariance for image scale, light rotation changes.

In the case that the key point is centered, the direction of the pixel gradient in the neighborhood of the key point is rotated by an Angle θ. The main purpose is to make it invariant to rotation change. Rotate the key point neighborhood to the main direction and divide it into 4 × 4 small blocks of 3 sigma. Compute eight directions in each small block, then sum the results of the directions and get the seed points. The gradient histogram is divided by each sub-block into eight directions with a size of 45°, which is different from the calculation of the main direction of the key point at this point, so the gradient information of the seed point has eight different directions. Since there is a seed point containing eight directional gradients in each of the 16 sub-blocks, we can obtain the SIFT feature vectors with a total dimension of 16 × 8 = 128, and use a Gaussian function with a variance of 6σ to obtain the 128-dimensional feature vectors.

Kd-Tree is a binary Tree structure. Kd-Tree algorithm has many advantages: (1) since the segmentation superplane will keep changing in the process of building KD tree, data points of different clusters can be easily distinguished. (2) Since the resolution of the image can be adjusted when cutting the KD tree space, the height of the tree can be adjusted through the distribution of data points. (3) The cutting surface of KD tree can be adjusted according to different situations. In the process of building KD tree, a data point can be calculated by using the partition hypersurface equation, and the point can be judged to be in the right subtree or the left subtree.

The core of BBF algorithm is to add the search priority part, the rest of the process is similar to the standard KD tree. BBF algorithm is a kind of nearest neighbor optimization algorithm, which gives priority to the points with high matching possibility. Either when all the nodes in the priority queue have been searched or when the preset time runs out, it treats the optimal point of the search as the registration point because of its timeout setting. BBF algorithm can solve the problem that the nearest neighbor cannot handle multidimensional vector after extending KD tree. As a result, the speed of BBF searches is increased, but the accuracy is greatly reduced. Therefore, the results obtained by BBF are relatively good but not the best.

You can see the appropriate benefit of BBF for procedural manipulation of query best-points. The search process can be interrupted at any time through the construction of the priority queue, which is very widely used in high-dimensional data search due to its good results.

Because BBF spends a lot of time searching the nodes of the priority queue, the number of nodes must be manipulated properly. As the number of nodes increases, the search time increases, so the purpose of setting the number of nodes to 200 in this article is to control the search time. In addition, the order of searching leaf nodes depends on the structure of KD tree, so the location of searching node may be ignored even if the location of node storage is known. This problem can be solved by using the distance between the node being searched and the current node to query the node. If the Euclidean distance between the target point and the query point is smaller than the distance between the target point and the nearest matching point, then the nearest matching point is the query point.

Before using RANSAC to remove the external points to calculate the transformation matrix parameters, a preliminary selection is made through the distance ratio between the nearest neighbor and the next nearest neighbor of the keypoint. This can reduce the workload of RANSAC and improve the purity of the matching point pairs using RANSAC. The K-nearest neighbor algorithm used in this paper can avoid large amounts of calculation caused by the exhaustive search method. The K value is taken as 2. Finally, each keypoint will get two matching points with the smallest and the next smallest distance. Let the eigenvector of the key points p be V_p, the vector of the closest matching point q in the first image is V_q and the vector of the next closest matching point q̄ in the first image is Vq̄, then the description matching pi(Vq,Vq̄),i=1,2,⋯ can be used. n is the number of key points extracted from the second image.

Assuming Vq̄ that are the adjacent key points Vq̄ in the first image. To make p and q become a set of accurate registration pairs, R needs to be smaller than the critical value T_R. At this time, it is represented by p ↔ q, and vice versa p ← | → q. Then the distance ratio is expressed as follows. (15)R=D(Vp,Vq)D(Vp,Vq̄) (16)Vq=argmin{D(Vp,Vi)|Vi∈B,i=1,⋯,m} (17)Vq̄=argmin{D(Vp,Vj)|Vj∈B,j=1,⋯,m,Vj≠Vq} (18){p↔q if R<TRp←|→qotherwise The experimental results show that when the critical value T_R of Equation (18) is 0.65, the accuracy and the matching logarithm of keypoints can get better results. If the set A is used to represent the final registration pair, then, A={pi(Vq,Vq̄)|i=1,2,…m1}. m₁ is the number of key point registration pairs in the corresponding second frame of the first frame.

Then, in the second frame image, the nearest and next closest registration points of all the key points matching in the first frame image are searched. According to the above process, the error points whose R larger than the threshold value T_R are eliminated. The final keypoint registration pair is represented by a set B, that is, B={qi(Vp,Vp̄)|i=1,2,⋯m2}, m₂ is the number of key point registration pairs in the corresponding first frame of the second frame image.

Finally, all the data are searched. For a certain registration pair (p_i, q_j) in the set A, if a matching pair (q_j, p_i) can be found in the set B, this point is accepted. The set C is used to store the keypoint pairs obtained through bidirectional registration, where C={pi(Vq,Vq̄)|i=1,2,⋯m3}. m₃ is the number of keypoint pairs obtained through bidirectional registration.

RANSAC is a robust parameter estimation algorithm (Triggs et al., 1999) proposed by Fischler and Bones. At first, people used this algorithm to estimate camera motion, and now it is widely used in parameter estimation. RANSAC can accurately get the correct registration relationship and eliminate the wrong matching points.

Firstly, the ratio of the nearest neighbor and the next nearest neighbor is used to express the registration results. The purpose of holography matrix detection is to find the most matching key points and eliminate those that do not meet the requirements. Therefore, we have the initial set of matching points.

For further purification of the key points after two-way registration, RANSAC's main approach is to assume a model applicable to all the key points that are correctly registered. The assumed interior points can be calculated from all the unknown parameters. All the data points are randomly sampled, and some points are set as interior points to form a random subset. Finally, random sampling is carried out until we find a transformation parameter model that can make the number of interior points maximum. Therefore, we need to remove the interference caused by mismatching points when we get interior points.

Generally, the transformation model with eight parameters is mainly used in splicing, the matrix transformation with 8 degrees of freedom. Because a pair of registration points can support two equations to establish the relationship, four non-collinear matching pairs are needed to calculate the parameters (Zhang et al., 2021).

The steps of RANSAC are as follows:

Four pairs are randomly selected in the registration set to solve the transformation matrix.

In solving the remaining registration pairs' matrix transformation, we need to use the transformation parameters and calculate the coordinate distance between points and point pairs. The inner point is the point within the error range, and the other points are called the outer point. In the first random sampling process, the maximum number of interior points is X, and we set MAX as MAX = X in the first arbitrary sampling process. If the ratio X is larger than MAX, then MAX = X.

This process is continued until the number of interior points is no longer increased. And the point set larger than the preset critical value of the number of interior points is regarded as the largest interior point set. The matrix parameter model is solved by using the least square method.

It is also necessary to set the maximum number of sampling, the critical value of distance, and the minimum number of points.

If there are n matching pairs in the registration, we can get the number of interior points is n(1 − ε). Generally, the number of interior points will not change after some iterations.