A crater composite feature matching method

Abstract

The present invention discloses a crater composite feature matching method, which relates to the field of image feature matching. The method comprises: first, extracting the crater center point set based on the crater identification result, using the Delaunay algorithm to obtain the crater triangle set, using the triangle and the crater corresponding to its vertex as the feature primitive, and constructing the crater geometric configuration feature network; secondly, using the geometric information of the feature primitive, constructing a shape descriptor that integrates the crater radius constraint; then, selecting the inscribed circle of the triangle as the support area within the feature primitive, and constructing a regional descriptor with scale, rotation and illumination invariance; then, performing weighted fusion on the shape descriptor and the regional descriptor to generate a composite feature descriptor; finally, using the nearest neighbor algorithm to match the generated feature descriptor. The present invention introduces geometric information constraints and regional gradient feature descriptions into the image feature matching method, maintains the matching accuracy when the image undergoes scale, rotation and illumination transformations, and can complete the task of accurately matching small celestial body crater images.

Description

Technical Field

[0001] The present invention relates to the technical field of image feature matching, and in particular to a meteorite crater composite feature matching method. Background Art

[0002] In small-body exploration missions, craters are widely distributed on their surfaces and therefore often used as navigation landmarks. In the field of visual navigation, detecting and identifying craters and matching their image features is a common technique. Traditional small-body image crater matching methods suffer from poor crater detection and recognition due to interference from scale and illumination variations in navigation images acquired during flyby and landing, which in turn affects image matching accuracy.

[0003] A deep learning-based crater detection and recognition method can reduce the impact of illumination variations on the recognition results, resulting in more accurate crater detection results. Crater geometry includes features such as shape and position. Compared to matching methods that rely solely on image features, incorporating geometric information constraints can reduce the impact of factors such as illumination variations on the matching results.

[0004] In view of this, designing a composite feature matching method that can achieve accurate matching under scale, rotation and illumination changes is of great significance to the development and application of image feature matching technology and the development of fields such as autonomous visual navigation. Summary of the Invention

[0005] The technical problem solved by the present invention is that in view of the fact that it is difficult to obtain accurate matching of crater images in the existing technology, an image feature matching method that combines geometric information constraints with regional gradient feature matching algorithms is proposed to effectively improve the matching accuracy of crater images, thereby improving the accuracy of visual navigation methods in small celestial body exploration missions.

[0006] The present invention proposes a crater composite feature matching method, comprising:

[0007] Step A: Based on the crater identification results, extract the crater center point set P = {p1, ..., p n}, use Delaunay triangulation to obtain a triangle set T = {t1, t2, ...}, and use the triangles and the craters corresponding to their vertices as feature primitives to construct a crater geometric configuration feature network;

[0008] In the triangle set T, the kth triangle t k The corresponding characteristic primitive is centered at O i (i=1, 2, 3), radius is r i The three craters C i and the triangle O1O2O3 with the center of the crater as the vertex;

[0009] Step B: Construct a shape descriptor using the geometric information of feature primitives;

[0010] Step C: Select the support region within the feature primitive and construct the region descriptor;

[0011] Step D: Fuse the shape descriptor and the region descriptor to generate a composite feature descriptor;

[0012] Step E: Use the nearest neighbor algorithm to match the generated feature descriptors.

[0013] Furthermore, in step B, the shape descriptor is constructed using the geometric information of the feature primitives, specifically in the following manner:

[0014] The shape descriptor uses S D =[δ1,δ2,δ3,δ4,δ5] T Indicates that:

[0015]

[0016] r1, r2, and r3 are three meteorite craters C ri The radius of the triangle O1O2O3 is r1≥r2≥r3; the vertices O1, O2, and O3 correspond to three sides a, b, and c, respectively, satisfying a≥b≥c; each item is normalized using the parameter a, so that the shape descriptor has scale invariance; at the same time, according to the stable geometric relationship between the side length of the triangle and the radius of the crater, the order of the vectors in the above formula can ensure that the shape descriptor has rotation invariance.

[0017] Furthermore, in step C, a support region is selected within the feature primitive and a region descriptor is constructed, specifically in the following manner:

[0018] Step C1: Select the inscribed circle of the characteristic primitive triangle as the support area G;

[0019] Step C2: Calculate the gradient amplitude m(x, y) and direction α(x, y) of the pixel points in the support region. The total number of pixels in G is n. Divide the pixels in the support region into 8 regions according to their gradient direction α(x, y). Sum the gradient amplitudes of the pixels in each region according to the gradient direction and calculate the gradient histogram. The peak value is used as the main direction of the region descriptor. Traverse all intervals counterclockwise. The sum of the gradient amplitudes in the i-th direction interval can be expressed as:

[0020]

[0021] Step C3: Starting from the interval where the main direction of the region descriptor is located, use M in the order of interval traversal i Construct feature description vector I = [Mt ,...,M8,M1,...,M t-1 ] T , thereby maintaining the rotation invariance of the region descriptor, where the tth interval is the interval where the main direction is located;

[0022] Step C4: Normalize the feature description vector I to ensure that the region descriptor has scale invariance, and finally generate the region descriptor B D :

[0023]

[0024] Among them, M max is the maximum value in I, M min is the minimum value among I.

[0025] Furthermore, in step D, the shape descriptor and the region descriptor are fused to generate a composite feature descriptor, specifically in the following manner:

[0026] Introducing the shape descriptor weight vector w S (x T ,x C ) and the region descriptor weight vector w B (x T ,x C ), realize the weight allocation of shape descriptor and region descriptor in the composite feature descriptor, where x T is a set containing the lengths of the three sides of the triangle, x C is a set of three crater radius values:

[0027] w S (x T ,x C )=[w S1 (x T ,x C ),w S2 (x T ,x C ),...,w S5 (x T ,x C )] T

[0028] w B (x T ,x C )=[w B1 (x T ,x C ),w B2 (x T ,x C ),...,w B8 (xT ,x C )] T

[0029] And the sum of the two sets of weights is equal to 1, that is:

[0030]

[0031] Combine the shape descriptors S obtained in step B and step C D and region descriptor B D , generate a composite feature descriptor D, where the tth interval is the interval where the main direction is located:

[0032]

[0033] Furthermore, in step D, the shape descriptor weight vector w S (x T ,x C ) and the region descriptor weight vector w B (x T ,x C ) is obtained in the following way:

[0034] Define the geometric difference function G S (x T ,x C ) and shape feature saliency quantization function Q S (x T ,x C ):

[0035]

[0036] Among them, x T is a set containing the lengths of the three sides of the triangle, x C is a set of three crater radius values, G(x T ,x C ) has a value range of

[0037] Q S (x T ,x C ) is normalized and the result is used to assign the shape descriptor weight vector:

[0038]

[0039] Let w S (x T ,x C ) Each item takes the same value, w B (x T ,x C) take the same value, that is:

[0040] w S1 (x T ,x C )=w S2 (x T ,x C )=...=w S5 (x T ,x C )

[0041] w B1 (x T ,x C )=w B2 (x T ,x C )=...=w B8 (x T ,x C )

[0042] Combined with the steps in D Get w S (x T ,x C ) and w B (x T ,x C ) values, so that when the geometric features are more significant, w S (x T ,x C ), the greater the weight assigned to each item.

[0043] Compared with the prior art, the advantages and positive effects of the present invention are:

[0044] The present invention proposes a crater composite feature matching method, which performs weighted fusion of geometric information and regional gradient information, and realizes the weight distribution of shape descriptors and regional descriptors according to the prominence of the crater shape, which can improve the matching accuracy of crater images.

[0045] The method proposed in the present invention can obtain accurate matching results when the image has scale, rotation and illumination changes, and can complete the precise matching task of small celestial body images. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] After reading the specific implementation examples with reference to the accompanying drawings, readers will have a clearer understanding of various aspects of the present invention.

[0047] Figure 1 This is an overall flow chart of an embodiment of the present invention;

[0048] Figure 2 Schematic diagram of a triangular characteristic primitive of a meteorite crater in an embodiment of the present invention;

[0049] Figure 3 A schematic diagram of a region construction descriptor in an embodiment of the present invention;

[0050] Figure 4 This is a matching result diagram under scale transformation in an embodiment of the present invention;

[0051] Figure 5 : This is a matching result diagram under rotation transformation in an embodiment of the present invention;

[0052] Figure 6 This is a matching result diagram under illumination transformation in an embodiment of the present invention; DETAILED DESCRIPTION

[0053] The following will be combined with the accompanying drawings to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of them. Therefore, the detailed description of the embodiments of the present invention provided below in the accompanying drawings is not intended to limit the scope of the invention claimed.

[0054] like Figure 1 The overall flow chart of the present invention is shown, which specifically includes the following steps:

[0055] 1. Geometric configuration feature network construction

[0056] The Delaunay algorithm is used to construct a geometric configuration feature network for the obtained crater center: the crater center point set P = {p1, ..., p n} Perform Delaunay triangulation to obtain the crater triangle set T = {t1, t2, ...}, with triangles and corresponding craters as feature primitives, where p i is the center of the i-th crater, and satisfies:

[0057] (1) Maximizing the minimum angle property: The minimum angle of the triangle in the triangle set T is the largest.

[0058] (2) Minimize the circumscribed circle property: The maximum value of the circumscribed circle radius of the triangles in the triangle set T is the smallest.

[0059] 2. Shape descriptor construction

[0060] Figure 2 The schematic diagram of the crater triangle characteristic primitive. In the triangle set T, the kth triangle t k The corresponding characteristic primitive is centered at O i (i=1, 2, 3), radius is r i The three craters C iThe shape descriptor is constructed by the ratios of r1, r2, r3, b, c and a to the vertices O1, O2 and O3.

[0061] The shape descriptor uses S D =[δ1,δ2,δ3,δ4,δ5] T Indicates that:

[0062]

[0063] The crater triangles to be matched in the small-body scene image are similar to the corresponding matching crater triangles in the database, and their radii are proportional to each other. To accurately represent navigation landmark information, a crater radius constraint is introduced into the shape descriptor. Each term is normalized by the parameter a, making the shape descriptor scale-invariant. Furthermore, due to the stable geometric relationship between triangle side length and crater radius, the order of the vectors in the above equation ensures that the shape descriptor is rotationally invariant.

[0064] 3. Region Descriptor Construction

[0065] The construction of the region descriptor includes the following steps, specifically:

[0066] 3.1 Select support area

[0067] Figure 3 A schematic diagram is constructed for the region descriptor. Within the range covered by the feature primitives, the inscribed circle of △ABC is selected as the region descriptor support region G.

[0068] 3.2 Statistical gradient histogram and obtaining the main direction

[0069] Calculate the gradient magnitude m(x,y) and direction α(x,y) of the pixel point (x,y) in the support region G. Divide the pixels in the support region into 8 regions according to their gradient directions. Sum the gradient magnitudes of the pixels in each region according to the gradient directions and calculate the gradient histogram. The peak of the gradient histogram represents the main gradient direction of the pixel point in the support region. Traverse all intervals counterclockwise, and the sum of the gradient magnitudes in the i-th direction interval can be expressed as:

[0070]

[0071] 3.3 Traversing the gradient histogram to construct feature description vector

[0072] Starting from the interval where the main direction of the region descriptor is located, follow the interval traversal order and use M i Construct feature description vector I = [M t ,...,M8,M1,...,Mt-1 ] T , thereby maintaining the rotation invariance of the region descriptor, where the tth interval is the interval where the main direction is located;

[0073] 3.4 Constructing weighted region descriptors

[0074] Normalize the feature description vector I to ensure that the region descriptor has scale invariance, and finally generate the region descriptor B D :

[0075]

[0076] Among them, M max is the maximum value in I, M min is the minimum value among I.

[0077] 4. Fusion generates composite feature descriptors

[0078] Introducing the shape descriptor weight vector w S (x T ,x C ) and the region descriptor weight vector w B (x T ,x C ), realize the weight allocation of shape descriptor and region descriptor in the composite feature descriptor, where x T is a set containing the lengths of the three sides of the triangle, x C is a set of three crater radius values:

[0079] w S (x T ,x C )=[w S1 (x T ,x C ),w S2 (x T ,x C ),...,w S5 (x T ,x C )] T

[0080] w B (x T ,x C )=[w B1 (x T ,x C ),w B2 (x T ,x C ),...,w B8 (x T ,x C )]T

[0081] And the sum of the two sets of weights is equal to 1, that is:

[0082]

[0083] Shape descriptor weight vector w S (x T ,x C ) and the region descriptor weight vector w B (x T ,x C ) is taken as follows:

[0084] Define the geometric difference function G S (x T ,x C ) and shape feature saliency quantization function Q S (x T ,x C ):

[0085]

[0086] Among them, x T is a set containing the lengths of the three sides of the triangle, x C is a set of three crater radius values, G(x T ,x C ) has a value range of

[0087] Q S (x T ,x C ) is normalized and the result is used to assign the shape descriptor weight vector:

[0088]

[0089] Let w S (x T ,x C ) Each item takes the same value, w B (x T ,x C ) take the same value, that is:

[0090] w S1 (x T ,x C )=w S2 (x T ,x C )=...=w S5 (x T ,x C )

[0091] w B1 (x T ,x C )=w B2 (x T ,x C )=...=w B8 (x T ,x C )

[0092] Combined with the steps in step 4 Get w S (x T ,x C ) and w B (x T ,x C ) values, so that when the geometric features are more significant, w S (x T ,x C ), the greater the weight assigned to each item.

[0093] Combine the shape descriptors and region descriptors generated in steps 2 and 3 to generate a composite feature descriptor D. The tth interval is the interval where the main direction is located:

[0094]

[0095] 5. Descriptor Matching

[0096] When performing descriptor matching, the registration effect is measured by calculating the distance between the two descriptor vectors. When registering descriptor D1, if its best registration item is D2, the condition for judging whether the match is established is:

[0097]

[0098] Among them, γ=0.9.

[0099] If the above matching conditions are met, the searched feature descriptors are used as matching pairs. Finally, the RANSAC algorithm is used to eliminate false matches to ensure the accuracy of crater navigation landmark matching.

[0100] The algorithm proposed in this paper is verified using the small celestial body images from NASA's official website. The matching results under scale, rotation and brightness transformation are shown in the following figure: Figures 4 to 6 As shown in Figure 3, the correct matching rate is 88.24% under scale change, 87.88% under rotation change, and 82.14% under illumination change.

[0101] The above specific description further illustrates the purpose, technical solutions and beneficial effects of the invention in detail. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A crater composite feature matching method that integrates geometric information constraints and regional gradient features, the method comprising: Step A: Based on the crater identification results, extract the crater center point set P = {p1, ..., p n }, use Delaunay triangulation to obtain a triangle set T = {t1, t2, ...}, and use the triangles and the craters corresponding to their vertices as feature primitives to construct a crater geometric configuration feature network; In the triangle set T, the kth triangle t k The corresponding characteristic primitive is centered at O i (i=1, 2, 3), radius is r i The three craters C i and the triangle O1O2O3 with the center of the crater as the vertex; Step B: Construct a shape descriptor using the geometric information of feature primitives; Step C: Select the support region within the feature primitive and construct the region descriptor; Step D: Fuse the shape descriptor and the region descriptor to generate a composite feature descriptor; Step E: Use the nearest neighbor algorithm to match the generated feature descriptors; In step B, the shape descriptor is constructed using the geometric information of the feature primitives, specifically in the following manner: The shape descriptor uses S D =[δ1,δ2,δ3,δ4,δ5] T Indicates that: r1, r2, and r3 are the three craters C ri The radius of the triangle O1O2O3 is r1≥r2≥r3; the three sides of the vertices O1, O2, and O3 are a, b, and c respectively, satisfying a≥b≥c; In step C, a support region is selected within the feature primitive and a region descriptor is constructed, specifically in the following manner: Step C1: Define the inscribed circle of the selected feature primitive triangle as the support area G; Step C2: Calculate the gradient magnitude m(x,y) and direction α(x,y) of the pixel point (x,y) in the support region G, where the total number of pixels in G is n; divide the pixels in the support region into 8 regions according to their gradient direction α(x,y), sum the gradient magnitudes of the pixels in each region according to the divided gradient direction and calculate the gradient histogram. The peak value is used as the main direction of the region descriptor, and all intervals are traversed counterclockwise. The sum of the gradient magnitudes of the j-th direction interval can be expressed as: Step C3: Starting from the interval where the main direction of the region descriptor is located, use M in the order of interval traversal j Construct feature description vector I = [M t ,...,M8,M1,...,M t-1 ] T , where the tth interval is the interval where the main direction is located; Step C4: Normalize the feature description vector I and finally generate the region descriptor B D : Among them, M max is the maximum value in I, M min is the minimum value in I; In step D, the shape descriptor and the region descriptor are fused to generate a composite feature descriptor, specifically in the following manner: Introducing the shape descriptor weight vector w S (x T ,x C ) and the region descriptor weight vector w B (x T ,x C ), realize the weight allocation of shape descriptor and region descriptor in the composite feature descriptor, where x T is a set containing the lengths of the three sides of the triangle, x C is a set of three crater radius values: w S (x T ,x C )=[w S1 (x T ,x C ),w S2 (x T ,x C ),...,w S5 (x T ,x C )] T w B (x T ,x C )=[w B1 (x T ,x C ),w B2 (x T ,x C ),...,w B8 (x T ,x C )] T And the sum of the two sets of weights is equal to 1, that is: Combine the shape descriptors S obtained in step B and step C D and region descriptor B D , generate a composite feature descriptor D, where the tth interval is the interval where the main direction is located: In step D, the shape descriptor weight vector w S (x T ,x C ) and the region descriptor weight vector w B (x T ,x C ) is obtained in the following way: Define the geometric difference function G S (x T ,x C ) and shape feature saliency quantization function Q S (x T ,x C ): Among them, x T is a set containing the lengths of the three sides of the triangle, x C is a set of three crater radius values, G(x T ,x C ) has a value range of Q S (x T ,x C ) is normalized, and the result is used to assign the shape descriptor weight vector; Let w S (x T ,x C ) Each item takes the same value, w B (x T ,x C ) take the same value, that is: w S1 (x T ,x C )=w S2 (x T ,x C )=...=w S5 (x T ,x C ) w B1 (x T ,x C )=w B2 (x T ,x C )=...=w B8 (x T ,x C ) Combined with the above Get w S (x T ,x C ) and w B (x T ,x C ) values, so that when the geometric features are more significant, w S (x T ,x C ), the greater the weight assigned to each item.

Images