Improved automatic impact crater detection on Mars based on morphological image processing and template matching

ABSTRACT Impact craters help scientists to understand the geological history of planetary bodies. The aim of this paper is to improve the existing methodology for impact craters detection in images of planetary surfaces using a new approach based on morphological image processing (MIP). The improved methodology uses MIP followed by template matching based on fast Fourier transform (FFT). In this phase, a probability volume is generated based on the correlation between templates and images. The analysis of this probability volume allows the detection of different size of impact craters. We have applied the improved methodology to detect impact craters in a set of images from Thermal Emission Imaging System onboard the 2001 Mars Odyssey Space probe. The improved methodology has achieved a crater detection rate of 92.23% which can be considered robust, since results were obtained based on geomorphological features, different illumination conditions and low spatial resolution. The achieved results proved the viability of using MIP and template matching by FFT, to detect impact craters from planetary surfaces.


Introduction
Exploration of the solar system has been intensified in recent years. More than 100 probes were sent in the space for planetary exploration; recent probes were sent to Mars. Besides Earth, Mars has been the most studied planet, and so huge and wide-range space-imagery data-sets dedicated to Mars are available (Alves & Vaz 2007). These images have shown the patterns, distributions and morphological structures that compose the Martian surface. Detecting and counting of these structures, especially impact craters, allows us to understand the geological history of Mars (Molloy & Stepinski 2007;Capitan & Wiel 2012;Kim et al. 2012;Luo et al. 2015).
Due to important geological information that impact craters carry, they represent the most studied features in the solar system. The study of this structure offers clues about the local composition of the surface, and the detailed analysis of their morphology and distribution can provide relevant information on the geological history of a planetary body (Wetzler et al. 2005). Counting craters in remotely sensed images is the only tool that provides relative dating of remote planetary surfaces (Urbach & Stepinski 2009). Unfortunately, these investigations are monotonous and time-consuming tasks for humans, because they need to examine manually a wide range of information available in the constantly growing imagery data-sets. Thus, there is the need of efficient and reliable algorithms to detect these structures automatically.
CONTACT Miriam Maria Pedrosa pedrosammp@fab.mil.br Automatic methods of detecting impact craters help experts rapidly obtaining information about the planet when new data are acquired. The knowledge of the position, size and depth of impact craters can assist in the mission planning and even at the time of rovers landing (Yu et al. 2014). In addition, it can facilitate the process of characterization of craters and the establishment of a degradation degree. Marques and Pina (2015) have developed an algorithm to automatically delineate the contour of impact craters using polar coordinates. Salamuniccar et al. (2014) have developed a methodology to detect craters on Phobos, while Pedrosa et al. (2015) developed an algorithm to detect craters on Mercury and Xie et al. (2013) proposed a methodology to detect craters on Moon. In the case of Mars, various crater detection algorithms (CDAs) have been proposed over the last decades to address the challenges related to the automation of impact crater detection process. Each CDA adopts a different approach such as the methodology presented by Martins et al. (2009) in which they use the boosting approach. This technique was initially used by Viola and Jones (2004) to face detection and presented relevant results when it was used to detect impact craters. Efforts have been made to combined boosting approach with other techniques for the purpose of detecting impact craters. Ding et al. (2011) used boosting and transfer learning; Bandeira et al. (2012) jointed mathematical morphology (MM) and texture features in combination with the boosting ensemble supervised learning algorithm to classification, and Wang et al. (2015) integrated an improved sparse kernel density estimator into the Boost algorithm to detect sub-kilometre craters. Urbach and Stepinski (2009) developed a CDA based on MM and supervised machine leaning to detect sub-km craters in high resolution images. Burl and Wetzler (2011) employed some supervised learning algorithms to detect craters in Viking Orbiter images. Cohen and Ding (2014) presented a method by using genetic algorithm to find a high performing subset of features for a given classifier. Jin and Zhang (2014) developed a modified adaboosting approach to detect small size craters on Mars. Liu et al. (2012) combined active learning with semi-supervised learning to build an adaptive learning system to automatically detect craters from high resolution panchromatic planetary images. Yin et al. (2015) used Gist features to detect craters on the images from Mars Orbiter Camera.
For the CDA validation, most of the time a region in the image is selected, while many other scientists use the global coverage and consequently generate crater catalogues. Salamuniccar and Loncaric (2008), for instance, proposed a method for integration of several existing catalogues to obtain a unique craters catalogue which contains 57,633 craters. Stepinski et al. (2009) have presented an automated system to catalogue impact craters using data from a digital elevation model (DEM) of Mars. Their CDA allowed to generate a catalogue containing craters with size over 5 km; Salamuniccar and Loncaric (2010) developed a methodology to detect Martian craters based on digital topography data which employs fuzzy edge detection and Hough transform (Hough 1962;Rodionova et al. 2000). Their algorithm uses several parameters to achieve greater accuracy, such as techniques based on value and direction of gradient, and morphometric measurements, which enable to define a diameter-depth relation; Salamuniccar et al. (2011) used a CDA to generate a catalogue with 125,130 craters greater than 2 km diameter.
CDAs are developed to search craters in optical and non-optical images. Non-optical image-based methodologies have been suggested to avoid certain drawbacks presented by optical images. However, DEM-based methodologies encounter problems with the low spatial resolution of the data-sets (Bandeira et al. 2007). Hence, optical image-based methodologies still present advantages over nonoptical image-based methodologies since the data provided by optical range offer the widest amount of information about Mars and other planetary bodies from solar system. This advantage has stimulated researchers to pay more attention to the development of methodologies that deal with optical images. Bandeira et al. (2007) proposed a methodology for impact crater recognition on Mars based on a probability volume created by template matching. This methodology presents a satisfactory procedure to reduce limitations caused by the low spatial resolution of images and shows efficient computational implementation. The first phase in Bandeira et al. (2007) significantly contributes to achieve success. This phase identifies possible impact craters based on the analysis of shadows in the images.
The presence and the length of these shadows depend on the angle of illumination used while acquiring an image and the degradation degree of the craters as well. Moreover, mountains, rocks, dunes and other structures also produce shadows. This limits the automated analysis of shadows, and therefore, hampers greater accuracy in the detection of impact craters.
Although degraded impact craters do not produce shadows, it still preserves their circular contours. Hence, it is appropriate to use this feature to find impact craters. So a first phase based on morphological image processing (MIP), in contrast with shadow analysis, avoids the drawbacks presented when dealing with shadows. In addition, MIP remains underexplored in pre-processing phases, and consequently, its potential for impact crater detection is almost unknown.
In this paper, we propose the use of MIP as an improved methodology for automatic detection of impact craters in optical images of planetary surfaces. Unlike other methods, we have achieved success in automatic detection of craters larger than 500 m in diameter in low resolution mosaic images, which have not been catalogued yet. The most complete catalogue from Mars contains only craters over 1 km generated manually (Robbins 2012). Therefore, our CDA has the ability to deal with large areas covered by mosaic images, which holds different environmental conditions since each image was taken individually under assorted time. Consequently, this improved methodology provided better accuracy in the detection of impact craters, and apart from that, the MIP is found to be more computationally efficient.

Data-set
The improved methodology was applied to detect craters from digital images obtained in 2002 by the Thermal Emission Imaging System camera. These images are part of a mosaic that covers a total area of about 83,300 km 2 of the Utopia plain. The images of the mosaic contain 256 grey levels with spatial resolution of 100 m/pixel. The size of each image is 500 £ 500 pixels. Our methodology was tested using 48 images which comprises several image conditions, such as difference in contrast and sharpness, saturation and blurring, presence of noise and different angles of illumination leading varied shadows. Additionally, the data-set holds different reliefs with varying crater sizes and forms, elevations and geological features. Figure 1 shows four images to illustrate some of the conditions found in the 48 images that compose data-set. The mosaic is composed by images captured at different time, so it presents difference related with illumination, which is evidenced by different shadow positions inside the craters. These variations are clearly seen in the image (Figure 1(a)). Some of images show difference in contrast ( Figure 1(b)). Besides these variations, different sizes of craters show different degradation levels. Some craters are perfectly bowl-shaped, others have an oval form, showing different amounts of shadow, but others are smaller and shallower. In addition, there are differences in relation to the sharpness within the same image, because there are some regions with noise, saturation and blurring. These features can be seen from Figure 1(c,d). These variations are important and allow us to evaluate the most efficient and independent algorithm for the inherent variations in the image.

The improved methodology
The improved methodology employs digital image processing and features recognition techniques to detect impact craters from planetary surface images. This methodology comprises three phases: candidate selection based on MIP which removes noises and finds probable craters; template matching based on fast Fourier transform (FFT) that establishes correlations among candidates and templates of various dimensions to create a probability volume (three-dimensional (3D) matrix); crater detection which analyses the probability volume to find dimensions and locations of impact craters. The improved methodology differs from that presented by Bandeira et al. (2007) since we have replaced the original candidate selection phase with our new phase using MIP. Figure 2 shows different components of our methodology. The image processing was performed using the set of morphological operators available in the MM toolbox in Matlab.

Candidate selection based on MIP
The candidate selection is based on MIP, and is intended to convey that different components of the improved methodology apply concepts of MM to perform digital image processing. We refer to MM the set of methods developed by George Matheron and Jean Serra in 1964, which transforms and analyses structures within digital images (Matheron 1975;Serra 1982;Soille 2003).
The basic principle of MM is to use a geometric form known as structuring element (SE) to probe a digital image, quantifying the mode in which this element fits within this image (Dougherty & Lotufo 2003). This element is a mask which probes the image to be processed. Usually the size adopted for SE is 3 £ 3, 5 £ 5 or 9 £ 9. Thus, it is possible to select structural information from the image, marking locations at which the SE fits within the digital image. By choosing the size and the shape of the SE, we can control the type of structural information to be selected. In other words, the SE is a known set which is compared to an unknown image set, through some transformations which take into account the neighbourhood of the pixel and the connectivity given by the SE.
MM has two fundamental operations: Erosion: The erosion of an image containing a given set A, through a SE B, is the minimum value between the translations (displacement) of the image by the SE. Mathematically, the erosion of a set A by a set B (SE) is denoted by AQB and is defined by The opposite of erosion, the dilation of an image containing a set A by a given SE B, is the maximum value between the translation of the image by the SE. The dilation of a set A by a set B (SE) is denoted by A È B and defined by A È B ¼ ðA c QB À Þ c , where c denotes the complement of the set and B À ¼ Àb; b 2 B f g . However, to improve MIP applications, erosion and dilation are often used together to compose more complex operations, such as opening, closing, morphological gradient, watershed transform, etc.
In our pre-processing phase, we apply MIP to fit SEs within digital images to select structures that seem to be craters (candidates). The selection of candidates depends on the application of a MIPbased algorithm described in four steps as follows.

Highlight
Our pre-processing phase starts smoothing images by applying a maximum and a minimum volume, defined, respectively, as v max and v min operators. Volume is a criterion that combines shape and contrast. This criterion takes into account the volume level of a component, and the sum of all areas of the level components above it, including itself. It can also be demonstrated by the integral of the umbra of the signal above the level component (Dougherty & Lotufo 2003). V max removes peaks with a volume smaller than a defined value. V min removes valleys with a volume smaller than a defined value. V max is defined as 120 and v min = 30. We have chosen these values since they allow the elimination of geological features which are not of our interest in this work, such as small dunes, stones, small holes caused by secondary impacts and so on. The example of smoothing applied is presented in Figure 3(a).

Edge detection
For greyscale images, the gradient increases the difference among grey levels and enhances the edges of the craters. The arithmetic difference between the dilation and the erosion by a SE B defines the morphological gradient (Soille 2003). The morphological gradient is denoted by r: where d represents dilation and e represents erosion. We applied the gradient (Figure 3(b)) using a disk with 5 pixels of radius as a SE, since it was the most suitable SE to enhance the edges of the craters in our images.

Filtering
The goal of this step is to eliminate objects that do not correspond to impact craters. The area-closing connects areas below a stipulated threshold. We used 300 pixels as the threshold for this filter, since this value enables filling the interior of the targets. The area-closing is defined as: where f represents the closing. We reduce noises in the images by applying the opening filter. The opening g of an image f by a SE B is denoted by g B ðf Þ and is defined as the erosion of f by B followed by the dilation with the reflected SE B : This filter uses a disk with 5 pixels of radius as the SE to eliminate small noise-grain radii. The opening filter reduces noise, but also degrades the edges of large craters. We applied the closing filter to connect these edges and recover the original shape of the craters.
The closing f of an image f by a SE B is denoted by f B ðf Þ and is defined as the dilation of f with a SE B followed by erosion with the reflected SE B : This filter uses a disk with 3 pixels of radius as the SE to better round the edges of the craters. Some groups of candidates were grouped as single candidates, since a relatively small distance among them results from the filtering. These groups were separated in the segmentation step.

Segmentation
Morphological watershed transform from markers separates one candidate from another. It starts by obtaining the markers from regional maxima which are determined by the distance function. Regional maxima consist of maximum values surrounded by lower values of grey. These maxima are determined from the distance function which calculates the distance of each object's pixel in relation to the closest pixel belonging to the background. The result of the distance function allows us to determine maximum values due to its grey scale image. The higher values correspond to the most distant pixels in the background. Next, the distance function is inverted to convert regional maxima into new regional minima. These new images containing the regional minima were used as input by the morphological watershed transform.
Hence, the watershed associated with the set of regional minima M ¼ [ i2R m i of an image f is the complement of the union of all the retention basins C f ðm i Þ (Pr eteux 1993): The watershed considers the image as a topographic surface, to see which markers work as holes punched in each regional minimum. While performing the watershed, this topographic surface is flooded from below, by letting the water rise from each hole, at a uniform rate across the image. Dams are built to avoid the merging due to the rising water coming from distinct holes. When the only visible structures above the waterlines are the top of these dams, the required watershed lines have been located and the watershed transform is finished (Frigo & Johnson 1998;Dougherty & Lotufo 2003). The overlap of each step of the watershed application is shown in Figure 3(c). The final results of pre-processing phase are presented in Figure 3(d).
The MIP step applied also removes craters that are on the boundary of the image, since only whole candidates are selected for the next step. The white structures represent all candidates in the image. The goal of this step is to eliminate objects that do not correspond to impact craters. The area-closing connects areas less than a stipulated threshold. We used 300 pixels as the threshold for this filter, since this value allows filling the interior of the targets.

Template matching based on FFT
This phase seeks circular features within image, because impact craters have roughly circular shapes. Therefore, circular shapes are used as templates to search and match candidates within images. Each template is presented as a white crown surrounded by a black square background. Figure 4 shows three examples of used templates. The diameters of the templates vary from 5 to 100 pixels, considering the external edge of the crown, since crowns of different sizes represent the distinct natural dimensions of the impact craters. The internal radius of the crown is half the external radius for all templates. There is a template for each radius. Each template must be cross-correlated with the scene of the image, from the template with 5 pixels to the 100 pixels template. This phase evaluates correspondence between templates and the scene of the image of Figure 3. There are consequently correlations between templates and selected candidates within the image. These correlations are calculated by FFT, multiplying the FT of the image by the complex conjugate of the FT of each template (Frigo & Johnson 1998;Gonzalez & Woods 2010). The FFT allows low computational effort when compared with the Hough transform (Bandeira et al. 2007). This is a very important consideration when large amounts of images need to be examined. Figure 5 presents the results achieved from correlations between the image of Figure 3 and a template with radius of 10 pixels ( Figure 5(a)), and a template with a radius of 20 pixels ( Figure 5(b)).
The results of correlations were normalized between 0% and 100%. This is performed by dividing each achieved result by the number of white pixels in the template. Each new value represents the probability p of a pixel (u, v) to be in the middle of an impact crater of radius r. Therefore, Figure 5 reveals achieved probabilities of correlations between an image and a template with a radius of  10 pixels (Figure 5(a)) and a template with a radius of 20 pixels (Figure 5(b)). Thus, from Figure 5, the higher the probability values are, the brighter are the pixels.
All information generated by the second phase process is stored in a 3D array, called probability volume (Bandeira et al. 2007), consisting of n plans, which contains u £ v elements. Each plan is generated by a different template. Thus, the probability volume consists of planes that represent the probability of each pixel being the centre of a crater with radius r. Through the probability volume, it is possible to identify the pixels with a high chance of being the centre of a crater.

Crater detection
The analysis of the probability volume created from the results of the template matching enables detection of impact craters in images. Circular shapes produce unique signatures in the probability volume. Therefore, the internal structure of this volume can be used to detect impact craters and determine their dimensions.
Once again, MIP was employed for impact craters identification using regional maxima. According to Soille (2003), a regional maxima M of an image f is a set of connected pixels having the same value t, and whose external boundary pixels have a value less than t. The regional maxima can be achieved from the morphological reconstruction by dilation R d f , which is defined as the geodesic dilation of a marker image g with respect to g iterated until stability is reached. This is shown in Equation (6) (Soille 2003): where i is such that d ðgÞ. So, the regional maximum is defined as follows (Soille 2003): The regional maxima show impact crater features and also other structural features, which can be wrongly interpreted as craters. In order to minimize this problem, extreme values are extended. Adjacent elements of the regional maxima may show probability of being in the middle of a crater or combined with other features. These extended maxima are calculated in two steps, first, in the maxima, all heights less than a given value h are suppressed, applying the h-maxima transform as defined by Equation (8) (Soille 2003): The regional maxima of the h-maxima transformation is called extended maxima EMAX h and is defined as follows (Soille 2003): 3D connectivity is considered for this calculus and each element of the probability volume is compared with a neighbour's 26 pixels. There are signatures of impact craters in the horizontal planes of the probability volume for each radius. Each radius holding extended maxima have been analysed.
The area-opening transform (Soille 2003) removes connected components whose area, in number of pixels, is smaller than a given threshold value λ, so this operator eliminates small objects to be markers of craters, based on a given proportion of r: Equation (10) shows that the area opening can be achieved by the union of all openings with connected SEs whose size in number of pixels equals λ.
The classic circularity index (CI) allows to eliminate non-circular objects (Bandeira et al. 2007): All weak crater candidates whose roundness is outside the range 0.9-1.7 were discarded. Each centre of mass of the remaining objects represents a central point with coordinates (x, y) of a crater, radius r and probability p (p corresponds to the confidence level of each detection). This phase analyses every plane of the volume where the maxima shows a probability of at least 20%, since craters are rarely expressed for probabilities lower than 20%. Figure 6 shows the final result having an image containing varying craters.

Results
The improved methodology was applied to detect impact craters automatically in an imagery dataset. The results were evaluated establishing a comparison with a manually built Ground Truth (GT). The GT enabled us to carry out a quantitative evaluation of the improved methodology. The evaluation is based on three indicators, presented, respectively, by Equations (12)-(14).
From Equation (12), TDR P shows true detection rate, TD P is the number of true detections for a minimum probability P and GT is the Ground Truth. From Equation (13), FDR P means false detection rate, FD P is the number of false detections for a minimum probability P , TD P the number of true detections and FD P the number of false detections. In addition, Equation (14) considers the false negative (FN) pixels that are the number of craters present in GT which was not detected by methodology. Bandeira et al. (2007) used Equations (12) and (13) the first time for detection of craters, and Shufelt and McKeown (1993) used Equation (14) which allowed to evaluate the quality of CDA performance. Based on the aforementioned matrix, our methodology shows detection rates for maxima with probability of at least 20%: TDR 20 = 92.23% and FDR 20 = 23.57%. In other words, the improved methodology shows detection of 3462 craters in the set of 48 images from a mosaic. It is important to mention that the CDA found craters with sizes larger than 500 m in diameter. Figure 7 summarizes the measure reached by the performance of the algorithm. The result is relevant considering the most complete catalogue from Mars which contains craters over 1 km long (Robbins & Hynek 2012).
There are 3753 GT craters in the selected area for the test. The improved methodology automatically detected 3462 impact craters with 71.82% algorithm performance showing better performance of algorithm compared with previous CDA. Considering the difficulty of detecting craters in mosaic data due to different environmental conditions and the low resolution of the data, our CDA shows improved performance for crater detection. Figure 6 shows the results of our methodology applied to four images from our data-set. Among the data-set used, the best results achieved were those shown in Figure 6(d) which achieved around 97.10% and 12.86% for true and false detection rates, respectively. The rates for a probability of at least 20% for these four images are presented in Table 1. The correct choice of probability level is important, since if the level is increased, it can lead to the elimination of false detections, but it can also reduce the true detections. Thus, the value 20 was the best value to balance TDR p and FDR p .

Discussion
In this paper, we have presented an improved methodology for automatic impact crater detection on Mars, based on MIP and template matching. We proposed the application of a pre-processing phase based on MIP to improve an available methodology for automatic detection of impact craters in optical images of planetary surfaces.
We have also shown the results achieved by the improved methodology when applied to an optical imagery data-set. Although the optical images present different conditions of illumination, varying geomorphological settings and low spatial resolution, the improved methodology shows very high detection rates. The improved methodology reached a crater detection rate of 92.23%, in contrast it shows a false detection rate of 23.57%, which is considered reasonable for this type of data, since the data-sets are from mosaic composition holding individual images taken at different times and thus they have different environmental conditions. These results are promising and it's in line with crater impact detection literature (Bandeira et al. 2007;Bue & Stepinski 2007;Martins et al. 2009;Salamuniccar & Loncaric 2010;Bandeira et al. 2012). We recognize the importance of comparing our results to the state of art in crater detection methodologies. Unfortunately, it is hard to establish appropriate comparisons between the improved methodology and other published crater detection methodologies since each one has demonstrated its performance using a different imagery data-set. However, we have used an imagery data-set similar to the data-set used by Bandeira et al. (2007) to establish a closer comparison. For this comparison, we reached better results since the improved methodology enables to detect craters for maxima below a probability of 30%, which was not found by Bandeira et al. (2007). This results mainly from the use of our pre-processing phase based on MIP. In contrast with the shadow analysis used by Bandeira et al. (2007), our pre-processing phase avoids the drawbacks presented while dealing with shadows, such as irregular sizes and shapes of shadows that limit better results.
The success of MIP as part of the improved methodology results from a suitable choice of morphological operators, SE and threshold values from among those tested. The benefits of using MIP to increase crater detection rates have been previously suggested by Bandeira et al. (2007).
On the other hand, template matching is efficient when used to find circular features within an image. This technique is important in the detection process, since craters have roughly circular shape. Moreover, FFT increased efficiency by calculating the correlation between images and templates in the frequency domain which is better compared to other techniques that carried out the same correlation in the space domain. The FFT also reduces computational efforts, which is important when examining large number of images. As a result, we are able to perform a fast, automatic and non-tedious impact craters detection algorithm. Moreover, with a single methodology, we carried out the detection of a large number of impact craters in images with great range of terrains on the surface of Mars.

Conclusions
We conclude that our work has achieved its goal, since the improved methodology efficiently applies MIP and template matching using FFT to detect impact craters, and enables to achieve important final results. By analysing the results achieved from the improved methodology, we also conclude that MIP is an efficient tool, when applied, first, to select crater candidates and to help the detection of impact craters.
Finally, our improved methodology can be applied for mapping of Martian craters in order to contribute to a greater understanding of the geological history of Mars and other planets.