CRATER DETECTION USING TEXTURE FEATURE AND RANDOM PROJECTION DEPTH FUNCTION

In this paper, a novel automatic crater detection algorithm (CDA) based on traditional texture feature and random projection depth function has been proposed. By using traditional texture feature, mathematical morphology is used to identify crater initially. To further reduce the false detection rate, random projection depth function is used. For this purpose, firstly, gray level co-occurrence matrix and a novel grade level co-occurrence matrix are both used to further obtain the texture features of these candidate craters. Secondly, based on the above collected features, random projection depth function is used to refine the crater candidate detection results. LRO Narrow Angle Camera (NAC) mosaic images (1 m/pixel) and Wide-angle Camera (WAC) mosaic images (100 m/pixel) are used to test the accuracy of proposed method. The experimental results indicate our proposed method is robust to detect craters located in different terrains.


INTRODUCTION
Craters have been used as important landmarks for high-precise landing of lander, autonomous spacecraft and rover navigation and control (Yu et al., 2014;Wang et al., 2015). Furthermore, craters also play an important role in the study of planets chronology (Barlow, 2015). The size frequency distribution of primary craters can provide the primary mechanism for establishing chronology of planetary surfaces (Head et al., 2010). In this way, an automatic crater detection method is necessary.
Over time, a number of automated crater detection algorithms (CDA) have been developed. By assuming that the shape of craters is circular, techniques such like Hough transform (Kiryati, 1991), circle fitting (Salamuniccar et al., 2011) and template matching (Flores-Mendez, 2003;Bandeira et al., 2007) are always used. Considering under certain light conditions, the sunward side of craters always present as locally brightest feature and the backside of craters present as locally darkest feature, the technique of highlight-shadow region matching is also widely used (Urbach, 2007). And the highlight-shadow feature has been proved to be an effective method for detecting craters located in various terrains (Urbach and Stepinski, 2009). In addition, with the continuous development of computer vision technology, machine learning techniques are widely used (Lienhart et al., 2002;Ratsch et al., 2001), including support vector machines (SVMs) (Suykens et al., 1999), boosting approach (Bandeira et al., 2012), decision tree (Mishra et al., 2013) and CNNs (Palafox et al., 2017). These methods do not rely on expert's domain knowledge, but depend on learning the features based on training samples, which are more efficient for detecting multi-scale caters. However, these methods request a lot of labelled data for training and their performance depends on the quality and number of training data, in addition, the model parameters are complex.
In this paper, a novel automatic CDA has been proposed. Firstly, the candidate detection method is based on the algorithm proposed in (Urbach et al., 2009) with our modifications which lead to a more adaptive results to our whole proposed method. Secondly, candidate regions are further classified as crater or non-crater regions using gray level co-occurrence matrix (Gadelmawla, 2004) and grade level co-occurrence matrix and random projected depth function. In contrast to other CDAs, the proposed method does not rely on strong assumptions about crater edge shape, and does not rely on large amount of labelled data, however, it focuses on the current scene, taking all candidate crater region as input samples, and by using the idea of anomaly detection, the candidate crater regions will be refined. Using this strategy, the real craters can be detected while the noncrater regions are rejected.

Candidate Craters Detection
The crater detection method proposed by Urbach et al. (2009) has been proved as a convenient and efficient method to detect candidate region of craters. This method can detect almost all small-scale craters but also include many false detected craters. In this way, these detected crates, including correctly detected and false detected, can be further classified by random projection depth function. We use this method with our modifications in the proposed CDA and it will be described briefly. The key insight of Urbach et al. (2009) is to consider that the crater candidates can be recognized as a pair of crescent-like highlight and shadow regions. Under the lower solar elevation angle condition, the above described feature is significant and the crater detection result is well, however, those degraded craters might be missed as they may not show clear crescent-like illuminated and shadowed patterns. In this way, the proposed method does not used shape filters as Urbach et al. (2009) did, and the main idea is to segment highlight and shadow regions from the image which satisfy the area constrain. The method identifies and processes highlight and shadow regions simultaneously by using both the original image and inverted image. Firstly, the oversize shape features are removed by using median filter. Secondly, the shape features lack of sufficient contrasts are removed by using power filter. Thirdly, the undersize shape features are removed by using mathematical morphology operators. And then, shape filters are used to identify highlight and shadow shapes which have geometries consistent with being parts of craters. Finally, the remaining highlight and shadow regions are matched to each other to form the candidate craters.

Texture Features Extraction of Candidate Craters
For these candidate craters regions, we first extract their texture features by using gray level co-occurrence matrix (Gadelmawla, 2004) and grade level co-occurrence matrix.

Gray level co-occurrence matrix
The gray level co-occurrence matrix (GLCM) is a reliable method to describe the texture of a gray image by quantifying the spatial distribution between the pairs of pixels. Given a gray level image I = {I (xi, yi), i = 1,2,…, n, (xi, yi)∈ D}, where i is the index of pixels, n is the number of pixels, (xi, yi) denotes the spatial coordinate of pixel i, D ∈R 2 is the image domain, Ii is the gray intensity of pixel i. Given an element Gd, θ (Ni, Nj) of the GLCM of an image, where Ni is the ith quantized gray level of pixel i, d is the distance of pixel i and pixel j, calculated as ((yjyi) 2 +(xj -xi) 2 ) 1/2 , and θ is the orientation of the pair of pixel, calculated as arctan((yj -yi)/ (xj -xi)). In order to generate a GLCM for an image, three parameters should be confirmed in advance, that are the distance d, the orientation θ and the quantized gray level Ng. The size of the GLCM is affected by the Ng, the larger Ng, the larger the dimension of GLCM. In this way, GLCM matrix of a gray image can be expressed as follows, , , where G d, θ (Ni, Nj) is calculated by counting how often pairs with specific quantized gray levels of Ni and Nj and in a specified spatial relationship. The Gd, θ is then normalized to make the sum of all the elements equal to one. Based on G d, θ (Ni, Nj), GLCM indices can be computed. In this paper, according to our a large number of experimental results, it is shown that for these crater candidate regions, our GLCM indices can describe the texture feature efficiently. 1) Energy (Gadelmawla, 2004) The energy is valued by the sum of squares of GLCM elements, measuring the texture uniformity, or pixel pair repetitions. When f1 is large, the texture is depth and the energy is large; otherwise, when f1 is small, the texture is fine and the energy is small. High energy occurs when the distribution of GLCM elements is constant or periodic.
2) Contrast (Gadelmawla, 2004) 2 The contrast measures the drastic change of gray level between contiguous pixels, reflecting the depth of the texture. If the texture dramatically changed, the contrast is large and the effect is clear; otherwise, if the contrast is small, the texture is smooth and the effect is fuzzy. The more pairs of pixels with high contrast, the greater the value.
3) Homogeneity (Gadelmawla, 2004) ( ) The homogeneity measures the change of the texture of local area. It is sensitive to the presence of near diagonal elements in a GLCM, which reflects the similarity of gray level between adjacent pixels. The higher the value of homogeneity is, the smaller the changes of the texture happen in the different areas of an image. 4) Correlation (Gadelmawla, 2004) , The correlation measures the consistency of image texture, reflecting the similarity of elements of GLCM in row or column direction. When the matrix element values are equal, the correlation value is large; on the contrary, if the matrix pixel values differ greatly, the correlation value is small. If there is a horizontal texture in the image, the correlation value of the horizontal matrix is greater than that of the rest of the matrix.

Grade level co-occurrence matrix
After carefully analyzing the textures of craters and non-craters images, the textures of craters can be further effectively represented by grade of each pixel within each image block, both the magnitude and direction. In this way, gradient co-occurrence matrix is proposed to extract the texture features of the candidate crater image blocks.
In order to calculate the gradient co-occurrence matrix, the gradient of each pixel within each block should be calculated firstly to form gradient magnitude image Im = { Im (xi, yi), i = 1,2,…, n, (xi, yi)∈ D} and gradient direction image Iα= { Iα(xi, yi), i = 1,2,…, n, (xi, yi)∈ D}. The gradient of each pixel is calculated by using the convolution template as dx = [-1 0 1], dy = [-1 0 1] T , and for pixel i(xi, yi), its gradient is calculated as, where Ix(xi, yi), Iy(xi, yi) is the gradient value of pixel i in horizontal direction and vertical direction respectively. And then for pixel i, the magnitude Im(xi, yi) and direction Iα(xi, yi) of gradient are calculated as follows, Based on Im and Iα,, the construction of gradient co-occurrence matrix is straightforward. Similar with the method to construct GLCM, the gradient magnitude image Im is first quantized into pre-divided grade levels, and different grade values are mapped into corresponding grade magnitude level. In the same way, the gradient direction image Iα is also mapped into corresponding grade direction level. Then, the elements of gradient cooccurrence matrix are calculated by counting how often pairs of pixels with specific gradient magnitude levels and specific gradient direction levels, that is given Nm and Nα,, how often the pairs of pixels satisfying Im(xi, yi) ∈ Nm,i, Im(xj yj) ∈ Nm,j, Iα(xi, yi) ∈ Nα, Iα(xj yj) ∈ Nα. The gradient co-occurrence matrix is then normalized to make the sum of all the elements equal to one. For a given Nα, the generic element of the matrix is noted Gm, α (Nmi, Nmi) and based on them, gradient co-occurrence matrix indices can be computed. In this paper, two gradient co-occurrence matrix indices have been selected.

1) Contrast of gradient distribution
, , 2) Homogeneity of gradient distribution Based on the above recognitions, for calculating GLCM matrix, the orientation θ is along the illumination direction and perpendicular to illumination direction(quantized into integral multiple of 45°), and GLCM matrix of each orientation has been calculated, and four statistical indices from each GLCM matrix have been extracted. For calculating gradient co-occurrence matrix, the gradient direction α is along the illumination direction and perpendicular to illumination direction (quantized into integral multiple of 45°). Finally, 12-dimensional feature vector vf (4 features of each GLCM for 2 orientations, 2 features of each gradient co-occurrence matrix for 2 orientations) of each crater candidate region has been extracted.

Random Projection Depth Function for False Craters Elimination
Analogous to linear order in one dimension, projection depth function provides an ordering for each point in a multidimensional data set from the 'center' of to 'outlier' values, that is, with respect to a data cloud or a probability distribution it assigns to each point its degree of centrality. The definition of a projection depth function was firstly proposed by Donoho (1992), and it is based on the definition of outlier. Outlier refers to the data that is far greater than or less than other data in the data set, showing obvious deviation. In general, outliers are different from gross errors. They are real and normal data in the data set, but they only show some extremes. For a given dataset X = {Xi ; i = 1, ..., n}, where i is the index of data point, n is the total number of data points in the dataset; Xi = (xi1 , xi2 , ... , xid) Τ is the data point, and T is the transpose operator. The outlier function is defined as follows (Zuo, 2003), where u∈R d ，||u||2 = 1, MED is the median and MAD is the median absolute deviation. MAD is a robust estimator of variability and is defined as follows (Zuo, 2003), Compared with traditional (µ(⋅), σ(⋅)), the pair of robust estimators (MED(⋅), MAD(⋅)) are not seriously affected by outliers. In this way, with respect to a multivariate data set, projection depth function measures the centrality of a data point (a vector Xi) as the maximum outliers of the one-dimensional scale functional in any one-dimensional projection, and is defined as (Zuo, 2003), where sup(⋅) is supremum operator.
Considering that in (14), an infinite set of random projections need to be calculated to obtain PD (Xi, X), however, it is impossible in practice. In our case, following the suggestion of (Zuo, 2006) that replacing the supremum in (14) by a maximum over a finite number of randomly chosen projection. Therefore, a stochastic approximation of the random projection depth PD(Xi, X) can be calculated by using m random projections uniformly distribute in R d as follows (Zuo, 2003), As discussed above, random projection depth function provides a vector ordering in the vector space, from the 'deepest' to 'outward' ordering. In our paper, given a vector data set V = { vf1, vf 2, …, vf i…, v f n}, where n is number of candidate crater regions, vf i is a 12-dimensional feature vector described by gray level cooccurrence matrix and grade level co-occurrence matrix of ith candidate crater region. In this way, random projection depth value for each candidate crater region can be calculated as, where u∈R 12 , m = 1000.
The random projection depth function can provide an order from the centre to the outside of the dataset, that is, the closer the data points are to the centre of the dataset, the greater the random projection depth value, and vice versa. The data set centre is the median of the data set, with the maximum random projection depth value. The order of any two data points in the dataset is: where ↔ is the equivalent operator, the relationship of equation (17) indicates the centre degree of vf i is less than that of vf j. If the vector space V satisfies the following two conditions: (1) V = {VB, VF}. that is the data set has two components which are target data VB with certain characteristics and outlier data VF with different characteristics comparing with data.
(2) VB ∩ VF = ∅ and #{VB } > #{VF}, where represents the number of elements in the dataset. Obviously, our candidate crater regions vector data set V satisfies the above two conditions, therefore, random projection depth function is an efficient way to refine the crater detection result.

Data Set
The first kind data set we selected is Lunar Reconnaissance Orbiter Camera (LROC) mages. The second kind data used for crater detection is Wide-angle Camera (WAC) "morphologic" mosaic images, which was created by mosaicking monochrome (643 nm) WAC images with an average 60° solar altitude angle.
For elevation purpose, the ground truth craters were manually identified in ArcMap software by using CraterTool. To guarantee the reliability of the ground truth detection, two operators had identified crates from these images and the detected results were combined in the end. Three indexes were used to evaluate the performance of the proposed approach which are the integrity rate (TDR), the accuracy rate (FDR), and the detection rate (DR). These three indexes are defined as follows (Umenyiora et al., 2012), where TP is the true positive and is defined as the number of correctly detected, FP is the false positive and is defined as the number of incorrectly detected, FN is the false negative and is defined as the number of missed detected.

Detection Results and Analysis
For the first data set, the selected area is located between longitudes of 199.607°E ~ 199.609°E, and the latitudes of 45.156°S ~ 15.158°S. The resolution of this image is 1m/px. This selected test areas contain diverse types and sizes of craters. It can be seen from the image the illumination is non-uniformity. In addition, the terrain of the selected area is very rugged. All these factors contribute to challenges to the proposed method to detect the crater exactly. Figure 1 shows detected craters, and the detection result is show as red circles. The total extraction is 2138, the maximum diameter of these extracted crater is 203m and the minimum diameter of these extracted crater is 5m. The detailed quantitative analysis result is presented in Table1.

Figure 1. The detection result of craters in the LROC image
From above detected result, it can be seen that our proposed method can detect different scales craters, especially these smallscale craters. At the same time, it can be seen that our proposed method is not affected by light direction and nonuniformity of light. The incorrectly detected craters are always those next to the ravine, whose texture features are similar with craters. The miss detected craters are those texture features have been distorted largely.
For the second data set, the test site is located between longitudes of 49°E ~ 69°E, and the latitudes of 41°N ~ 45°N. The resolution of this image is 100-m/px. The size of the image is 1000×1000 pixel. Most craters in this image appear as 'bowl-like'. The ground truth craters used for verifying the accuracy of proposed method is the data set made by Robbins (Robbins, 2019). The data set is a new, global database of lunar impact craters, estimated to be a complete census of all craters with diameters larger than 1-2 km. The Figure 2 shows both the detected result and craters of Robbins' data set. Given that there are more than ten thousand craters, in order to show them more clearly, all these craters are presents as their centres. Craters of ground truth craters are indicated by red dots, and craters of detected results are indicated by blue dots.
It can be seen from the Figure 2 that the amounts of craters detected by proposed method are larger than the Robbins' data set. That is because our propose method can detected multi-scale craters, for this area, the maximum diameter of detected crater is 21600m and the minimum diameter of detected crater is 500m, while the diameters of craters from Robbins' data set are all >1000m. In this way, according to diameter, we further divided our detected results into two parts that are craters with diameters > 1000m and craters with diameters < 1000m. In Figure 3, we present craters with diameter > 1km, where craters of ground truth craters are indicated by red dots, and craters of detected results are indicated by blue dots. We can find that for some of the centres of craters detected by proposed method and Robbins are totally overlapped, and some of them although not completely coincident, the location is very close. According to the detection result, the miss detected craters are those have seriously degraded, and the whole area almost has no shadow and highlight region. The quantitative evaluation of craters with diameters > 1000m is shown in Table 1. Furthermore, we cropped an area where the craters with diameter of <1000m are more concentrated, as presented in Figure 2 with red rectangle, and the detected result are presented in detail in Figure 4. The crater detection result indicates that for these terrain relatively flat areas, the proposed can detect almost all craters.

CONCLUSION
In this paper, we proposed a crater detection method based on GLCM、 grade level co-occurrence matrix and random projection depth function. By using highlight-shadow feature, candidate crater regions are detected. However, the false detective ratio is high, therefore, the random project depth function has been used to reduce the false detective ratio, and the method is efficient. In this way, we were able to verify average greater than 80% of ground truth craters. Specifically, the propose method is good at detecting multi-scale craters especially the small-scale craters. In addition, the proposed method is robust to illumination direction and illumination nonuniformity and complexity of terrain.