# Robust Crater Detection Algorithm Based on Maximum Entropy Threshold Segmentation

Authors: Peng Wu; Rongjun Mu; Yanpeng Deng
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

For future lunar exploration and planetary missions, the digital elevation model (DEM) of the target object may not be well prepared before the mission, so developing a new robust crater detection algorithm (CDA) without prepared high-precision DEM is needed to meet the requirements of a high-reliability and high-precision detection and navigation system. In this paper, we presented a new robust lunar CDA method based on maximum entropy threshold segmentation. By calculating the entropy distribution of the ternary image, the threshold for retaining the maximum amount of image information is selected adaptively, a variety of evaluation indicators are proposed, and a multiple-indicator constraint matrix is constructed to realize the extraction and fitting of the craters. The proposed method has the following advantages: (1) it has strong robustness and is capable of extracting complete craters under multiple illumination conditions, which makes it suitable for the extraction of large-scale planetary and lunar images; (2) the extracted crater edges are clear and complete and do not merge with the surrounding environment edge; and (3) it avoids the problem of parameter sensitivity that is present in a traditional CDA algorithm. The proposed method was verified using an image taken by the Chang’e-2 lunar probe, and a comparison with the traditional method based on morphology and adaptive Canny edge detection shows that the number of craters detected increases by more than 35%, while the computational efficiency is improved by more than 40%.


---

## Body

## 1. Introduction

Planetary landing exploration is a key approach to study the formation and evolution of the solar system, the development and utilization of space resources, asteroid impact defense, and other major scientific issues. In the past decades, many countries have carried out a series of successful exploration missions for the Moon (Apollo [1, 2], Chang’e [3–7]), Mars (Tianwen-1 [8], Curiosity [9], and Opportunity [10]), asteroids (Rosetta [11]), and other celestial bodies. In general, terrain with high scientific value is more complex, and landing is more difficult; therefore, crater detection algorithms are important for lunar exploration and planetary missions. Chang’e-4 is the only exploration mission in a complex area that has successfully landed on the far side or back of the Moon [12, 13]. Compared with the traditional lunar sea area on the front of the Moon, the terrain on the back of the Moon is more complex and lacks high-precision terrain elevation data. Furthermore, traditional communication technology based on deep space networks presents greater risks for a mission due to communication delay problems [14], and therefore, algorithms need to be more robust.To overcome the limitation of the traditional deep space network, it is necessary to develop a robust autonomous detection method for lunar craters under any illumination and construct autonomous hazard detection and avoidance (HDA) technology based on optical information for accurate navigation and autonomous landing-point selection for the lunar surface high-precision descent landing guidance, navigation, and control (GNC) system.To detect craters, many crater detection algorithms (CDAs) have been studied by researchers all over the world. Some methods have focused on positive detection rates, while other methods pay more attention to speed or robustness. Studies by Stepinski et al., Delatte et al., Silburt et al., and Emami et al. [15–18] used machine learning and convolutional neural networks (CNNs) separately to process digital elevation models (DEMs). Studies by Wang et al. [19] presented a high generalization performance transfer learning-based CDA. Some traditional methods are also used in this field [20], such as the Hough transform, Canny edge detection [21, 22], cascade decision forest [23], solar illumination direction, and texture analysis [24]. Jia et al. extracted more all craters larger than 200 m of Chang’e-5 landing area based on digital orthophoto map generated from more than 700 Lunar Reconnaissance Orbiter Camera, Narrow Angle Camera [25]. Solarna et al. proposed a marked point process CDA method using graph cuts and decimated wavelets [26]. Bandeira et al. detected subkilometer craters and achieved a high detection rate through texture analysis and shape information [27]. Galloway et al. [28] used a Hough transform and Canny edge detection for processing high-resolution images, which required a large number of calculations. Lee and Hogan [29] proposed an automated CDA with human-level performance by multiple neural networks, but the method used not only DEM but also thermal infrared imagery, so the method did not have universal applicability. Francis et al. [30] presented a dataset of craters for CDA based on deep learning. Cui et al. [31] directly used the crater curve matching for the planetary landing navigation system and obtained superior accuracy. Chen et al. [32] presented a CDA based on terrain analysis and mathematical morphology applied to DEM with a resolution of 100 m. Lee et al. [33] presented a CNN-based CDA, but the method needs to determine detection probabilities in advance for mission planning.However, many existing CDAs have different limitations. Methods based on machine learning and CNNs require suitable train sets [34]. The Hough transform has poor adaptability because only circular craters can be detected. Furthermore, texture analysis cannot be adjusted to a variety of lighting conditions. Therefore, these methods have problems with robustness or speed; CDA based on image or optical information still have research value. In the view of CDA using unsupervised algorithms and CNNs, Emami and other researchers made a detailed research and comparative analysis [35]. From the study of Emami et al. [35], we can conclude that CDA based on convex grouping and interest points performed better in combining with CNNs for hypothesis verification. Some research results indicate that CDAs based on optical information processing still have great research values and have application prospects in the field of navigation. Maass and other researchers [36] used mainly image segmentation method to detect craters and applied research results to CNAV, which achieved good performance. Furthermore, Maass [37] firstly investigated image illumination direction in segmentation-based CDA for spacecraft navigation and then pointed that the navigation estimator can be considered robust and accurate enough in proposed CDA.To meet the needs of reliable on-orbit real-time navigation and HDA systems, a robust and fast CDA is proposed in this paper. First, the gray image is filtered by an adaptive Gaussian filter to remove image noise while retaining edge information. Then, the image is divided into highlighted, shadowed, and average gray areas via maximum entropy threshold segmentation to reduce the sensitivity of the results to light sources and to improve the reliability of Canny edge detection. Finally, a multiple-indicator constraint method is proposed to match two edges that belong to one crater, and an ellipse is used to fit the two edges.The structure of this paper is as follows. In Section2, an image ternary segmentation method based on a maximum entropy threshold is proposed to deal with the disadvantages of traditional and adaptive Canny detection in CDAs; in Section 3, a multiple-indicator constraint matrix is constructed to accomplish crater edge matching and ellipse fitting; in Section 4, simulations are presented, and a lunar image is used as an example to demonstrate the detection quality and computational cost of the proposed model; finally, our conclusions are drawn in Section 5.


## 2. Crater Edge Detection Based on Maximum Entropy Threshold

In this study, craters were detected by calculating the gradient of the image, the maximum entropy threshold was used to segment the image before Canny edge detection, and then the craters could be detected after the clear edges were extracted. Compared to traditional methods, the crater edges extracted by our method are clear and obvious and do not merge with the surrounding environment; they do not need morphological processing such as dilation and erosion and can be directly used for crater edge matching and ellipse fitting.


### 2.1. Adaptive Canny Edge Detection

The craters on the lunar surface form highlighted and shadowed areas that are quite different from the surrounding environment and form obvious gray gradient edges under direct sunlight. Therefore, craters can be effectively detected by edge detection. There are many methods of edge detection, such as the Roberts operator, Prewitt operator, and Canny operator, but the Roberts and Prewitt operators have two disadvantages. The first disadvantage is that they are easily affected by noise and can detect multiple edges; the other disadvantage is that the scale of these operators is fixed, which is not conducive to detecting edges of different scales. Compared with the Roberts and Prewitt operators, the Canny edge detection can effectively suppress multiresponse edges and can improve the accuracy of edge detection via nonmaximum suppression. Therefore, in this study, we used the Canny edge detection.First, for denoising the image and preserving the edge information of craters, a Gaussian filter is commonly used. The Sobel operator is then used to calculate the gradient of the image as follows:(1)Sx=−101−202−101Sy=121000−1−2−1,where Sx,Sy are Sobel operators along the x-axis and y-axis, respectively. The norm G and direction θ of the image gradient can be calculated as follows:
(2)G=Gx2+Gy2=Sx∗A2+Sy∗A2,θ=arctanSx∗A/Sy∗A,where A is the convolution window of the Canny edge detection. Finally, the strong and weak edges are extracted by nonmaximum suppression and dual-threshold detection. Dual-threshold, high-threshold τh, and low-threshold τl are the extraction standards of crater edges. The edges with gradient intensity bigger than high-threshold τh are sure to be crater edges, called strong crater edges, and the edges with gradient intensity smaller than low-threshold τl are sure to be noncrater edges. The edges whose gradient intensity lies between high-threshold τh and low-threshold τl are classified as weak crater edges or nonedges by their connectivity. As shown in Figure 1, the gradient intensity of edge A is bigger than τh, so edge A is sure to be a crater edge. The gradient intensity of crater B is between τh and τh, and edge B is connected to edge A, so it is a weak crater edge. Edge C is not connected to any crater edges, so it is not a crater edge.Figure 1
Crater edge extraction method by the Canny edge detection.It can be concluded that the selection of dual-thresholds, high-thresholdτh, and low-threshold τl will seriously affect the effect of the Canny edge detection. Various values of τh and τl were selected to analyze a 7 m resolution digital orthophoto map (DOM) image of the lunar surface taken by Chang’e-2. The influence of different τh and τl values is shown in Figure 2. Figure 2(a) shows the original lunar image including regions A and B; Figure 2(b) shows the original region A image and different Canny edge detection results for region A with different τl; Figure 2(c) shows the original region B image and different Canny edge detection results for region B with different τh; Figure 2(d) shows different Canny edge detection results for region B with different τl.Figure 2
Performance and results for various values of high-thresholdτh and low-threshold τl in the Canny edge detection.
(a)(b)(c)(d)We can infer from the results that the crater detection algorithm based on Canny has the following two disadvantages due to the parameter sensitivity: (i) the edges detected by Canny usually connect with surrounding edges, as shown in Figure2(b); (ii) it is hard to choose an appropriate dual-threshold autonomously under different light conditions, which cannot meet the demand of a lunar exploration mission. To address these problems, He [38] proposed an adaptive dual-threshold selection method in the Canny edge detection algorithm, i.e., through the histogram statistics of the gray image after nonmaximum suppression, determine the gradient corresponding to the maximum number of pixels gm and the gradient amplitude mean square deviation of the maximum number of pixels σm as follows:
(3)σm=∑i=1khgigi−gm2∑i=1khgi,where k is the number of histogram bins and hgi is the number of pixels that fall into bin whose gray value is gi.The high-thresholdτh is as follows:
(4)τh=gm+σm.Remove the pixels whose gray value is stronger thangm, and recalculate gm′ and σm′ using the steps above. The low-threshold τl is as follows:
(5)τl=gm′+σm′.This adaptive dual-threshold selection method was used to detect crater edges, as shown in Figure2(a), and the results are given in Figure 3. For the same shooting area, by selecting different regions to analyze, the dual-thresholds selected by the traditional Canny edge detection method are quite different, which were easily affected by the camera field of view and the shooting area.Figure 3
Performance and results of the adaptive dual-threshold selection method.
(a)(b)First, dual-thresholdsτl,τh were chosen to extract crater edges, when only region A or region B is analyzed, and the corresponding τl,τh values are different from the τl,τh of the whole image, as shown in Table 1, which leads to different results; this is not stable and robust, as shown in Figures 3(a) and 3(b).Table 1
Dual-thresholdsτl,τh of different images.
ImageτlτhRegion A0.2780.774Region B0.1180.728Whole mage0.3130.815Figure3(a) shows the Canny edge detection results for the whole image; Figure 3(b) shows the Canny edge detection results for regions A and B, where (i) is the enlargement of region A in Figure 3(a), (ii) is the Canny edge detection result for region A alone, (iii) is the enlargement of region B in Figure 3(a), and (iv) is the Canny edge detection result for region B alone.From Table1 and Figure 3 results, we can infer that although the adaptive dual-threshold selection method in the Canny edge detection algorithm by Jiang [38] can select dual-threshold autonomously and extract the crater edges, the edges are usually broken and discontinuous and are sometimes connected with the surrounding area, which is not conducive to further processing. In those circumstances, the edges of craters cannot be extracted effectively for ellipse fitting; therefore, an image ternary segmentation method based on the maximum entropy threshold is proposed, which can effectively solve this problem.


### 2.2. Image Ternary Segmentation Method Based on Maximum Entropy Threshold

Under sunlight, the shadowed area, the highlighted area, and the surrounding environment can be seen on the lunar surface. The ternary segmentation method based on maximum entropy threshold is proposed to distinguish these three parts. After obtaining a filtered image, examining all gray colors, and calculating image entropy, the image can be ternary segmented by two thresholds,α,βα≤β, which correspond to maximum entropy.The entropy of the imageH is
(6)H=−∑i=0255pilogpi,where pi is the probability density of the pixels whose gray value is i.We selectα,β to segment the image and divide the image into shadowed areas C1, background areas C2, and highlighted areas C3. The corresponding cumulative probability is as follows:
(7)P1α=∑i=0αpi,P2α,β=∑i=α+1βpi,P3β=∑i=β+1255pi.Pi is satisfied as follows:
(8)P1α+P2α,β+P3β=1.The entropies of those three areas are as follows:(9)H1α=−∑i=0αpiP1αlogpiP1α,,H2α,β=−∑i=α+1βpiP2α,βlogpiP2α,β,H3β=−∑i=β+1255piP3βlogpiP3β.The total entropy of the image is the sum ofHi as follows:
(10)Hα,β=H1α+H2α,β+H3β.We obtain a maximumHα,β by keeping the image information via selecting appropriate values of α,β.
(11)α,β=argmaxα,β∈0,255Hα,β.By choosing differentα,β values to segment the image, as shown in Figure 4(a), this leads to the distribution of the segmented image entropy, as shown in Figure 5. By finding the thresholds α,β corresponding to the maximum entropy, the ternary segmentation image result based on the maximum entropy threshold is as shown in Figure 4(b).Figure 4
Ternary segmentation image result.
(a)(b)Figure 5
Distribution of segmented image entropy.We then detect the edges of craters after getting a ternary segmentation image. For the two different thresholdsτl,τh, the Canny edge detection produces results similar to those shown in Figure 6. From the results, we can infer that this method is robust in terms of dual-thresholds τl,τh. No matter how the dual-thresholds τl,τh are chosen, clear and continuous edges can be detected after image ternary segmentation, which helps to match semicrater edges and ellipse fitting in the next step.Figure 6
Canny detection results with various values ofτl,τh.


## 2.1. Adaptive Canny Edge Detection

The craters on the lunar surface form highlighted and shadowed areas that are quite different from the surrounding environment and form obvious gray gradient edges under direct sunlight. Therefore, craters can be effectively detected by edge detection. There are many methods of edge detection, such as the Roberts operator, Prewitt operator, and Canny operator, but the Roberts and Prewitt operators have two disadvantages. The first disadvantage is that they are easily affected by noise and can detect multiple edges; the other disadvantage is that the scale of these operators is fixed, which is not conducive to detecting edges of different scales. Compared with the Roberts and Prewitt operators, the Canny edge detection can effectively suppress multiresponse edges and can improve the accuracy of edge detection via nonmaximum suppression. Therefore, in this study, we used the Canny edge detection.First, for denoising the image and preserving the edge information of craters, a Gaussian filter is commonly used. The Sobel operator is then used to calculate the gradient of the image as follows:(1)Sx=−101−202−101Sy=121000−1−2−1,where Sx,Sy are Sobel operators along the x-axis and y-axis, respectively. The norm G and direction θ of the image gradient can be calculated as follows:
(2)G=Gx2+Gy2=Sx∗A2+Sy∗A2,θ=arctanSx∗A/Sy∗A,where A is the convolution window of the Canny edge detection. Finally, the strong and weak edges are extracted by nonmaximum suppression and dual-threshold detection. Dual-threshold, high-threshold τh, and low-threshold τl are the extraction standards of crater edges. The edges with gradient intensity bigger than high-threshold τh are sure to be crater edges, called strong crater edges, and the edges with gradient intensity smaller than low-threshold τl are sure to be noncrater edges. The edges whose gradient intensity lies between high-threshold τh and low-threshold τl are classified as weak crater edges or nonedges by their connectivity. As shown in Figure 1, the gradient intensity of edge A is bigger than τh, so edge A is sure to be a crater edge. The gradient intensity of crater B is between τh and τh, and edge B is connected to edge A, so it is a weak crater edge. Edge C is not connected to any crater edges, so it is not a crater edge.Figure 1
Crater edge extraction method by the Canny edge detection.It can be concluded that the selection of dual-thresholds, high-thresholdτh, and low-threshold τl will seriously affect the effect of the Canny edge detection. Various values of τh and τl were selected to analyze a 7 m resolution digital orthophoto map (DOM) image of the lunar surface taken by Chang’e-2. The influence of different τh and τl values is shown in Figure 2. Figure 2(a) shows the original lunar image including regions A and B; Figure 2(b) shows the original region A image and different Canny edge detection results for region A with different τl; Figure 2(c) shows the original region B image and different Canny edge detection results for region B with different τh; Figure 2(d) shows different Canny edge detection results for region B with different τl.Figure 2
Performance and results for various values of high-thresholdτh and low-threshold τl in the Canny edge detection.
(a)(b)(c)(d)We can infer from the results that the crater detection algorithm based on Canny has the following two disadvantages due to the parameter sensitivity: (i) the edges detected by Canny usually connect with surrounding edges, as shown in Figure2(b); (ii) it is hard to choose an appropriate dual-threshold autonomously under different light conditions, which cannot meet the demand of a lunar exploration mission. To address these problems, He [38] proposed an adaptive dual-threshold selection method in the Canny edge detection algorithm, i.e., through the histogram statistics of the gray image after nonmaximum suppression, determine the gradient corresponding to the maximum number of pixels gm and the gradient amplitude mean square deviation of the maximum number of pixels σm as follows:
(3)σm=∑i=1khgigi−gm2∑i=1khgi,where k is the number of histogram bins and hgi is the number of pixels that fall into bin whose gray value is gi.The high-thresholdτh is as follows:
(4)τh=gm+σm.Remove the pixels whose gray value is stronger thangm, and recalculate gm′ and σm′ using the steps above. The low-threshold τl is as follows:
(5)τl=gm′+σm′.This adaptive dual-threshold selection method was used to detect crater edges, as shown in Figure2(a), and the results are given in Figure 3. For the same shooting area, by selecting different regions to analyze, the dual-thresholds selected by the traditional Canny edge detection method are quite different, which were easily affected by the camera field of view and the shooting area.Figure 3
Performance and results of the adaptive dual-threshold selection method.
(a)(b)First, dual-thresholdsτl,τh were chosen to extract crater edges, when only region A or region B is analyzed, and the corresponding τl,τh values are different from the τl,τh of the whole image, as shown in Table 1, which leads to different results; this is not stable and robust, as shown in Figures 3(a) and 3(b).Table 1
Dual-thresholdsτl,τh of different images.
ImageτlτhRegion A0.2780.774Region B0.1180.728Whole mage0.3130.815Figure3(a) shows the Canny edge detection results for the whole image; Figure 3(b) shows the Canny edge detection results for regions A and B, where (i) is the enlargement of region A in Figure 3(a), (ii) is the Canny edge detection result for region A alone, (iii) is the enlargement of region B in Figure 3(a), and (iv) is the Canny edge detection result for region B alone.From Table1 and Figure 3 results, we can infer that although the adaptive dual-threshold selection method in the Canny edge detection algorithm by Jiang [38] can select dual-threshold autonomously and extract the crater edges, the edges are usually broken and discontinuous and are sometimes connected with the surrounding area, which is not conducive to further processing. In those circumstances, the edges of craters cannot be extracted effectively for ellipse fitting; therefore, an image ternary segmentation method based on the maximum entropy threshold is proposed, which can effectively solve this problem.


## 2.2. Image Ternary Segmentation Method Based on Maximum Entropy Threshold

Under sunlight, the shadowed area, the highlighted area, and the surrounding environment can be seen on the lunar surface. The ternary segmentation method based on maximum entropy threshold is proposed to distinguish these three parts. After obtaining a filtered image, examining all gray colors, and calculating image entropy, the image can be ternary segmented by two thresholds,α,βα≤β, which correspond to maximum entropy.The entropy of the imageH is
(6)H=−∑i=0255pilogpi,where pi is the probability density of the pixels whose gray value is i.We selectα,β to segment the image and divide the image into shadowed areas C1, background areas C2, and highlighted areas C3. The corresponding cumulative probability is as follows:
(7)P1α=∑i=0αpi,P2α,β=∑i=α+1βpi,P3β=∑i=β+1255pi.Pi is satisfied as follows:
(8)P1α+P2α,β+P3β=1.The entropies of those three areas are as follows:(9)H1α=−∑i=0αpiP1αlogpiP1α,,H2α,β=−∑i=α+1βpiP2α,βlogpiP2α,β,H3β=−∑i=β+1255piP3βlogpiP3β.The total entropy of the image is the sum ofHi as follows:
(10)Hα,β=H1α+H2α,β+H3β.We obtain a maximumHα,β by keeping the image information via selecting appropriate values of α,β.
(11)α,β=argmaxα,β∈0,255Hα,β.By choosing differentα,β values to segment the image, as shown in Figure 4(a), this leads to the distribution of the segmented image entropy, as shown in Figure 5. By finding the thresholds α,β corresponding to the maximum entropy, the ternary segmentation image result based on the maximum entropy threshold is as shown in Figure 4(b).Figure 4
Ternary segmentation image result.
(a)(b)Figure 5
Distribution of segmented image entropy.We then detect the edges of craters after getting a ternary segmentation image. For the two different thresholdsτl,τh, the Canny edge detection produces results similar to those shown in Figure 6. From the results, we can infer that this method is robust in terms of dual-thresholds τl,τh. No matter how the dual-thresholds τl,τh are chosen, clear and continuous edges can be detected after image ternary segmentation, which helps to match semicrater edges and ellipse fitting in the next step.Figure 6
Canny detection results with various values ofτl,τh.


## 3. Crater Edge Matching and Ellipse Fitting Based on Multiple-Indicator Constraints

After image ternary segmentation and crater edge detection, we need to match pairs of edges that belong to the same crater. To achieve this goal, we first use a series of feature points to represent the semicrater edge, as Figure7 shows.Figure 7
Feature points of semicrater edge.One semicrater edge can be represented by four feature points; they are the verticesA and B, geometric center C1, and midpoint of vertices C2. Then, the orientation vector of the semicrater edge r can be defined as follows:
(12)r=C1C2⟶.Taking into account the length, distance, orientation, brightness, and other constraints comprehensively, we match two semicrater edgesv1,v2 by constructing multiple constraints Ji. When Ji=1, two edges match best, and when Ji=0, two edges match worst; they do not match at all. The constraints are constructed as follows:
(1)
Normalized length similarity constraintJ1vi,vjThe lengths of two matched edges should not be widely different. The bigger the constraint, the better the two semicrater edges match.(13)J1vi,vj=1−li−ljli+ljJ1∈0,1.where li,lj are pixel lengths of the edges to be matched
(2)
Brightness constraintJ2vi,vjUnder sunlight, two semicrater edges belonging to the same crater should be in different highlighted and shadowed areas. The gray value of the geometric centerC1 is used to distinguish the two semicrater edges to be matched as follows:
(14)J2vi,vj=1lexi≠lexj0lexi=lexj,where lexi,lexj are the two geometric centers’ gray values of the edges to be matched. lexi≠lexj indicates that two edges belong to the highlighted area and the shadowed area, respectively
(3)
Normalized orientation constraintJ3vi,vjTwo semicrater edges belonging to the same crater should be located in opposite directions, which means two orientation vectorsri,rj are at an obtuse angle. The bigger the constraint is, the better the two semicrater edges match.
(15)J3vi,vj=rirj−ri⋅rj2rirjJ3∈0,1(4)
Distance constraintJ4vi,vj,J5vi,vjTwo semicrater edges belonging to the same crater should not be far apart. We use the geometric meanlilj to measure the distance between two semicrater edges. The distance is divided into horizontal distance J4vi,vj and vertical distance J5vi,vj as follows:
(16)J4vi,vj=1r1⋅s<η1lilj0r1⋅s≥η1lilj,J5vi,vj=1r1×s<η2lilj0r1×s≥η2lilj,where s is the sunlight vector and η1 and η2 are scale parameters. In general, the distance between two semicrater edges belonging to the same crater is shorter than half the perimeter of the crater. And the geometric mean lilj is usually shorter than half the perimeter of the crater. The relationship of ηi and lilj can be written as follows:
(17)ηkliljmaxli,lj=1k=1,2.In this paper, we choseη1=η2=2.5(5)
Normalized shape constraintJ6vi,vjAccording to the above constraintsJii=1,2,⋯,5, to match the edges of craters, there are two matching results, that is, X-shape results and O-shape results; O-shape results are real craters, and X-shape results usually come from stone shadow influence. To avoid this influence, the vector from the center of the shadowed area to the center of the highlighted area should be consistent with the direction of the sunlight vector as follows:
(18)J6vi,vj=C2iC2j⟶−lexi−lexjC2iC2j⟶⋅s2C2iC2j⟶,J6∈0,1.For any two semicrater edges, the multiple-indicator constraintJvi,vj can be calculated to represent the similarity between them as follows:
(19)Jvi,vj=Ji,j=∏kJkvi,vj,k=1,2,⋯,6.Forn semicrater edges, we construct a matching similarity multiple-indicator constraint matrix as follows:
(20)D=Jv1,v1Jv1,v2⋯Jv1,vnJv2,v1Jv2,v2⋯Jv2,vn⋮⋮⋱⋮Jvn,v1Jvn,v2⋯Jvn,vn.D is a symmetric matrix with diagonal elements that are zero; Ji,j means the similarity between edges vi and vj.First, we choose two largest elementsJmax1 in D, which are located in coordinates i,j and j,i. Then, we remove the other elements in D whose positions are row i and column j, and row j and column i. We name the new multiple-indicator constraint matrix D1. Then, we choose the new two largest elements Jmax2 in D1 and repeat the above steps until the largest element Jmaxm∈Dm−1 is lower than the matching threshold ε. After the semicrater edge matching is complete, we ellipse fit pairs of edges that belong to one crater.Craters can be described mathematically as follows:(21)FX=XTAX+bTX+c=0.We use the least-squares method (LSM) to calculateA,b,c. Since the projections of the craters are ellipse, it is necessary to exclude the case for which A has nonpositive eigenvalues. After ellipse fitting, the candidate for the crater center is the following:
(22)Xo=xoyoT=−12A−1b.


## 4. Results and Discussion

Crater detection is basic to optical navigation, crater reconstruction, and HDA. To verify the effectiveness and robustness of the method proposed, two experiments are carried out in this article. The first experiment verifies that the proposed CDA works well and has robustness under different illumination directions for the same area. The second experiment verifies that compared with the traditional method, the proposed CDA can detect more craters and has lower computational efficiency.


### 4.1. Simulation Verification on Sand Table under Different Illumination Directions

Crater detection is basic to optical navigation, crater reconstruction, HAD, and landing site selection. A2m×2m simulated lunar sand table is constructed to verify the robustness of the proposed CDA. The parallel light source is set at different positions of the sand table to simulate different sunlight illumination directions. The environment around the sand table is processed to avoid the influence on the detection of sand table craters. The experimental conditions are shown in Table 2.Table 2
Experience conditions.
ItemParameterSand table size2m×2mNumber of craters8Maximum depth of crater0.30 mReference height of sand table0.18 mParallel light source5600K×2MT9V034 camera45 fps@1280×960, 60 fps@1280×72080 fps@640×480, 160 fps@320×240Hardware-in-loop simulation experiment is constructed by a front-facing camera, sand table, parallel light source, and data switching exchanges. Figure8 shows the camera and the lunar sand box used in the experiment. Figure 9 shows the hardware-in-loop simulation environment.Figure 8
Front-facing camera and lunar sand table in experiment. (a) MT9V034 camera. (b) Lunar sand table.
(a)(b)Figure 9
Hardware-in-loop simulation environment.The parallel light sources are, respectively, set on the left and right sides of the lunar simulation sand table, with a height of about one meter, to simulate the different light directions of the Sun. After Gaussian smoothing of the sand table images under different light conditions, maximum entropy threshold segmentation is used to form the different highlight and shadow areas of the crater, and then the craters are detected. Detected results are as shown in Figure10.Figure 10
Detection results under different illumination directions.
(a)(b)(c)(d)(e)(f)Figure10 shows the following: (a) parallel light source is set on the left of sand table; (b) maximum entropy threshold segmentation of (a); (c) detected craters of (a); (d) parallel light source is set on the right of sand table; (e) maximum entropy threshold segmentation of (d); and (e) detected craters of (d).We can make some discussions from Figure10.
(1)
Seven craters are detected in each image. Crater No. 1~6 are all detected in two images because their shadows are easily segmented from the surrounding environment. Crater No. 7 and crater No. 8 are all covered by the shadow of crater No. 2, so the proposed CDA can detect only one of them in each direction of sunlight. This is caused by the shape and distribution of the craters on the sand table, not by our proposed method(2)
Although crater No. 1~6 are all detected in two images, the ellipse shapes of the craters are different. That is because shadows of the same craters under different sunlight directions are usually quite different as shown in Figures10(b) and 10(e)(3)
CDA proposed can work well under different illumination directions, which means our CDA method has strong robustness.


### 4.2. Simulation Verification in Lunar Image

According to the lunar surface image data released by Chang’e series of probes, we can obtain optical images of the lunar surface with different resolutions. In this paper, a DOM (7 m) gray image, detected by the Chang’e-2 lunar probe, namedCE2_GRAS_DOM_07m_C101_66N172W_A.tif was selected as the experimental data to extract large and small craters. We can download it from https://moon.bao.ac.cn/searchOrder_dataSearchData.search. The image information and PC parameters are shown in Table 3.Table 3
Image information and PC parameters.
Image informationPC parameterImage sizeImage resolutionCPURAM9638×124727 mIntel i9-1090032 GThere were 1337 craters counted in Chang’e-2 DOM (7 m) image, and the number 1337 was used as the ground truth (Table4). To detect craters fast and accurately, a sliding window sized 963×1247 is used to detect small craters.Table 4
Ground truth and detection results.
ItemTotal numberNumber of cratersD<70m70m<D<700mD>700mGround truth133796120437Detection results117077105637Percentage87.5%80.2%87.7%100%At a matching threshold ofε=0.5, 1170 craters can be detected in this image without false positives. The detection results 1170 craters are shown in Figure 11. The image resolution is 7 m per pixel, and we can calculate the average diameter of each detected crater. For different sizes of crater, our method has different extraction rates; the bigger the crater, the higher the extraction rate. Through analysis, the main reasons for the decline in the extraction rate of smaller craters are as follows: the slope of lunar surface with small craters is too small, the changes of light and darkness are not obvious, and there is a large degree of overlap with the surrounding craters.Figure 11
Crater detection results.To make the detection results more obvious, Figure12 shows the ground truth of all craters by manual counting. The red circles represent the center of the undetected craters, and the blue circles represent the crater detected by the method proposed.Figure 12
Ground truth including detected craters and undetected craters.Table5 shows the extracted information for craters whose average diameters are larger than 700 m. The recorded information includes the average diameter, center coordinates, and direction angle (θ).Table 5
D>700m crater extracted information.
Crater numberAverage diameter (m)Center coordinatesθ (rad)Crater numberAverage diameter (m)Center coordinatesθ (rad)1743.8(20593, 2888)−3.07222847.3(3258, 4001)3.105201036.8(10940, 8802)−3.1403784.0(233, 3486)2.620212678.5(660, 3470)−3.0794831.8(621, 5962)2.630224977.8(1060, 7760)2.6755983.9(833, 5802)2.445232922.7(2450, 1900)−1.33161070.3(919, 6135)2.773241321.5(2910, 5960)−2.9817784.1(10029, 6139)3.131257467.4(5240, 2290)−2.91981023.9(2016, 7452)3.046263782.9(5440, 5050)−0.9409966.5(2480, 7419)2.749273986.1(6270, 7120)−2.936103875.1(4856, 6660)−2.461282407.3(6160, 1370)2.538111137.6(7525, 6770)−1.293293973.4(6880, 1480)−2.424122342.1(9993, 6947)2.566308360.0(8720, 1000)−1.814132101.2(11755, 7191)−1.003311153.2(9230, 3740)−1.76914905.5(1063, 9242)−1.056325564.8(6040, 3940)−1.483151187.7(1699, 8675)2.774333521.4(11610, 940)−2.494162828.8(3555, 9513)2.802341777.8(4210, 5620)2.916171652.6(9167, 9006)−2.298353665.2(11120, 8600)2.735181102.9(9445, 8651)2.832362281.0(11210, 6670)−0.908192271.4(10194, 9394)−0.792376030.9(10430, 5980)3.014Many researchers [39, 40] proposed different parameters to evaluate the performance of the CDA. These methods work by counting craters detected in different states, i.e., true positives (TPs) corresponding to real detected craters and false positives (FPs) corresponding to detected craters that do not exist. A false negative (FN) represents the number of existing craters that are not detected. After obtaining the TP, FP, and FN, the quality factors can be calculated, which include the detection percentage (D), the quality percentage (Q), the branching factor (B), and the precision (P) as follows:
(23)Detectionpercentage=100%×TPTP+FN,Qualitypercentage=100%×TPTP+FP+FN,Branchingfactor=FPTP,Precision=100%×TPTP+FP.P shows that our method proposed is better than random guess; D, Q, and B can help us to choose the best threshold.Different thresholds were chosen, i.e.,ε from 0.1 to 0.9, and then the detected crater numbers and quality factors were calculated using our proposed method. The results are shown in Table 6. The symbol ∗ represents the optimal value.Table 6
Detected crater numbers and quality factors.
Matching thresholdεNumber of cratersQuality factorsDetection resultsDQBPTPFPFN0.1143089.8%∗76.7%0.1984.0%12012291360.2127689.7%84.8%0.0694.0%1199771380.3126489.7%85.5%0.0594.9%1199651380.4123088.6%85.6%0.0496.3%1184461530.5117087.5%87.5%∗0∗100%∗117001670.6107480.3%80.3%0∗100%∗107402630.796272.0%72.0%0∗100%∗96203750.890968.0%68.0%0∗100%∗90904280.920915.6%15.6%0∗100%∗20901128The detection percentage curve as shown in Figure13 was drawn using the data in Table 6.Figure 13
Detection percentage (D) curve.We can make the following conclusions from Table6 and Figure 13.
(1)
Whenε=0.1, although the number of detected craters is even more than manual counting ground truth, TP and detection percentage do not reach 100%. There are two reasons, the first reason is that some structures of old craters are broken, and the other reason is that some semicrater edges mismatch to other craters(2)
When the matching thresholdε increases, the constraint for matching semicrater edges becomes stronger, and therefore, the number of extracted craters decreases; TP and FP decrease at the same time(3)
When the matching thresholdε increases, FN increases along with it, because of irregularly shaped, incomplete crater edges, or craters that did not have an obvious change in gray values(4)
The CDA achieves the best result only by choosing the appropriate threshold. When the matching thresholdε increases, quality factors D and B become better but the number of detected craters decreases. When ε=0.5, there is no FP which matters with navigation, and Q and B are all the best value, which means ε=0.5 is an appropriate threshold(5)
Design an evaluation parameterE to comprehensively consider various quality factors.(24)E=D+Q+1−B+P.Whenε=0.5, E=3.75, which is larger than other parameters that correspond to other ε.The area under the curve (AUC) in Figure13 could be calculated to evaluate the performance of the proposed method; the higher the AUC, the better the proposed method.
(i)
AUC=1 means that the proposed CDA is perfect, and it can detect all craters without any FP and FN(ii)
0.5<AUC<1 means the proposed CDA is better than random guess. If the matching threshold ε is set appropriately, the proposed method can work very well(iii)
AUC=0.5 means that the proposed method works the same with random guess(iv)
AUC<0.5 means that the proposed method works even worse than random guess.The AUC of the proposed method is 0.9192, which is much higher than 0.5, which means that our method works well and has robustness for crater detection.Table7 shows a comparison of the quality factors of our proposed method with He’s CDA [38]. For the same lunar surface picture, our method obtained a better detection percentage, quality percentage, and branching factor, due to a larger number of TPs, smaller number of FPs, and smaller number of FNs.Table 7
Comparison with other crater detection algorithms.
Method (ε=0.5)Crater numberQuality factorsDetection resultsDQBTPFPFNOur approach117087.5%87.5%011700167He [38]97964.6%60.0%0.12865104472Hough transform91745.8%38.5%0.42613255724Highlight/shadow region104759.2%52.5%0.22792171545Compared to He’s method and other traditional methods, the proposed method improves TP by 35.3%, and FP and FN are much smaller than He’s method. Quality factors are also better. It can be concluded that the proposed method is a more effective large-scale image CDA.In Table8, we compared the quality factors using a published lunar crater database [24, 41]; there are 4333 craters at the landing site.Table 8
Comparison with He’s CDA.
Methods (ε=0.5)Crater numberQuality factorsDetection resultsDQBTPFPFNOur approach117088.5%88.1%0.005383519498He [38]97958.4%44.2%0.55253214011801By dividing the crater detection algorithm into substeps (image preprocessing, ternary segmentation, adaptive dual-thresholdτlτh selection, edge detection, edge matching, and ellipse fitting), we can compare the difference between the proposed method and He’s method. Therefore, we selected a subgraph in Figure 11 of 963×1247, calculated the time for each step, and analyzed the computational cost; the result is shown in Figure 14 and Table 9.Figure 14
Computational cost of different methods.Table 9
Computational cost comparison with He’s CDA.
Method (ε=0.5)Image preprocessingTernary segmentationτl,τh selectionEdge detectionEdge matchEllipse fittingTotal timeOur approach0.395 s0.257 s0 s0.353 s0.701 s0.025 s1.831 sHe’s approach0.364 s0 s1.112 s0.664 s0.842 s0.015 s3.057 sCompared with He’s method, the computational cost of the proposed CDA is significantly reduced. The time for image preprocessing is basically the same; He’s method takes a long time to adaptively select two thresholds, while ternary segmentation based on the maximum entropy threshold, as proposed in this paper, takes a short time. In addition, as compared with He’s method, using our method, the edges detected are continuous, smooth, and without broken edges, and the time for edge matching is also shorter. Because our method has less mismatching and strong robustness, the ellipse fitting time is also slightly less than for He’s method. Overall, compared to He’s method, the computational efficiency of He’s method can be improved by 40.1%.


## 4.1. Simulation Verification on Sand Table under Different Illumination Directions

Crater detection is basic to optical navigation, crater reconstruction, HAD, and landing site selection. A2m×2m simulated lunar sand table is constructed to verify the robustness of the proposed CDA. The parallel light source is set at different positions of the sand table to simulate different sunlight illumination directions. The environment around the sand table is processed to avoid the influence on the detection of sand table craters. The experimental conditions are shown in Table 2.Table 2
Experience conditions.
ItemParameterSand table size2m×2mNumber of craters8Maximum depth of crater0.30 mReference height of sand table0.18 mParallel light source5600K×2MT9V034 camera45 fps@1280×960, 60 fps@1280×72080 fps@640×480, 160 fps@320×240Hardware-in-loop simulation experiment is constructed by a front-facing camera, sand table, parallel light source, and data switching exchanges. Figure8 shows the camera and the lunar sand box used in the experiment. Figure 9 shows the hardware-in-loop simulation environment.Figure 8
Front-facing camera and lunar sand table in experiment. (a) MT9V034 camera. (b) Lunar sand table.
(a)(b)Figure 9
Hardware-in-loop simulation environment.The parallel light sources are, respectively, set on the left and right sides of the lunar simulation sand table, with a height of about one meter, to simulate the different light directions of the Sun. After Gaussian smoothing of the sand table images under different light conditions, maximum entropy threshold segmentation is used to form the different highlight and shadow areas of the crater, and then the craters are detected. Detected results are as shown in Figure10.Figure 10
Detection results under different illumination directions.
(a)(b)(c)(d)(e)(f)Figure10 shows the following: (a) parallel light source is set on the left of sand table; (b) maximum entropy threshold segmentation of (a); (c) detected craters of (a); (d) parallel light source is set on the right of sand table; (e) maximum entropy threshold segmentation of (d); and (e) detected craters of (d).We can make some discussions from Figure10.
(1)
Seven craters are detected in each image. Crater No. 1~6 are all detected in two images because their shadows are easily segmented from the surrounding environment. Crater No. 7 and crater No. 8 are all covered by the shadow of crater No. 2, so the proposed CDA can detect only one of them in each direction of sunlight. This is caused by the shape and distribution of the craters on the sand table, not by our proposed method(2)
Although crater No. 1~6 are all detected in two images, the ellipse shapes of the craters are different. That is because shadows of the same craters under different sunlight directions are usually quite different as shown in Figures10(b) and 10(e)(3)
CDA proposed can work well under different illumination directions, which means our CDA method has strong robustness.


## 4.2. Simulation Verification in Lunar Image

According to the lunar surface image data released by Chang’e series of probes, we can obtain optical images of the lunar surface with different resolutions. In this paper, a DOM (7 m) gray image, detected by the Chang’e-2 lunar probe, namedCE2_GRAS_DOM_07m_C101_66N172W_A.tif was selected as the experimental data to extract large and small craters. We can download it from https://moon.bao.ac.cn/searchOrder_dataSearchData.search. The image information and PC parameters are shown in Table 3.Table 3
Image information and PC parameters.
Image informationPC parameterImage sizeImage resolutionCPURAM9638×124727 mIntel i9-1090032 GThere were 1337 craters counted in Chang’e-2 DOM (7 m) image, and the number 1337 was used as the ground truth (Table4). To detect craters fast and accurately, a sliding window sized 963×1247 is used to detect small craters.Table 4
Ground truth and detection results.
ItemTotal numberNumber of cratersD<70m70m<D<700mD>700mGround truth133796120437Detection results117077105637Percentage87.5%80.2%87.7%100%At a matching threshold ofε=0.5, 1170 craters can be detected in this image without false positives. The detection results 1170 craters are shown in Figure 11. The image resolution is 7 m per pixel, and we can calculate the average diameter of each detected crater. For different sizes of crater, our method has different extraction rates; the bigger the crater, the higher the extraction rate. Through analysis, the main reasons for the decline in the extraction rate of smaller craters are as follows: the slope of lunar surface with small craters is too small, the changes of light and darkness are not obvious, and there is a large degree of overlap with the surrounding craters.Figure 11
Crater detection results.To make the detection results more obvious, Figure12 shows the ground truth of all craters by manual counting. The red circles represent the center of the undetected craters, and the blue circles represent the crater detected by the method proposed.Figure 12
Ground truth including detected craters and undetected craters.Table5 shows the extracted information for craters whose average diameters are larger than 700 m. The recorded information includes the average diameter, center coordinates, and direction angle (θ).Table 5
D>700m crater extracted information.
Crater numberAverage diameter (m)Center coordinatesθ (rad)Crater numberAverage diameter (m)Center coordinatesθ (rad)1743.8(20593, 2888)−3.07222847.3(3258, 4001)3.105201036.8(10940, 8802)−3.1403784.0(233, 3486)2.620212678.5(660, 3470)−3.0794831.8(621, 5962)2.630224977.8(1060, 7760)2.6755983.9(833, 5802)2.445232922.7(2450, 1900)−1.33161070.3(919, 6135)2.773241321.5(2910, 5960)−2.9817784.1(10029, 6139)3.131257467.4(5240, 2290)−2.91981023.9(2016, 7452)3.046263782.9(5440, 5050)−0.9409966.5(2480, 7419)2.749273986.1(6270, 7120)−2.936103875.1(4856, 6660)−2.461282407.3(6160, 1370)2.538111137.6(7525, 6770)−1.293293973.4(6880, 1480)−2.424122342.1(9993, 6947)2.566308360.0(8720, 1000)−1.814132101.2(11755, 7191)−1.003311153.2(9230, 3740)−1.76914905.5(1063, 9242)−1.056325564.8(6040, 3940)−1.483151187.7(1699, 8675)2.774333521.4(11610, 940)−2.494162828.8(3555, 9513)2.802341777.8(4210, 5620)2.916171652.6(9167, 9006)−2.298353665.2(11120, 8600)2.735181102.9(9445, 8651)2.832362281.0(11210, 6670)−0.908192271.4(10194, 9394)−0.792376030.9(10430, 5980)3.014Many researchers [39, 40] proposed different parameters to evaluate the performance of the CDA. These methods work by counting craters detected in different states, i.e., true positives (TPs) corresponding to real detected craters and false positives (FPs) corresponding to detected craters that do not exist. A false negative (FN) represents the number of existing craters that are not detected. After obtaining the TP, FP, and FN, the quality factors can be calculated, which include the detection percentage (D), the quality percentage (Q), the branching factor (B), and the precision (P) as follows:
(23)Detectionpercentage=100%×TPTP+FN,Qualitypercentage=100%×TPTP+FP+FN,Branchingfactor=FPTP,Precision=100%×TPTP+FP.P shows that our method proposed is better than random guess; D, Q, and B can help us to choose the best threshold.Different thresholds were chosen, i.e.,ε from 0.1 to 0.9, and then the detected crater numbers and quality factors were calculated using our proposed method. The results are shown in Table 6. The symbol ∗ represents the optimal value.Table 6
Detected crater numbers and quality factors.
Matching thresholdεNumber of cratersQuality factorsDetection resultsDQBPTPFPFN0.1143089.8%∗76.7%0.1984.0%12012291360.2127689.7%84.8%0.0694.0%1199771380.3126489.7%85.5%0.0594.9%1199651380.4123088.6%85.6%0.0496.3%1184461530.5117087.5%87.5%∗0∗100%∗117001670.6107480.3%80.3%0∗100%∗107402630.796272.0%72.0%0∗100%∗96203750.890968.0%68.0%0∗100%∗90904280.920915.6%15.6%0∗100%∗20901128The detection percentage curve as shown in Figure13 was drawn using the data in Table 6.Figure 13
Detection percentage (D) curve.We can make the following conclusions from Table6 and Figure 13.
(1)
Whenε=0.1, although the number of detected craters is even more than manual counting ground truth, TP and detection percentage do not reach 100%. There are two reasons, the first reason is that some structures of old craters are broken, and the other reason is that some semicrater edges mismatch to other craters(2)
When the matching thresholdε increases, the constraint for matching semicrater edges becomes stronger, and therefore, the number of extracted craters decreases; TP and FP decrease at the same time(3)
When the matching thresholdε increases, FN increases along with it, because of irregularly shaped, incomplete crater edges, or craters that did not have an obvious change in gray values(4)
The CDA achieves the best result only by choosing the appropriate threshold. When the matching thresholdε increases, quality factors D and B become better but the number of detected craters decreases. When ε=0.5, there is no FP which matters with navigation, and Q and B are all the best value, which means ε=0.5 is an appropriate threshold(5)
Design an evaluation parameterE to comprehensively consider various quality factors.(24)E=D+Q+1−B+P.Whenε=0.5, E=3.75, which is larger than other parameters that correspond to other ε.The area under the curve (AUC) in Figure13 could be calculated to evaluate the performance of the proposed method; the higher the AUC, the better the proposed method.
(i)
AUC=1 means that the proposed CDA is perfect, and it can detect all craters without any FP and FN(ii)
0.5<AUC<1 means the proposed CDA is better than random guess. If the matching threshold ε is set appropriately, the proposed method can work very well(iii)
AUC=0.5 means that the proposed method works the same with random guess(iv)
AUC<0.5 means that the proposed method works even worse than random guess.The AUC of the proposed method is 0.9192, which is much higher than 0.5, which means that our method works well and has robustness for crater detection.Table7 shows a comparison of the quality factors of our proposed method with He’s CDA [38]. For the same lunar surface picture, our method obtained a better detection percentage, quality percentage, and branching factor, due to a larger number of TPs, smaller number of FPs, and smaller number of FNs.Table 7
Comparison with other crater detection algorithms.
Method (ε=0.5)Crater numberQuality factorsDetection resultsDQBTPFPFNOur approach117087.5%87.5%011700167He [38]97964.6%60.0%0.12865104472Hough transform91745.8%38.5%0.42613255724Highlight/shadow region104759.2%52.5%0.22792171545Compared to He’s method and other traditional methods, the proposed method improves TP by 35.3%, and FP and FN are much smaller than He’s method. Quality factors are also better. It can be concluded that the proposed method is a more effective large-scale image CDA.In Table8, we compared the quality factors using a published lunar crater database [24, 41]; there are 4333 craters at the landing site.Table 8
Comparison with He’s CDA.
Methods (ε=0.5)Crater numberQuality factorsDetection resultsDQBTPFPFNOur approach117088.5%88.1%0.005383519498He [38]97958.4%44.2%0.55253214011801By dividing the crater detection algorithm into substeps (image preprocessing, ternary segmentation, adaptive dual-thresholdτlτh selection, edge detection, edge matching, and ellipse fitting), we can compare the difference between the proposed method and He’s method. Therefore, we selected a subgraph in Figure 11 of 963×1247, calculated the time for each step, and analyzed the computational cost; the result is shown in Figure 14 and Table 9.Figure 14
Computational cost of different methods.Table 9
Computational cost comparison with He’s CDA.
Method (ε=0.5)Image preprocessingTernary segmentationτl,τh selectionEdge detectionEdge matchEllipse fittingTotal timeOur approach0.395 s0.257 s0 s0.353 s0.701 s0.025 s1.831 sHe’s approach0.364 s0 s1.112 s0.664 s0.842 s0.015 s3.057 sCompared with He’s method, the computational cost of the proposed CDA is significantly reduced. The time for image preprocessing is basically the same; He’s method takes a long time to adaptively select two thresholds, while ternary segmentation based on the maximum entropy threshold, as proposed in this paper, takes a short time. In addition, as compared with He’s method, using our method, the edges detected are continuous, smooth, and without broken edges, and the time for edge matching is also shorter. Because our method has less mismatching and strong robustness, the ellipse fitting time is also slightly less than for He’s method. Overall, compared to He’s method, the computational efficiency of He’s method can be improved by 40.1%.


## 5. Conclusions

In this paper, a robust crater detection algorithm based on a maximum entropy threshold was proposed. The proposed method can effectively improve the robustness of crater detection algorithms for Mars, the Moon, asteroids, and other celestial bodies and can enhance the quantity and reliability of crater extraction. It can also be used for the vision navigation tasks of Mars, the Moon, asteroids, and other celestial bodies, hazard avoidance, and path planning. In our method, the entropy distribution is first calculated using two thresholds of the detected image, and the maximum entropy threshold of the image is calculated to make the image ternary. Then, the edge of the ternary image is detected to complete the edge extraction of the crater; a multiple-indicator constraint matrix is constructed, and the extracted edge of the crater is matched and fitted by comprehensively considering the length similarity, brightness, and other normalized factors, so as to complete the extraction of the crater.A simulation lunar sand table and a lunar surface image captured by the Chang’e probe were used to verify the proposed algorithm. The simulation results show that, compared with traditional crater detection algorithms, the extraction results and quality factors of the proposed method are significantly improved, the robustness is good, TPs increase by more than 30%, and the computational efficiency is improved by more than 40%. The method therefore has good development prospects for exploration missions to Mars, the Moon, asteroids, and other celestial bodies.



---

*Source: 4872248-2022-01-24.xml*