ASTEROID ( 21 ) LUTETIA : SEMI-AUTOMATIC IMPACT CRATERS DETECTION AND CLASSIFICATION

The need to develop an automated method, independent of lighting and surface conditions, for the identification and measurement of impact craters, as well as the creation of a reliable and efficient tool, has become a justification of our studies. This paper presents a methodology for the detection of impact craters based on their spectral and spatial features. The analysis aims at evaluation of the algorithm capabilities to determinate the spatial parameters of impact craters presented in a time series. In this way, time-consuming visual interpretation of images would be reduced to the special cases. The developed algorithm is tested on a set of OSIRIS high resolution images of asteroid Lutetia surface which is characterized by varied landforms and the abundance of craters created by collisions with smaller bodies of the solar system.The proposed methodology consists of three main steps: characterisation of objects of interest on limited set of data, semi-automatic extraction of impact craters performed for total set of data by applying the Mathematical Morphology image processing (Serra, 1988, Soille, 2003), and finally, creating libraries of spatial and spectral parameters for extracted impact craters, i.e. the coordinates of the crater center, semimajor and semi-minor axis, shadow length and cross-section. The overall accuracy of the proposed method is 98%, the Kappa coefficient is 0.84, the correlation coefficient is ~0.80, the omission error 24.11%, the commission error 3.45%. The obtained results show that methods based on Mathematical Morphology operators are effective also with a limited number of data and low-contrast images.


INTRODUCTION
Formation of impact craters is one of the physical and geological processes that enable understanding of the evolution of planets.Determination of spatial parameters of craters (shape, size, inclination) and their distribution on the surface of the considered objects is the basic method of remote determination of its age (Soderblom et al., 1974).The problem of developing a reliable algorithm for impact craters detection raises great interest among the researchers who work on: autonomous landings on asteroids (Leroy et al., 2001) and planets (Johnson et al., 2005), the proper identification of rovers location using surface data (Li et al., 2002) and their autonomous navigation (Maimone et al., 2004), co-registration of data projected in different coordinate systems (Michael, 2003), searching for craters on the surface of the Earth (Chicarro et al., 2003), and advanced statistical analyses (Salamunićcar, 2004;2009).Automation of the crater detection method is a practical solution for surfaces which are heavily impacted by craters, i.e. objects where the manual mapping of smaller craters is an extremely complex and time-consuming task.The need for such a solution is expressed in numerous publications addressing this subject.A set of 73 algorithms for crater detection was presented in G. Salamunićcar, S. Lončarić, 2008a and then supplemented by another 12 by the same authors (Salamunićcar, G. and Lončarić, 2008b).Crater detection studies can be classified by the type of the method used, i.e. supervised and unsupervised (autonomous).The initial step in the unsupervised crater detection process is the application of edge detection algorithms (Leroy et al., 2001;Honda et al., 2002;Cheng et al., 2003;Stepinski et al., 2007)0, image texture calculation (Barata et al., 2004), or combination of edge detection techniques with texture calculations (Bandeira et al., 2011).
In addition, techniques for combining edge detection with their spatial directional analysis (Kim, Muller, 2003) and patterns matching (Magee et al., 2003) are widely used.In publication Knežević et al., 2008 an overview of the crater detection algorithms based on the gradient edge filters is presented.The processed image, might be automatically analysed using Hough Transformation (Hough, 1962), genetic algorithms (Honda et al., 2002), or radial coherence algorithms defining areas of rotational symmetry (Earl et al., 2005).Even though unsupervised crater detection methods are widely used, they might provide worse results than supervised/semi-supervised methods which are based on the concept of machine learning.The learning phase of these methods is based on a training dataset containing craters marked in a process of visual interpretation performed by an expert in the field.In the detection phase, the created classification model is applied on whole dataset, i.e. the image or the series of images, in order to extract structures which are not defined in the training dataset.It is worth noting that applying both unsupervised and supervised (or semi-supervised) crater detection approaches includes the step of successive narrowing of the initial dataset of the potential candidates.In unsupervised methods from the initial candidates dataset the erroneous results are removed by applying a threshold value defining the alignment of the detected crater to the circle or to the ellipse.In supervised methods, the filtering of the initial candidates dataset is performed by applying additional classification with the predefined threshold value indicating the positive results.Therefore, regardless of the method chosen, a very important step is to select properly the features that describe all possible candidates and the criteria that accurately extract the structures of interest from other structures present in the candidate dataset.
The growing need to develop an automatic method of identification and measurement of impact crater parameters as well as the creation of a reliable and efficient tool gave rise to address this issue in our research.In this paper we present the methodology for detecting the sub-kilometer impact craters based on their spectral and spatial characteristics.The aim of the analysis is to develop and evaluate the satellite data processing algorithm.In this way, time-consuming visual analysis would be reduced to the special cases.The developed approach used Mathematical Morphology (MM) operators both in the selection of candidates dataset and in the results filtering.The developed crater detection algorithm based on the approach presented by Urbach and Stepinski, 2009.The authors described the possibility of detecting the Martian craters thanks to the lighting conditions prevailing on the planet.With the constant direction of illumination and the negligible atmosphere, the craters on the Mars are projected on the panchromatic images in the form of shaded and illuminated paired areas.This observation gives a base for the proposed crater detection algorithm.Authors proposed the identification of shaded and illuminated craters by using Mathematical Morphology technique, and eliminating erroneous results from the candidates dataset using 1 st Hu's Moment Invariant.The approach presented by Urbach and Stepinski was developed for planetary imaging with a five times higher spatial resolution than the one used in the following data analysis.Moreover the surface of Mars is characterised by a constant inclination close to zero gradient.Therefor the method proposed by the authors requires taking into account the characteristics of the images acquired during the "fly-by"of the spacecraft.The dynamic change of observation geometry, the object strongly varied surface configuration, and the highest resolution of the data used (i.e.~ 60m), impose the modification of selected method and adjustment to the new research area by implementing the Mathematical Morphology operators in initial phase of image processing.Innovation of the proposed approach consists in using the operator of the Mathematical Morphology already in the initial phase of image processing, in order to visualize the elements of darker and lighter images characterized by high contrast with the surroundings.In addition, the proposed approach enable to classify the extracted craters in refer to their sizes, the surface parameter and to characterise the topology of analysed structures by applying the containment criterion in morphological processing chain.

TEST SITE AND INPUT DATA
The objects of interest are impact craters.The study is based on data acquired by the OSIRIS camera on board Rosetta (Keller et al., 2007).On its way to comet 67P/Churyumov-Gerasimenko, the spacecraft flew by two asteroids, one of them was (21) Lutetia, the largest asteroid ever visited by a spacecraft.During the fly-by, the OSIRIS imaging system took more than 400 pictures through its wide (WAC) and narrow (NAC) angle cameras, covering approximately 50% of the asteroid surface.For our analysis we used only NAC images (Table 1), acquired in a time period covering the closest approach (3170 km).Three images from the input dataset were selected as the training data in a process of algorithm development.
All data are panchromatic of the spectral range covering the orange (590-620nm) and red band (620-750) suitable for reflectance analysis.The asteroid under investigation is very heavily covered by impact craters, and the density of spatial distribution of craters varies significantly in different regions.The selected input dataset is characterised by the highest spatial resolution in the observed hemisphere, revealing the very complex morphology of the asteroid and its craters in great details.The phase angle varying from 51 o to 87 o allowed us to measure precisely the crater parameters in different viewing conditions, for better accuracy.The analysis was narrow to the Achaia -Ac1 (Figure 1) geomorphological region as it is one of the two oldest region of the asteroid and thus the most heavily covered by the impact craters.There are craters of various types and sizes.In this region (Ac1, Ac2) a group of researchers from the OSIRIS camera team has visually identified 153 craters over 2800 km2 (Sierks et al., 2011).Most of the detected craters have a diameter of about 10km, while craters with a diameter below 1km are structures rarely found in this region.

CRATER DETECTION ALGORITHM
The proposed methodology consists of three main steps: characterisation of objects of interest on limited set of data, semi-automatic extraction of impact craters performed for total set of data, and finally, creating libraries of spatial and spectral parameters for extracted impact craters.

Reference dataset creation and craters features description
The reference dataset is created in the process of preliminary visual interpretation of available satellite imagery.The sub-kilometer craters are manually indicated by the expert.Visual interpretation is a supervised part of proposed method.It is performed once by two independent interpretations.The Achaia region was overlaid with a grid of 3km×3km in order to select randomly representative group of data.Figure 2 illustrates the location of 49 randomly selected cells, i.e. 35% of the total number of cells.The visual counting of impact craters was based on the caters characterization described in the next Subsection (3.2).As a final step, each individual crater was represented by a single central point and counted with the aid of GIS software.Based on available sources of information, the following types of craters are distinguished in the analysed area: -Regular-type craters, characterised by the shape of a cylindrical depression with a variable-height rim; -Craters with flat bottom, filled with material left after the collision; -Clusters of craters, an example of craters whose rims have been deformed as a result of post-collision; -Craters affected by landslides.Additionally, the following remarks were made regarding the characteristics of the analysed area and the objects of interest: -Lutetia is an asteroid heavily covered with the craters; -The clusters of larger craters are regularly distributed in the study area; -Landslides, clusters of boulders or rock outcrops were not identified in this region; -Numerous linear structures of fixed orientation are noticeable in the area of the Achaia region; -Craters in the analysed area do not show significant differences compared to craters visible on other rocky planets; -The test area is approximately a flat area.Asymmetrical craters occur within the external borders of the distinguished region, i.e. on the "slopes", where the surface of Lutetia strongly changes its curvature; -Impact craters can be classified based on their size, shape, length of the shadow line (depth) and diameter; -The minimum diameter of easily identifiable structures is ~600m.It is difficult to identify smaller structures due to image spatial resolution limits.The same threshold value was set for the visual interpretation described in Sierks et al., 2011; -The maximum diameter of analysed structures does not exceed 10km; -The analysed structures appear in the images as objects consisting of bright regions on a darker background (illuminated part) and dark regions on a lighter background (shaded part).Such dependence is directly related to the nature of the spatial construction of the crater and physical phenomena occurring within its surface; -The shape of shaded and illuminated part of craters is similar to the half-moon or semi-circle and it is often subjected to transformation due to distortions resulting from the geometry of the observation; -The shaded and illuminated parts of crater are located on a straight line with the direction corresponding to the light source; -The shaded part is always closer to the light source; -Areas should have a similar shape and a defined maximum distance depending on their surface.

Candidates selection
The supervised part, described in Subsection 3.1, is performed only once on the limited dataset, therefore the whole method is referred to be a semi-automatic.In the next stage of the analysis, i.e. unsupervised (automatic), the spatial and spectral characteristics of reference dataset are used.Based on the impact craters feature description two processing guidelines are created, i.e. general and detailed.The first one enables to identify the structures that most probably could be craters, while the second one allows proceed the filtering of craters candidates, ensuring more precise identification of the analysed structures.
To the general processing guidelines the following observations are included: (a) The area of the crater is identified with two forming it and highly contrasting with the surrounding regions: the darker region -a result of the shading effect, and the brighter regiona result of the illumination effect.These areas are adjacent to each other and they are highly contrasting with the surrounding.Such dependence is directly related to the nature of the spatial construction of the crater basin and physical phenomena occurring within its surface.For further description the following symbols will be used: "B"for brighter region and "D"for darker region.The starting point in the developed algorithm is the application of input data decomposition into a set of binary images Vh where the foreground values are pixels satisfying the condition described by Eq. ( 1) for a given gray level value h.
In that sense image can be also represented by two groups of connected components, i.e. "peaks" and "valleys", pixels with a value higher or lower than a given threshold h.
In order to detect the "peak" and "valleys" which are fulfilling the minimal contract condition, we applied the morphological reconstruction by dilation Eq. ( 2) in parallel to the decomposed image Vh and to its complement Vh C , respectively.
The use of the MM operator already in the initial phase of the analysis, allows to extract B and D regions characterised by a high contrast with the surroundings, which is an important information in the process of detecting impact craters.The obtained objects were then subjected to analysis using the characteristics of impact craters described in Subsection 3.1.The fulfillment of the given area and containment criterion is performed with the use of MM operators listed in Table 2.This morphological processing results in B and D region mask obtaining.The detected structures are in line with the requirements for the minimal contrast, the minimal and maximum area and the shape containment.

. Craters characteristics and the corresponding Mathematical Morphology operators
The automatic classification of detected structures is carried out by applying the operations described in Table 2 using two different thresholds for the maximum and minimum area.This approach allowed to obtain two independent sets of detected structures (B and D regions) that meet the following characteristics: (a) large structures, with an area of 100-600 pixels; (b) medium structures, with an area of 15-99 pixels.
A more detailed classification based on the area of detected structures takes place in the next stage of analysis.

Candidates filtering
Results of initial identification of objects of interest, called candidates, must then be filtered following the detailed processing guidelines, which are a continuation of the general description of objects described in Subsection 3. The firs requirements, the shape criterion, is verified by applying the specific measure of object elongation -1 st Hu's Moment Invariant (Hu, 1962).With an appropriate selection of the threshold value, it allows the elimination of linear structures (Flores-Méndez, A., Suarez-Cervantes, 2009).All objects with the 1 st Hu's Moment Invariant less than specified threshold are eliminated from further analysis, an example is presented in Figure 3.  3, Section 5.

IMPACT CRATERS SPATIAL PARAMETERS DETERMINATION
Determination of spatial parameters, shape and size of impact craters will enable the creation of a rough characterization of the examined surface, thus constituting a significant contribution to the geomorphological analysis of the object under consideration.
The process of determining the spatial parameters is an automatic process.The preliminary identification of objects of interest on the full dataset characterised by the highest spatial resolution (Table 1) is based fully on the methodology described in Section 3.Each of the 14 images was analysed in order to identify 8 main areas of interest, presenting selected impact craters (Figure 4).The selected set represents the terrain under different lighting conditions changeable in the time series data.

Carter mask creation
Due to the spatial character of the analysed structures, it is extremely important to determine the entire area representing the crater in the time series data (Table 1).For this purpose, the algorithm of automatic, quantitative crater detection algorithm was applied, followed by a series of mathematical operations.Both morphological operations and further mathematical calculations were carried out in the MATLAB R2012a programming environment.This process consists on the following steps: 1. extraction of the corresponding B and D regions (Figure 5a); 2. determination of carter geometry elements (Figure 5b and 5c): extremes of B and D regions, considering the direction of incident angle; centroids of B and D regions; crater's center as the geometrical center of the line connecting the centroids of B and D regions.3. calculation of the parameters of the inscribed ellipse; 4. detection of image pixels corresponding to the area determined by the entered ellipse (Figure 5d).

Crater profiling
Using automatically determined extremes of B and D regions, it is possible to identify the total area of the crater.The inscribed ellipse determines not only the approximate boundary of the crater but also its diameter based on the semi-minor and semi-major axis length.
Provided results enable to determinate the length of the shadow line in the image plane and local changes in the detected reflectance both in a single image scene and in a time series dataset.In addition, the minimum, maximum and average reflectance values are calculated for the region B and D cross-sections defining the shading area (shadow line) and the lighting area.The following profile are presented below.Based on the criterion of the minimal distance and mutual orientation, the obtained results, i.e.B and D regions, were saved in a form of matrix.Each of the detected pair of B and D regions has a separate, unique label.This record makes it easy to identify shaded and illuminated areas of the same crater.The geometrically determined center of the detected crater has been design as the geometrical center of the line connecting the centroid of the B and D region.The point, with coordinates corresponding to the geometrical center of the detected crater is simply its marker in the resulting database.This marker is counted for the entire analysed area, creating a database of the automatically detected impact craters.Each of the detected manually and automatically craters was assigned automatically to one of three size groups that meet the parameters described in Table 3.The final results of the semi-automatic detection of impact craters with respect to the size of detected structures are presented in Table 4. Additionally, counting of crater's markers takes place only inside randomly selected cells over the analysed area.These data will enable the perform the comparative analysis between the automatically obtained results and the reference dataset.

Diameter
In Subsection 3.1, the visual interpretation process for reference dataset collecting was discussed.In order to ensure the credibility of the created reference dataset, the visual interpretation was performed by two independent interpreters.The correlation coefficient between those two independent results for 49 randomly selected cells is equal to ~0.95.
In the comparative analysis, a statistical regression model was used -Generalized Linear Model -GLM, (Lavrakas, 2008).The number of structures determined in the process of visual interpretation (dependent variable) was compared with the number of structures determined in the automatic process (independent variable) for each of 49 randomly selected cells.On this basis, a model was created showing the relationship between the variables used.The correlation coefficient determined for automatic results and visual interpretations for 49 cells is equal to ~0.80.The obtained values of the correlation coefficient for two samples are statistically significant at the level of α = 0.05.Additionally, in order to ensure the accuracy of the developed method, the correctness of determining the illuminated and shaded areas (B and D regions) was verified.In total, 12250 measurement points randomly scattered in the reference area were verified.Which is 10% of the total number of pixels in the reference set (49 cells with an area of 2,500 pixels).In addition, the condition of a fixed number of points per cell, i.e. 250 measurement points, has been preserved.In that way the set of measurement points is a statistically representative set (Congalton, 2009;Tillé, 2006).The results of the analysis are presented in Table 5.The detailed analysis of the obtained results showed the significant occurrence of errors of omission, false negative (24.11%) in relation to errors of commission, false positive (3.45%).The reason for that is the lack of sufficient image contrast.Such a situation takes place mainly due to a change in the curvature of the observed surface, in the shaded area, i.e. on the edges of the asteroid or due to the shape of the surface of the crater itself, i.e. flattening inside the basin.Another reason for this type of errors is incorrect pairing of the illuminated and shaded areas that create the crater.
The criterion of compliance with the direction of lighting was not been met.The obtained results show very good compliance with reference data (Overall Accuracy: 98,5%, Kappa: 0.84) (Congalton, 2009).

CONCLUSIONS
The obtained results show that the methods using the Mathematical Morphology operators are effective also with a limited number of data and low contrast of the image resulting from unfavorable lighting or physical surface properties.In comparison to object-oriented analyzes, the use of MM reduces the number of parameters, is characterised by easy adaptation to various environmental conditions and a short process time.It should be noted that the MM operators used in our studies provide the reliable results being applied on panchromatic imageries.Mathematical Morphology proved to be an extremely useful tool for impact craters identification and analysis.It has been demonstrated that the approach focused on the spatial characteristics of the objects of interest using Mathematical Morphology is a comprehensive one.It is used in the analyses relating to both the geometry and reflectance values of investigated objects.The method based on this type of solution is independent of the aforementioned parameters, computationally efficient and open to modification using a variable parameter such as: Structural Element.This element can take any shape and size, which allows detecting specific features of the image.In contrast to typical point transformations, morphological operations do not transform the whole image, but the part which surroundings correspond to the pattern of the constructed Structural Element.An important feature of the chosen method is the fact that the basic morphological transformations can be combined, which gives the basis for building complex tools of image analysis, enabling their easy extension or modification.A properly constructed Structural Element enables automatic detection of craters, and then determination of their shapes.The developed method aims at providing information on craters spatial features, which indirectly constitute an input in the study of geophysical parameters of the analysed object (e.g.reasoning about the nature of surface craters).The analysis of the illuminated ad shaded areas allows determination of important spatial parameters of craters, i.e. diameter, shadow line length and profiles.In addition, automatic detection of crater mask in the time series of 14 images of Rosetta probe allows to track changes in the projection of crater image.
The obtained results can be applied as an input data to surface optical properties determination task, the Hapke function parameters calculation (Hapke, 1993(Hapke, , 2002)), i.e. phase angle, angle of incidence, emission angle and coordinates of the position of individual image pixels in the Lutetia reference system.This would be an important contribution to the process of determining a number of interesting physical parameters, such as the local normal vector to the surface in the middle of the crater, the depth of the crater, or parameters defining light scattering in the crater.This topic is foreseen as a further direction of our research.

Figure 2 .
Figure 2. Randomly selected cells and the results of the visual interpretation carried out in Achaia region (b) B and D regions are characterised by a similar, compact, and non-linear shape.(c) Based on the area of region B and D and their distance, it is possible to classify the detected structures into two classes: large craters and other craters (small and medium).This approach allows to reject false results by taking into account the threshold distance of the detected areas.This mutual distance would depend on the area of the detected regions B and D (at constant observation conditions).
2. The detailed processing guidelines refer to the detected regions B and D. They enable to eliminate the false positive results from the set of craters candidates.To the detailed processing guidelines the following observations are included: (a) Regions B and D are most often crescent-or semicircle-shaped, distorted due to the geometry of the observation; (b) Regions B and D should meet the condition of mutual orientation consistent with the source of radiation.Region D will always be closer to the Sun.The centroid coordinates of region D must satisfy the condition of spatial location in relation to the centroid coordinates of region B.

Figure 3 .
Figure 3. Example of objects which did not fulfill the shape criterion, (a) "illuminated" regions; (b) 'shaded' regions Referring to the second observation presented in this Subsection, the spatial relation of the detected B and D regions is examined.For each of the region D, the nearest region B is automatically identified.For such a pair of objects, the control of the distance and the directionality, i.e. the angle of inclination of the straight line connecting these two areas, is carried out.The prerequisite for the proper detection of the corresponding areas is the condition of the minimum distance and mutual location consistent with the direction to the lighting source.At the candidates filtration stage the preliminary classification of region B and D is maintained.The objects of interest are divided into the large and medium group (Subsection 3.2).Only the final step of the analysis separates from the set of medium objects those which have the smallest area.In this way the group of "small" B and D regions are created.This task is performed as follows: medium structures are processing with the morphological dilation operator.The results contain corresponding pairs of objects overlapping each other and those which are still disjointed.The first group of object remain the medium structures but the second group become a new class of small structures.The final result of the implemented approach are binary images, i.e. masks, representing automatically identified objects classified into the three main groups presented in Table3, Section 5.

Figure 5 .
Figure 5.The example of automatic crater mask detection for the area of interest: (a) B (green) and D (red) regions; (b) and (c) crater's center (red), centroids (blue) and extrema (green) of B and D regions; (d) ellipse inscribed in the detected extrema (crater's mask)Figure6shows the detected crater mask (Figure4are no. 7).This example perfectly presents changes in the geometry of the crater mapped on the image plane (including shadow line changes) resulting from the variable conditions of observation (the geometry of the OSIRIS camera position relative to the Lutetia asteroid and the light source).

Figure 6 .
Figure 6.The example of the automatic crater mask detection proceed for 4 successive images from a time series dataset -Table 1: (a) image no.1; (b) image no.5; (c) image no.6; (d) image no.14.

Figure 7 .
Figure 7.The example of detected crater profile plot

Table 1
. Input dataset (the training dataset is presented with a gray)

Table 5 .
Accuracy assessment of developed methodology