Method for small-scale lunar impact crater extraction based on image spatial domain and frequency domain features

Abstract

The application relates to a small-scale moon surface impact crater extraction method based on image space domain and frequency domain features, which comprises the following steps: acquiring LROC NAC images, and performing impact crater seed detection by using mathematical morphology to obtain impact crater seed points; based on the impact crater seed points, impact crater candidate region extraction is performed by using space domain features and based on statistical morphology constraint path set; based on the impact crater seed points, impact crater candidate region extraction is performed by using frequency domain features and based on structure consistency constraint multi-path; cross verification of the impact crater candidate regions extracted based on the space domain and the frequency domain is performed based on a difference factor, false extraction impact craters are removed, and the final impact crater extraction result is determined. Compared with the prior art, the application has the advantages of high detection result accuracy.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of small-scale lunar impact crater detection, and in particular to a small-scale lunar impact crater extraction method based on image spatial domain and frequency domain features. BACKGROUND

[0002] Impact craters are the main features of the lunar surface and play an important role in lunar exploration. For example, impact craters can be used as an effective tool to estimate the absolute age of the lunar surface and help understand the morphological features and geological processes of the lunar surface. Large-scale impact craters with a diameter greater than 10 kilometers shape the basic morphology of the moon, while small lunar craters with a diameter less than 100 meters are constantly forming and changing the lunar surface morphology. Accurate extraction of small-scale impact craters on the lunar surface has important scientific significance and is also important for engineering design decisions. The present application focuses on detecting small impact craters from Lunar Reconnaissance Orbiter Camera Narrow Angle Camera (LROC NAC) images.

[0003] At present, many effective impact crater detection methods have been proposed. As for these methods, various types of features are used. The first widely used feature is the geometric feature of the impact crater. Geometric features refer to the edges and shapes of impact craters. Generally, the shape of an impact crater is always modeled as a circle or an ellipse. Impact crater detection is performed using shape information based on template matching methods, however, these methods fail when the impact crater edges are blurred or severely degraded. Utilizing texture information features is another important means of detecting impact craters. Typical texture features include highlight-shading features, Haar-like features, histogram of oriented gradients (HOG), etc. Among different features, the highlight-shading feature is an important feature and has been used in various methods. The highlight-shading feature describes the gray level distribution of the shadow and sunlit surface under certain lighting conditions. For machine learning methods, the detection of impact craters is usually based on the idea of classification, and widely used classifiers include support vector machines, AdaBoost-based models, decision trees, etc. In addition, methods based on deep convolutional neural networks play an important role in impact crater detection. Commonly used deep networks include region-based convolutional neural networks such as Faster R-CNN and Mask R-CNN. However, in general, at least hundreds or thousands of impact crater samples are needed to train machine learning or deep learning models. In terms of small-scale impact crater detection, both machine learning and deep learning-based methods face the challenge of insufficient training samples.

[0004] Generally, under certain lighting conditions, a crater can be considered as composed of highlight and shadow features. However, due to the influence of terrain, lighting and degradation, there is no consistent geometric feature or reflectivity contrast between the highlights and shadows of different craters. In order to overcome these challenges, CraterDANet: A Convolutional Neural Network for Small-Scale Crater Detection via Synthetic-to-Real Domain Adaptation proposes an efficient highlight-shadow feature descriptor which can adaptively detect craters with different morphologies and highlight-shadow contrasts. The disadvantage of this method is that only one path operator is used when detecting highlight (or shadow) features, which is easily affected by noise, resulting in incomplete or incorrect detection. SUMMARY

[0005] The purpose of the present application is to provide a small-scale lunar surface crater extraction method based on image spatial domain and frequency domain features, and to avoid the problem of being easily disturbed by noise caused by using only one path operator by proposing a statistical morphology constrained path set, thereby improving detection accuracy. At the same time, in order to fully utilize the highlight-shadow feature to detect craters, structural features are introduced to detect the highlight-shadow features of small craters with adaptive scales, thereby improving the robustness to different highlight-shadow contrasts.

[0006] The purpose of the present application can be achieved by the following technical solutions:

[0007] A small-scale lunar surface crater extraction method based on image spatial domain and frequency domain features, comprising the following steps:

[0008] Obtaining LROC NAC images;

[0009] Based on the LROC NAC images, using mathematical morphology to detect crater seeds, obtaining crater seed points;

[0010] Based on the crater seed points, using spatial domain features, and based on a statistical morphology constrained path set, extracting crater candidate regions;

[0011] Based on the crater seed points, using frequency domain features, and based on structural consistency constrained multi-path, extracting crater candidate regions;

[0012] Based on the difference factor, cross-verification of the crater candidate regions extracted by spatial domain and frequency domain, eliminating false extraction craters, and determining the final crater extraction result.

[0013] The crater seed points are pixels with the largest local gradient and gradient direction along the lighting direction.

[0014] The impact crater seed detection using mathematical morphology includes the following steps:

[0015] The gradient size and direction of the current pixel point are obtained, denoted as G mag and G dir , respectively.

[0016] The pixel with local maximum gradient is detected based on the morphological dilation operation, denoted as C ini '.

[0017]

[0018] Wherein, G mag dilate is the gradient amplitude graph after morphological dilation operation, B = strel ('square', r) is a structural element, r represents the size of B, represents the dilation operation.

[0019] The extraction of impact crater seed points based on illumination direction constraint:

[0020] S ini = {(x i ,y i ) | (x i ,y i ) ∈ C i ' ni ,G dir (x i ,y i )- θ ≤ τ}

[0021] Wherein, τ is defined as the preconfigured difference between the light direction and the gradient direction, θ is the light angle, S ini represents the pixel that satisfies the gradient size and direction constraint at the same time, i.e. the impact crater seed point.

[0022] The impact crater candidate region extraction based on the statistical morphological constraint path set includes the following steps:

[0023] The image space domain and its path representation are defined based on graph theory, wherein the image space domain is defined as a directed and weighted acyclic graph, the path extension direction is along the light direction, perpendicular to the light direction and 45° direction with both, and the path is defined as a pixel sequence, which constitutes a path only when and if the adjacent elements in the path satisfy the path operator construction condition.

[0024] With the current impact crater seed point as the starting point, the pixel points are processed along the light direction and the inverse light direction based on the dilation and erosion operations of statistical morphological, respectively.

[0025] The path operator construction condition is determined based on the statistical morphological processing result, a corresponding path operator is constructed, bright features and dark features of the impact crater are detected respectively, and path sets of the bright features and the dark features are obtained;

[0026] The candidate region of the impact crater is determined based on the path sets of the bright features and the dark features.

[0027] The statistical morphological dilation and erosion operation on the pixel point is represented as:

[0028]

[0029] Wherein, I'(x i ) represents the statistical morphological result of the pixel x i . I(x i ) represents the image gray value of the pixel x i , x i r is the neighborhood pixel of the pixel x i , r = 1, 2, …, 8, I(x i r ) represents the image gray value of the pixel x i r , Z(x i ) is a normalization constant; δ r is an indicator function, that is, for the morphological dilation operation, if I(x i r ) ≥ I(x i ), δ r = 1, otherwise δ r = 0, for the morphological erosion operation, if I(x i r ) ≤ I(x i ), δ r = 1, otherwise δ r = 0.

[0030] The path operator construction condition is:

[0031] For the detection of the bright features of the impact crater, the pixel points satisfying I'(x j )-I(x i ) ≥ 0 are classified into the bright feature path;

[0032] For the detection of the dark features of the impact crater, the pixel points satisfying I'(x j )-I(x i ) ≤ 0 are classified into the dark feature path.

[0033] The method for determining the candidate region of the impact crater based on the structural consistency constraint multi-path includes the following steps:

[0034] A path representation in frequency domain is defined based on a directed and unweighted graph, wherein the path extending directions are along the illumination direction, the inverse illumination direction and the orthogonal direction of the two directions, and the path is defined as a pixel sequence, and the path is formed when and only when the adjacent elements in the path satisfy the structural consistency constraint multi-path operator construction condition;

[0035] The image domain information of the image is converted into frequency domain information based on a log-Gabor filter;

[0036] A local phase is determined based on the frequency domain information;

[0037] A structural consistency constraint multi-path operator construction condition is determined based on the local phase with the current seed point as the center, a corresponding path operator is constructed, and the bright feature and the dark feature of the impact crater are detected respectively, so as to construct a structural consistency multi-path set of the bright feature and the dark feature;

[0038] A candidate region of the impact crater is determined based on the structural consistency multi-path set, and a structural consistency index is calculated based on the Euclidean distance between the paths, and if the structural consistency index is less than a preconfigured threshold value, it is indicated that the candidate region covers the whole.

[0039] The log-Gabor filter is defined as:

[0040]

[0041] wherein w is the frequency of the current pixel point, w0 is the center frequency of the filter, and σ0 is the filter bandwidth;

[0042] The frequency domain information processed by the log-Gabor filter is denoted as f M (x) M (x)=[f e (x),f o (x)] T wherein f e (x) represents an odd component with an asymmetric property in the frequency domain, and f o (x) represents an even component with a symmetric property in the frequency domain.

[0043] The local phase is:

[0044]

[0045] If the local phase satisfies the structural consistency constraint multi-path operator construction condition the corresponding pixel point is attributed to the path of the bright feature;

[0046] If the local phase satisfies the structural consistency constraint multi-path operator construction condition the corresponding pixel point is attributed to the path of the dark feature.

[0047] The difference factor is defined as:

[0048]

[0049] wherein C i I = (v i I , u i I , r i I ) represents a crater candidate region extracted based on statistical morphological constraint path set, (v i I , u i ) is the crater center coordinate thereof, and r i I is the radius; C i F = (v i F , u i F , r i F ) represents a crater candidate region extracted based on structural consistency constraint multi-path, (v i F , u i F ) is the crater center coordinate thereof, and r i F is the radius; R max and R min respectively represent the maximum and minimum radius of C i I and C i F ; d(v i (C i I ~ C i F , u i (C i I ~ C i F )) represents the distance between the crater centers of C i I and C i F .

[0050] If the difference factor CDF is less than a preconfigured threshold value, the current crater candidate region is considered as a correctly extracted crater, otherwise, the current extracted candidate region is eliminated.

[0051] Compared with the prior art, the present application has the following beneficial effects:

[0052] The present application provides a method for detecting small-scale impact craters on the moon, which involves two impact crater descriptors based on image spatial domain and frequency domain: statistical morphology constraint path set and structure consistency constraint multi-path. The statistical morphology constraint path set is constructed in the image spatial domain and consists of two main parts: statistical morphology is used to resist local gray abnormal value to ensure the integrity of the detected impact craters, and the design of the path set is to ensure the uniformity of the detected highlights (or shadows) through the mutual restriction between multiple paths, thereby avoiding false detection caused by edge blur; the structure consistency constraint multi-path operator is a structure feature constructed based on the local phase derived from morphological signals, which can improve the robustness to different highlight-shade contrast. Therefore, the present application can combine the advantages of image features and structure features, to some extent, utilize the complementarity of the two features, effectively detect small-scale impact craters on the moon surface, and improve the detection accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0053] Figure 1 The flow chart of the method of the present application is shown in the figure;

[0054] Figure 2 The directed weighted graph and the path extension schematic diagram are shown in the figure;

[0055] Figure 3 The LROC NAC images of data (A) and data (B) in the embodiment of the present application are shown in the figure;

[0056] Figure 4 The impact crater extraction results of the LROC NAC images of data (A) and data (B) in the embodiment of the present application are shown in the figure. DETAILED DESCRIPTION

[0057] The present application will be described in detail below in combination with the drawings and specific embodiments. The present embodiment is implemented on the premise of the technical solution of the present application, and gives a detailed implementation manner and specific operation process, but the protection scope of the present application is not limited to the following embodiments.

[0058] The extraction of small-scale impact craters on the moon has important scientific and engineering significance. Compared with large-scale impact craters, small-scale impact craters are more in number and more complex in shape. Therefore, the extraction of small-scale impact craters on the moon is still a challenging task. The present invention proposes a method for extracting small-scale impact craters on the moon by combining image spatial domain and frequency domain features, which mainly consists of three key steps: first, a series of mathematical morphology operations are used to locate impact crater seed points, wherein each seed point represents a possible location of an impact crater. Then, image spatial and frequency domain features are used to extract impact crater candidate regions. In the image domain, based on a specific directed weighted graph, a path set operator under statistical morphology constraint is constructed to detect impact crater regions by using statistical morphology constraint. The operator can effectively adapt to the complex morphological features of impact craters to detect complete and correct impact crater regions. In the frequency domain, based on the defined directed acyclic graph, a structure consistency constraint multi-path operator is constructed to extract impact craters by using the local phase of Gabor-Log filter detection. Finally, at each seed position, an impact crater difference factor operator is defined to cross-validate the impact craters extracted based on spatial and frequency domain features to obtain the final impact crater detection result.

[0059] Specifically, the present embodiment provides a method for extracting small-scale impact craters on the moon based on image spatial domain and frequency domain features, as shown in the following steps: Figure 1

[0060] (1) Obtain LROC NAC images.

[0061] (2) Based on the LROC NAC images, use mathematical morphology to detect impact crater seeds to obtain impact crater seed points.

[0062] In the present embodiment, the impact crater seed point is defined as the pixel with the largest local gradient and gradient direction along the light direction. This definition is because for an impact crater, the pixel with the largest gradient is always located in the high-light-shadow transition area or its edge. However, their gradient directions are opposite. Specifically, the gradient direction of the high-light-shadow transition area is along the light direction. Therefore, the impact crater seed can be determined by considering the size and direction of the local gradient. The present embodiment uses a series of mathematical morphology operations to detect these pixels, and the specific extraction process is as follows:

[0063] (21) Obtain the gradient size and direction of the current pixel point, denoted as G mag and G dir , respectively.

[0064] (22) Detect the pixel with the largest local gradient based on the mathematical morphology dilation operation, denoted as C ini ’:

[0065] ​

[0066] where G mag dilate is the gradient magnitude image after morphological dilation operation, B = strel('square', r) is the structuring element, and r represents the size of B, represents the dilation operation;

[0067] (23) Extraction of impact crater seed points based on illumination direction constraint:

[0068] S ini = {(x i ,y i ) | (x i ,y i ) ∈ C i ' ni ,G dir (x i ,y i )- θ ≤ τ}

[0069] where τ is defined as the preconfigured difference between the illumination direction and the gradient direction, and θ is the illumination angle, which ensures that the gradient direction of the corresponding pixel and the illumination direction are within a certain range, S ini represents the pixels that satisfy both the gradient size and direction constraints, i.e., the impact crater seed points.

[0070] (3) Based on the impact crater seed points, spatial domain features are used to extract impact crater candidate regions based on statistical morphological constraints on path sets.

[0071] The highlight-shade feature of impact craters in the image is essentially the manifestation of the morphological features of impact craters under certain lighting conditions. Morphological features refer to the fact that impact craters are a kind of landform with "concave features". In the same image, although the contrast of highlights and shadows is different between different impact craters, their gray level change trend is similar, i.e., from the seed point, along (or against) the illumination direction, the gray value increases (or decreases). In terms of the gray value change of adjacent pixels, there may be local anomalies, but the overall gray level change trend is not affected. Local anomalies are caused by topographic changes, morphological degradation, or being covered by other impact craters, etc. Based on these considerations, this embodiment utilizes statistical morphological and path set strategies. On the one hand, the path can adapt to the changes in the highlight and shadow areas, and on the other hand, the design of statistical morphological and path set can improve the robustness of the detector to noise.

[0072] Specifically, the following steps are included:

[0073] (31) Define the image spatial domain and its path representation based on graph theory.

[0074] The image space domain D is defined as a directed and weighted acyclic graph, and the path extension directions are along the light direction, perpendicular to the light direction and 45° to both directions, as shown in Figure 2 Fig. 1. The dashed arrows represent the directions in which the path can be extended. A path is defined as a pixel sequence π(x i ) = (V(x i ), V(x i+1 ), V(x i+2 ),..., V(x i+k ),..., V(x i+m )), which can be constructed if and only if the adjacent elements V(x i ) and V(x i+1 ) satisfy the path operator construction condition. The path set is the collection of multiple paths obtained by different rules with the same pixel point as the starting point.

[0075] (32) Taking the current impact crater seed point as the starting point, the pixel points are processed along the light direction and the reverse light direction based on the dilation and erosion operations of statistical morphology, respectively;

[0076]

[0077] where I'(x i ) represents the statistical morphology result of pixel x i ; I(x i ) represents the image gray value of pixel x i ; x i r is the neighborhood pixel of pixel x i , r = 1, 2,..., 8, I(x i r ) represents the image gray value of pixel x i r , Z(x i ) is a normalization constant; δ r is an indicator function, that is, for the morphological dilation operation, if I(x i r ) ≥ I(x i ), δ r = 1, otherwise δ r = 0, for the morphological erosion operation, if I(x i r ) ≤ I(x i ), δ r = 1, otherwise δ r = 0.

[0078] (33) Determine path operator construction condition based on statistical morphological processing result, construct corresponding path operator, detect bright feature and dark feature of impact crater respectively, and obtain path set of bright feature and dark feature.

[0079] The path operator construction condition is:

[0080] For detection of bright feature of impact crater, pixel points satisfying I'(x j )-I(x i )≥0 are classified into bright feature path.

[0081] For detection of dark feature of impact crater, pixel points satisfying I'(x j )-I(x i )≤0 are classified into dark feature path.

[0082] Let the path set of extracted bright feature and dark feature of impact crater be denoted as Π H I (S i ), Π S I (S i ) respectively, then for each impact crater seed point, there is path set Path I (S i )={Π H I (S i ), Π S I (S i )}.

[0083] (34) Determine candidate region of impact crater based on path set of bright feature and dark feature.

[0084] According to path operator, candidate region of impact crater can be extracted, specifically:

[0085]

[0086]

[0087]

[0088]

[0089]

[0090] wherein, X Path I (Si) and Y Path I (Si) represent image coordinates of each element in horizontal and vertical directions of path respectively.

[0091] Candidate impact craters extracted based on statistical morphological constraint path sets are denoted as follows: in, The coordinates of the impact crater center are: Let be the radius.

[0092] (4) Based on the seed points of the impact crater, the candidate regions of the impact crater are extracted by using frequency domain features and multi-path based on structural consistency constraints.

[0093] This embodiment uses Gabor-log filters to analyze structural features obtained from the frequency domain to extract candidate impact craters. Structural features have been used in various image analysis tasks, such as image matching and object detection. A key advantage of structural features is that they primarily describe the topological and geometric information of an image and are independent of image intensity. For impact crater detection, topological and geometric information can correspond to the distribution of specular-shadow features and the shape of the impact crater, respectively. Therefore, this embodiment proposes using Gabor-log filtering results to construct structural features for impact crater detection. The proposed structural feature impact crater detection operator is named the structural consistency-constrained multipath operator, which is robust to different specular-shadow contrasts. Specifically, it includes the following steps:

[0094] (41) A frequency domain path representation is defined based on a directed unweighted graph, where the path extension directions are along the illumination direction, against the illumination direction, and the orthogonal direction of these two directions, respectively. The path is defined as π(x i |s k )=(V(x i ),V(x i+1 ),V(x i+2 ),...,V(x i+k ),...,V(x i+m The adjacent element V(x) is true if and only if... i ) and V(x i+1 When the structural consistency constraint multipath operator construction conditions are met, a path can be formed:

[0095] (42) Converting image domain information into frequency domain information based on log-Gabor filtering:

[0096]

[0097] Where w is the current pixel frequency, w0 is the center frequency of the filter, and σ0 is the filter bandwidth.

[0098] The frequency domain information after log-Gabor filtering is denoted as f. M (x), f M (x)=[f e (x),fo (x)] T where f e (x) represents odd components with asymmetric properties in the frequency domain, and f o (x) represents even components with symmetric properties in the frequency domain.

[0099] (43) Determine the local phase based on f e (x) and f o1 (x):

[0100]

[0101] (44) Determine the structure consistency constraint multi-path operator construction condition based on the local phase with the current seed point as the center, construct the corresponding path operator, and detect the bright features and dark features of the impact crater respectively, so as to construct the structure consistency multi-path set of the bright features and dark features;

[0102] If the local phase satisfies the structure consistency constraint multi-path operator construction condition , the corresponding pixel point is assigned to the path of the bright feature;

[0103] If the local phase satisfies the structure consistency constraint multi-path operator construction condition , the corresponding pixel point is assigned to the path of the dark feature.

[0104] The constructed multi-path is represented as: Π(x i |s k ) = {π1 H (x i |s k ), π2 H (x i |s k ), π1 S (x i |s k ), π 2’ S (x i |s k )}, and the path set of the bright features and dark features obtained in the frequency domain is respectively denoted as Π H F (S i ), Π S F (S i ). Then, for each impact crater seed point, there is a path Path F (S i ) = {Π H F (S i ), Π S F (Si )}.

[0105] (45)Based on the structure consistency multi-path set, a candidate area of the impact crater is determined, and a structure consistency index is calculated based on the Euclidean distance between the paths. If the structure consistency index is less than a preconfigured threshold value, it indicates that the candidate area covers the whole.

[0106] According to the path operator, a candidate area of the impact crater can be obtained. Specifically,

[0107] l i I =max(X Path( F Si) )-min(X Path F (Si) )

[0108] w i I =max(Y Path F (Si) )-min(Y Path F (Si) )

[0109] v i I =(max(X Path F (Si) )+min(X Path F (Si) ) / 2

[0110] u i I =(max(Y Path F (Si) )+min(Y Path F (Si) )) / 2

[0111] r i I =(l i I +w i I ) / 2

[0112] wherein, X Path F (Si) andY Path F (Si)respectively represent the image coordinates in horizontal and vertical directions of each element in the path.

[0113] The candidate impact crater extracted using the structural consistency constraint multi-path is denoted as C i F = (v i F , u i F , r i F ), wherein (v i F , u i F ) is the impact crater center coordinate, and r i F is the radius.

[0114] After the candidate region is extracted, the structural consistency index ISC = D (Π i F (x i | s k ) ~ Π j F (x j | s k )) is calculated, wherein D (·) is the Euclidean distance operator. It is judged whether ISC is less than or equal to the preconfigured threshold value. If yes, it indicates that the candidate region is complete, otherwise, it indicates that it is not complete, and the missed extraction impact crater will affect the final extraction result. In this embodiment, the candidate region should be as complete as possible.

[0115] (5) Cross verification of the impact crater candidate region extracted in the spatial domain and the frequency domain based on the difference factor, elimination of the wrong extraction impact crater, and determination of the final impact crater extraction result.

[0116] The difference factor is defined as:

[0117]

[0118] wherein R max and R min respectively represent the maximum and minimum radius of C i I and C i F ; d (v i (C i I ~ C i F , u i (C i I ~ C i F )) represents Ci I and C i F the distance between the centers of the impact craters.

[0119] If the difference factor CDF is less than a preconfigured threshold th, the current impact crater candidate region is considered as a correctly extracted impact crater, otherwise, the current extracted candidate region is rejected.

[0120] To verify the robustness and universality of the method described in the present application, the LROC NAC images of Chang'e-3 landing area and Chang'e-4 landing area are selected as experimental images, which are denoted as data (A) and data (B) respectively. Chang'e-3 probe landed in the northern part of Imbrium basin at west longitude 19.51° and north latitude 44.12°, which is a typical flat area of the moon's internal topography. Chang'e-4 probe landed in the Von Kármán impact crater in the South Pole-Aitken basin (SPA) on the far side of the moon at east longitude 177.59° and south latitude 45.46°, which has more complex terrain and more dense impact craters. The image of each selected test point is shown in Fig. 1, and the image information used is shown in Table 1. Figure 3

[0121] Table 1 Basic image information of data (A) and (B)

[0122]

[0123] To quantitatively evaluate the method of the present application, the impact craters manually identified in the ArcGIS software using the CraterTools plug-in are used as the ground truth to evaluate the proposed algorithm. The evaluation indicators used are as follows:

[0124]

[0125]

[0126]

[0127] Among them, true positive (TP) represents the number of correctly detected impact craters, false positive (FP) represents the number of false detections, and false negative (FN) represents the number of missed detections. The correct detection rate (TDR) represents the integrity of impact crater detection, and its optimal value is 1. The false detection rate (FDR) represents the correctness of impact crater detection, and its optimal value is 1. The detection rate (DR) represents the quality of the detection result, considering both integrity and correctness indicators, and its optimal value is 1.

[0128] Figure 4 ​(a) and (b) show the crater detection results of the proposed method on data (A) and (B), respectively, where the circles represent the extracted craters. Table 2 shows the quantitative evaluation results of the proposed method. For data (A), 2703 craters were detected from the image pyramid, of which 2602 were correct and 101 were false. The false detections were mainly caused by the gullies, which have similar image features to the craters. The missed craters were mainly large-scale craters with high degradation. Therefore, even though the image pyramid was built, the high-light-shadow feature could not be effectively expressed. For data (B), 13159 craters were detected from the image pyramid, of which 13063 were correct and 96 were false. The false detections were mainly caused by the concave grooves, which have similar image features to the craters. The missed craters were mainly covered by other smaller craters, which severely damaged their image features. The proposed method performed well in detecting small-scale craters, which is because the image features of these craters are highly consistent with the proposed feature descriptor. This indicates that the proposed feature descriptor is effective and robust for detecting small-scale craters. Figure 4 (a) can be seen that the proposed method can detect almost all the craters in the image, and it performs particularly well in detecting small-scale craters, with a TDR of 0.9513. The smallest crater detected by the proposed method is about 5.5 meters in diameter, and the largest crater is about 118.5 meters in diameter. From Figure 4 (a) can be seen that the terrain of data (A) is flat, and the density of craters is relatively low. As listed in the second row of Table 2, the total number of artificial ground truth craters in this area is 2654, and 2703 craters are detected, of which 2602 are correct and 101 are false. Most of the false detections are caused by gullies, which have similar image features to craters. The missed craters are mainly large-scale craters with high degradation. Therefore, even though the image pyramid is built, the high-light-shadow feature cannot be effectively expressed. From Figure 4 (b) can be seen that the craters in data (B) are densely distributed. The quantitative detection results are listed in the third row of Table 2. The number of artificial ground truth craters in this area is 13175, and the total number of detected craters is 13159, of which 13063 are correct and 96 are false. The excellent performance of the proposed method is due to two aspects. On the one hand, the flat terrain, uniform illumination, and small amount of non-crater topography reduce the false detection rate of craters. On the other hand, the obvious and consistent high-light-shadow feature of the craters improves the accuracy of crater detection. A small number of missed craters are due to being covered by other smaller craters, which severely damages their image features. These few false detections are mainly located in the grooves, which have similar image features to craters. According to the results, it can be found that the proposed method has an advantage in detecting small-scale craters. This is because the image features of these craters are highly consistent with the proposed feature descriptor, indicating that the proposed feature descriptor is effective and robust for detecting small-scale craters.

[0129] Table 2 Quantitative evaluation of the proposed method on data (A) and (B)

[0130]

[0131] In summary, the embodiment proposes a method for detecting small-scale impact craters on the moon, which involves two impact crater descriptors based on image spatial domain and frequency domain: statistical morphology constrained path set and structure consistency constrained multi-path. The statistical morphology constrained path set is constructed in the image spatial domain and consists of two main parts. Statistical morphology is used to resist local gray abnormal values to ensure the integrity of the detected impact craters, and the design of the path set is to ensure the uniformity of the detected highlights (or shadows) through the mutual restriction between multiple paths, thereby avoiding false detection caused by edge blur. The structure consistency constrained multi-path operator is a structure feature constructed based on the local phase derived from morphological signals, aiming to improve the robustness to different highlight-shade contrasts. According to the experimental results, it can be found that the image features and the structure features have their own advantages, and to some extent, the features of the two have complementarity, which can effectively detect small-scale impact craters on the moon surface.

[0132] The preferred embodiments of the present application are described in detail above. It should be understood that those of ordinary skill in the art can make modifications and variations without creative effort based on the concept of the present application. Therefore, any technical solutions obtained by logical analysis, reasoning or limited experiments based on the prior art according to the concept of the present application shall be within the protection scope defined by the claims.

Claims

1. A small-scale lunar impact crater extraction method based on image spatial domain and frequency domain features, characterized in that, The method comprises the following steps: acquiring LROC NAC images; based on the LROC NAC images, using mathematical morphology to detect impact crater seeds to obtain impact crater seed points; based on the impact crater seed points, using spatial domain features, and based on statistical morphology, constraining path sets to extract impact crater candidate regions; based on the impact crater seed points, using frequency domain features, and based on structural consistency, constraining multiple paths to extract impact crater candidate regions; based on a difference factor, cross-verification of the impact crater candidate regions extracted in the spatial domain and the frequency domain is performed to eliminate false extraction of impact craters, and a final impact crater extraction result is determined.

2. The method according to claim 1, wherein, The impact crater seed point is a pixel with the largest local gradient and a gradient direction along an illumination direction.

3. The method of claim 1, wherein, The impact crater seed detection using mathematical morphology comprises the following steps: Obtain the gradient size and direction of the current pixel point, respectively denoted as G mag and G dir ; Detecting pixels with local maximum gradient based on mathematical morphological dilation operation, denoted as C ini ’ : where G mag dilate is the gradient magnitude image after morphological dilation operation, B = strel('square', r) is a structure element, and r represents the size of B, represents the dilation operation; based on an illumination direction, the extraction of the impact crater seed point is constrained; S ini = {(x i ,y i ) | (x i ,y i ) e C ini , |G dir (x i ,y i ) - θ| < τ} where τ is defined as a preconfigured difference between the illumination direction and the gradient direction, θ is the illumination angle, S ini represents the pixel that satisfies both the gradient magnitude and direction constraints, i.e., the crater seed point.

4. The method of claim 1, wherein, The impact crater candidate region extraction based on the statistical morphology constraint path set comprises the following steps: based on graph theory, an image spatial domain and a path representation thereof are defined, wherein the image spatial domain is defined as a directed and weighted acyclic graph, the path extension directions are three directions along the illumination direction, perpendicular to the illumination direction, and 45° directions with respect to both directions, and a path is defined as a pixel sequence, which constitutes a path only when and if adjacent elements in the path satisfy a path operator construction condition; based on the current impact crater seed point as a starting point, along the illumination direction and the reverse illumination direction, pixels are processed based on statistical morphology dilation and erosion operations respectively; based on the statistical morphology processing result, the path operator construction condition is determined, the corresponding path operator is constructed, the bright features and the dark features of the impact craters are detected respectively, and the path sets of the bright features and the dark features are obtained; based on the path sets of the bright features and the dark features, the candidate region of the impact crater is determined.

5. The method according to claim 4, wherein, The processing of the pixels based on the statistical morphology dilation and erosion operations is represented as: where I'(x i ) represents the statistical morphological result of pixel x i ; I(x i ) denotes the image gray value of pixel x i ; x i r is the neighborhood pixel of pixel x i ; r = 1, 2, …, 8, I(x i r ) represents the image gray value of pixel x i r , and Z(x i ) is a normalization constant; δ r is an indicator function, that is, for the morphological dilation operation, if I(x i r ) ≥ I(x i ), δ r = 1, otherwise δ r = 0, for the morphological erosion operation, if I(x i r ) ≤ I(x i ), δ r = 1, otherwise δ r = 0.

6. The method of claim 5, wherein, The path operator construction condition is: For the detection of bright features of the impact crater, the pixel points satisfying I'(x j - I(x i ) ≥ 0 are classified into the bright feature path; For the detection of dark features of the impact crater, the pixel points satisfying I'(x j - I(x i ) < 0 are classified into the dark feature path.

7. The method of claim 1, wherein, The impact crater candidate region based on the structural consistency constraint multiple paths comprises the following steps: based on a directed and weighted graph, a frequency domain path representation is defined, wherein the path extension directions are along the illumination direction, the reverse illumination direction, and the orthogonal directions of the two directions, a path is defined as a pixel sequence, which constitutes a path only when and if adjacent elements in the path satisfy a structural consistency constraint multiple path operator construction condition; based on a log-Gabor filter, image domain information of the image is converted into frequency domain information; based on the frequency domain information, a local phase is determined; based on the local phase, the structural consistency constraint multiple path operator construction condition is determined with the current seed point as the center, the corresponding path operator is constructed, the bright features and the dark features of the impact craters are detected, and thus the structural consistency multiple path sets of the bright features and the dark features are constructed; based on the structural consistency multiple path sets, the candidate region of the impact crater is determined, a structural consistency index is calculated based on the Euclidean distance between the paths, and if the structural consistency index is less than a preconfigured threshold value, it is indicated that the candidate region covers completely.

8. The method according to claim 7, wherein, The log-Gabor filter is defined as: wherein w is a current pixel point frequency, w0 is a center frequency of the filter, and σ0 is a filter bandwidth; The frequency domain information after log-Gabor filtering is denoted as f M (x), f M (x) = [f e (x), f o (x)] T wherein f e (x) represents an odd component with asymmetric properties in the frequency domain, and f o (x) represents an even component with symmetric properties in the frequency domain.

9. The method according to claim 8, wherein, The local phase is: If the local phase satisfies the structural consistency constraint multipath operator construction condition The corresponding pixel point is classified into the bright feature path. If the local phase satisfies the structural consistency constraint multipath operator construction condition The corresponding pixel point is classified into the dark feature path.

10. The method of claim 1, wherein, The difference factor is defined as: where C i I = (v i I , u i I , r i I ) represents the impact crater candidate region based on statistical morphological constraint path set extraction, (v i I , u i I ) is the impact crater center coordinate, and r i I is the radius; C i F = (v i F , u i F , r i F ) represents the impact crater candidate region based on structural consistency constraint multi-path extraction, (v i F , u i F ) is the impact crater center coordinate, and r i F is the radius; R max and R min represent the maximum and minimum radii of C i I and C i F , respectively; d(v i (C i I ~ C i F , u i (C i I ~ C i F )) represents the distance between the impact crater centers of C i I and C i F . If the difference factor CDF is less than a preconfigured threshold, the current impact crater candidate region is considered as a correctly extracted impact crater, otherwise, the current extracted candidate region is rejected.

Images