An Liu,Donghua Zhou,?,Lixin Chen2,and Maoyin Chen

1.Department of Automation,Tsinghua University,Beijing 100084,China;2.Beijing Institute of Special Electronic-mechanical Technology,Beijing 100012,China

Saliency detection and edge feature matching approach for crater extraction

An Liu1,Donghua Zhou1,?,Lixin Chen2,and Maoyin Chen1

1.Department of Automation,Tsinghua University,Beijing 100084,China;
2.Beijing Institute of Special Electronic-mechanical Technology,Beijing 100012,China

Craters are salient terrain features on planetary surfaces,and provide useful information about the relative dating of geological unit of planets.In addition,they are ideal landmarks for spacecraft navigation.Due to low contrast and uneven illumination,automatic extraction of craters remains a challenging task.This paper presents a saliency detection method for crater edges and a feature matching algorithm based on edges information.The craters are extracted through saliency edges detection, edge extraction and selection,feature matching of the same crater edges and robust ellipse fitting In the edges matching algorithm, a crater feature model is proposed by analyzing the relationship between highlight region edges and shadow region ones.Then, crater edges are paired through the effective matching algorithm. Experiments of real planetary images show that the proposed approach is robust to different lights and topographies,and the detection rate is larger than 90%.

crater,automatic extraction,visual saliency,feature matching,edge detection.

1.Introduction

Impact craters are commonly geomorphic features standing out visually on any planetary surface.The data about craters including the distribution,pattern,morphology, quantity and size can yield important cues about the past and present geological processes of the planetary surface [1].Crater size-frequency distribution estimates the relative ages of terrain units[2].The knowledge of crater morphologies leads to a number of outstanding issues in planetary geomorphology,such as the nature of degradation processes,regional variations in geological material and distribution of subsurface volatile.In addition,craters are the main objects of hazard detection and avoidance for spacecraft safe landing[3].Because of their unique geometry and relatively established appearance under different conditions,craters are selected as ideal landmarks for spacecraft navigation and localization[4].However,unevenillumination,ruggedterrain and diverse characteristic of crater structures(such as degraded or incomplete and vague craters)lead to various difficultie in crater extraction from remote sensing images.Therefore,efficien detection of craters from planetary images remains a challenging problem[1,10].

Over the years,crater extraction and detection algorithms(CADs)have been widely addressed in the literature.The existing research efforts on crater detection can be mainly divided into four categories:manual detection, appearance contour fittin algorithms,machine learning methods,and geological information based analysis methods using digital elevation maps(DEM).Manual detection focuses on large craters with diameters larger than 5 km on planetary images[1].The appearancecontour fittin algorithms can be based on circular or elliptical edge fittin [3,4],including Hough transform and improved algorithm [5,6].Yang et al.developed a practical algorithm by using lighting source direction and edge convex properties, which achieved sub-pixel accuracy but lacked in adaptive ability to different imaging conditions[4].Sawabe et al. took the crater image as a shady and sunny pattern,and combined multiple ways to improve the adaptive ability to detect crater of different sizes[5].Supervised learning methods can be also applied in the crater detection,such as template matching[7],boosting algorithm[8],decision tree algorithm[9].DEM data can better represent three dimension forms of craters,including slop information[11]. However,the DEM-based crater detection algorithms fail to detect small size craters with low resolution[12].

In view of contour fittin algorithms,edge information is influence by illumination and image quality easily,and it is difficul to extract the real crater rims because of fragmentary and rough edges.Previous methods did not con-sider the visual saliency of craters morphology.For planetary images,salient regions usually contain medium and conspicuous impact craters.They are the most important visual contents according to human attention.The detection of saliency is very helpful in reducing the impact of uneven illumination and making the extracted features more distinctive.

In this paper,a saliency detection idea is introduced for planetary and lunar images.The algorithm is based on the spectral saliency detection approach while preserves a fast processing speed.Our approach differs from the existing works,in which we consider not only the human attention but also the specifi featuresofplanetaryandlunarimages.

The outline of the paper is as follows.In Section 2, a saliency detection method based on the crater feature modelis presented.Section3describesthe cratermatching algorithm in details.In Section 4,experiments using real planetary images testify the feasibility and validity of the proposed method.The main conclusions are summarized in Section 5.

2.Salient crater edge feature model

Salient object refers to the visual attended parts among the entire image,and is a powerful tool for object detection and image editing[14].For planetary and lunar images, visually attractive regions usually contain unique impact craters.Although craters may have a very wide range of variation(shown as Fig.1),a typical crater in an image has an elliptical rim and a bright to dark shading pattern.The edges of the craters are remarkable characteristic.

Fig.1 Crater images with different illumination

Since planetary images lack color information,in general,lowcontrastorfuzzyedgesaredifficul to distinguish, which may be potential crater rims.Moreover,crater rims are discontinuous and rough,which can be extracted by edge detection based on sharp grayscale gradient of the image,and also have gaps and holes.The detection of saliency is very helpful in reducing the impact of similar backgrounds and making the extracted features more distinctive.

2.1Salient edge detection method

Various contrast-based saliency detecting methods for images have been proposed,such as Itti’s model[13].A visual attention based model is used for ship detection[14]. These methods are designed for general images.They do not consider the specifi visual features of planetary and lunar images,therefore,the detection results are usually not satisfactory.Besides,CADs do not consider the human attention,and the results miss many insignifican small craters.

In our research,the salient regionbased image signature is used to get the highlight and shadow regions of craters for the candidate edges.We also generate a feature vector for each candidate edges and then use the feature matching to construct the crater rims.Finally,a robust least square fittin algorithm is applied for ellipse fittin to every matched crater rims.

If a crater is visually salient,usually both of its highlight and shadow regions are outstanding in the image and salient edges or rims can be detected.Edges informationof a crater can effectively retain the important structure of the crater,and eliminate the irrelevant information,especially, the highlight rim and the shadow edge appearing in pairs for craters.

Here,we use the image signature operator to extract salient regions,i.e.highlight regions and shadow regions [15],which are potential craters’regions.The image signature yields a saliency map through the sign information on the basis of the discrete cosine transform(DCT)for an image.It removes the spectrum amplitude of each DCT component,leaving only its symbol.According to our experiments,it can improvably detect the low contrast edges in uneven illumination and gross noise in craters’image. The concrete algorithm will be described in details in Section 3.

2.2Edge based feature model

Generally,the crater is a bowl shaped depression originated from collision of meteoroids on the planetary surface,and has a wide range of appearances.For example, some might have very vague rims and others might have broken rims.However,a classic crater in a gray-level image has a nearly circular or elliptical rim and a bright todark shading pattern,which is formed by the lighting azimuth and elevation as well as its own terrain.Shown as Fig.1,due to the angle between the sunlight and the planetary surface,only the part facing the sunlight is illuminated.Alongthesunlightingdirection,theilluminatedpart appears to be the bright region of a crater in the image,and the back part appears to be the dark region.Furthermore, since edge information just depends on the boundaries of highlight and shadow areas,they are becoming prominent. Thus,edge information description and extraction are critical in the crater detection here.

Observed from Fig.1(a)and Fig.1(b),for a salient crater,at least one highlight edge and one shadow edge can be detected,in spite of non-edge in pairs being detected.Therefore,the gathering edges in pair can describe the overview visual content of a planetary or lunar image. Thecombinationofnearbyhighlightandshadowedges according to certain rules may figur out one salient crater. And,the edges set of these craters form the saliency regions.

In this paper,an edge description of crater is proposed based on the crater illumination region feature and its relationship with edges in pair.Defin the region feature including a salient edge(shadow or highlight)as follows:

where,EC=(xce,yce)corresponds to the coordinates of the center of mass of the edge region,EOcorresponds to an orientation angle of the edge region,EPcorresponds to a perimeter of the edge,andEGrepresents the mean of all the intensity values in the region.If the coordinates of the geometry center of the edgeEC=(Xce,Yce)are calculated by

wherexe(i)andye(i)are the pixel point coordinates of the edge,respectively,andnis the number of pixels in the edge.

Then the crater feature model of edge(CME)information can be define as follows:

whereEDSandEDHrefer to the matched bright and shaded edges that construct the CME.The mass center point of the shadow region is(XceS,YceS),and (XceH,YceH)is the one of the bright regions.The other three factors are based on the feature vector of edge regionEDAtt,define as follows.

Dt:Distance factor.The distance between the centroid of pairing shaded and bright edges is calculated by

Az:Azimuth angle.The inclined direction of the paired highlightandshadowis approximatelyparalleltothe lighting direction.Letθbe the angle between the line from the mass center of shadow region(XceS,YceS)to the mass center of bright region(XceH,YceH)and the horizontal line,calculated according to the following equation:

Pm:Ratio of the perimeter of the shadow edgeEP(S) and the highlight edgeEP(H)of a crater is calculated by

Gy:Ratio of the mean intensity of the shadow edgeEG(S)and highlightedge regionEG(H)of a crater is calculated by

The purpose of CME is to match edges of the same crater.CME can describe the shape information,and take account of the image intensity profilbetween the two edgesofcraters.Thus,CME is notonlysuitableforlowsolarelevationangleillumination,butalsorobustformedium andhighsolarelevationangles.With the saliencydetection processing,CME can detect some of the vague and small craters that are more than ten pixels.

3.Crater extraction algorithm

Generally,crater extraction is based on the saliency detection with six main phases:preprocess,saliency regions detection,adaptive edge detection,candidate edge selection, edge matching and robust ellipse fitting

3.1Preprocess

In the preprocess phase,original planetary images are smoothed and contrast adjusted,which aims at removing the background noises to prevent these noise points from extracting false edge points in the next step.For this purpose,a geometric mean filte is applied to the input image forremovingbackgroundsuch as mountainsand ditch, which can be a better preservation of edge details.

3.2Saliency detection

The appearance of planetary and lunar images may be differentformanyreasons,suchas sunelevationandazimuth. In some cases,the highlight and shadow regions are not obvious.

Due to uneven grayscale change and noises contained in the image,it is difficul to detect the crater edge correctly.We introduce the image signature to robustly detect the conspicuous crater rims,which is as a simple yet powerful descriptor of the planetary image.The image signature has the ability to provide frequency analysis approach while preserving a fast processing speed.

The method for extracting salient regions is based on an image signature detector[15].The image signature is define as the symbolic information of DCT coefficient to extract salient regions of the image.

The image signature indicates salient craters whether it is highlight or shadow.A saliency map is obtained only by reconstructing the image signature,i.e,employing the inverse DCT transformation for the symbolic information of DCT.Formally,the image signature is define as

where I(i,j)is the given image,sign(?)is the sign operator,indicates the DCT of an image I.The reconstructed image takes the sign of the image I(i,j)in the transformed domain,and then inversely transforms it back into the spatial domain.

Fig.2 Saliency map of craters

3.3 Adaptive edge detection

Planetary images containing craters appear in different characteristics when conditions are different such as illumination(solar elevation angle),local terrain condition of crater(the size and depth of crater)and light shooting condition(see Fig.1).

Edge information is the most important information for the edge based crater detection algorithm.The edges of an image are extracted by applying the Canny edge detection algorithm[16].Otheredgedetectionmethodsareproposed in the transform domain[17].In general,the Canny operator has been mostly put into practice because of its excellent performance.The Canny edge detection operator uses a Gaussian filte function to smooth the image and suppress the noise,and the low and high thresholds have to be used to segment and connect edge pixels.Dual thresholds are critical to the performance of Canny detector.The traditional Canny detector gets a high threshold from the gradient cumulative histogram,the ratio of strong edge points to total pixels is 0.3,and the low threshold is usually 0.4 times of the high threshold[16].Thus,the mechanical selection can always result in that non-edges pixels may be wrongly detected as edge pixels in different images.Because the thresholds of Canny operator are set by experience,the traditional Canny detector does not have adaptive ability,and shows different performance to images in different conditions,which affects the crater detection in the end.

An adaptive threshold is proposed in edge detection algorithms with wide tolerance to imaging conditions for crater detection.The double thresholds are adaptively determined by using the gradient magnitude histogram (GMH)of the incoming image.

In a planetary image,the GMH usually exhibits single peak distribution and has an evident depression in low gradient magnitudeareas.Thus,obviousvalleys may not exist in these kinds of histograms.In such situations,the histogram concavity technique[18]is used to locate the high threshold Thand the Otsu’s method[19]is applied to determine the low threshold Tl.

The high threshold is selected to detect weak edges as much as possible,moreover,the low threshold is selected to extract the spurious edges as little as possible.If the gradient magnitude g(i,j)of pixel point(i,j)is greater than the high threshold Th,it is absolutely determined as an edge point.It is completely not an edge point when the gradient magnitude g(i,j)of point(i,j)is less than Tl. For these points whose gradientmagnitudesrange fromThto Tl,they are considered as the suspected edge points and their connectivityis examined.If their adjacent pixels have edge pixels,they are also considered as edge pixels;otherwise,they are non-edge pixels.

The Otsu’s method is a nonparametric approach selecting the thresholds in an optimum manner,which makes the separability of between-class variances maximum toautomatically determine the thresholdsToof the original image.The low thresholdTlis got fromTo.To be conservative,the low thresholdTl=τ?To.Hereτis a positive constant as the low threshold scale factor.In our research,τ=0.9.

The example of GMH is shown in Fig.3 with the original image Fig.2(b).

Fig.3 Threshold selection by histogram concavity

In Fig.3,we can see non-edge pixels has a peak and the edge area merged in the right fla shoulder.The high thresholdThis to search the maximum of the concavity in GMH corresponding magnitude indexTh.First, to compute the image GMHh(i)withkbins and normalize to a range from 0 to 1.The point at the largest binP1(x1,y1)(nonzero start bin)and the the minimum point corresponding to the nearby end of the histogramP2(x2,y2)(nonzeroend bin)are chosen to form a straight baselineP1P2.The concavityC(i)of histogram is define as the perpendicular distance between the points of bins histogram and the lineP1P2.If the equation of the baselineP1P2 isAx+By+C=0,and theith entry of the histogram is written ash(xi,yi),then

wherexiis the gradient magnitudegray level,andyiis the frequency of pixels.

The maximum of the concavityC(i)is obtained and the corresponding gradient magnitude indexg(i)is the high thresholdThwhichis usedto extractthe edges.Finally,the adaptive edge detection algorithm is completed by importing dual threshold vector[Tl,Th]into edge points connection processing.In Fig.3,the dual thresholdis[0.15,0.20].

3.4Pseudo edge elimination

Edge detection often produces some unexpected edges and we should select those belonging to crater rims.As Fig.4(a)shows,under the lighting of the solar with an inclination angle,all edges contain true rims,dividedboundary of two contrastive areas,and other terrain edges.In Fig.4(b),there is tow obvious pseudo edges in the middle of the crater area in edge detection results.

Fig.4 Edge selection

Here some constrains are used for the crater edge selection process.

(1)Illumination constrain.The gradient direction of crater rims and pseudo edges are opposite along the lighting source direction(see Fig.4).Assume the lighting source vector isη,the gradient vector of edge pixel isg, then the pixel will be reserved if the following condition is satisfied

whereβis an angle threshold(generally less than 45°).

(2)Curveconstraint.Allofcrateredgesareinanarcuate shape,and edges in the shape line should be deleted.

(3)Minimum number of pixels of edge.Considering the image noises and the reliability of ellipse fittin parameters,the size of crater landmarks cannot be too small. Hencethepixelnumberofeveryedgemust belargerthana threshold.Thethresholdofedgepixelnumberis six pixels.

Based on the foregoing principles,edges containing craterrims can be pickedout effectively.The result of edge selection is shown in Fig.4.

3.5Feature matching based on edges

The objectiveof feature matchingis to correctlyformpairs of edges that are in dark and bright areas belonging to thesame crater.Due to lighting and the bowl-shaped structure of craters,the crater rims show in pairs after the previous processing phases,and the pair of edges belonging to the same crater corresponds to the light area and the shade area,respectively.In the selected edges,there are some edges which are two parts of the same crater rims in shape but are not the same crater.In order to identify the crater,every two edges of the same crater must be matched successfully,and the non-paired edges should be removed. Edge feature matching is to search and combine edge pairs which are one to one corresponding to the same craters.

In edge feature matching phase,we search and couple the pairs of shadow and highlight edges into crater candidates by constructing the CME according to several constraints andrules.Accordingto the definitio of the model, we propose the following geometry dimensions and location constraints for edges matching.

Constraint 1Spatial adjacency.Matching of light and dark pairs is presented nearby.The distance between two centroid points of the edges should not be too close or too far.In this algorithm,the distance must be less than twice thelengthofthe longeredgeandlargerthanthe half-length of the shorter edge.Defin two distance thresholdsDmaxandDmin.The combinationmaysucceedwhen the following expression is satisfied

Constraint 2Light restriction.The edge in the shaded area is only matched by one in the light area.Lighting edge lies on the side of the light rays in the reverse direction,whereas the dark edge in the shadow area lies on the light source side.Considering the connection direction of centroid points in paired edges approximately parallel to the direction of light,given a small angle?between the connection direction and light direction.In our research,?=10°.EC(S)andEC(H)are denoted as the centroid in the shaded area and the light area,respectively.Lis denoted as lighting direction.The angle in direction of matching edges meets the following equation:

Constraint 3Orientation consistency.Considering the region properties of each edge,with a certain sunlight angle,there are the similar orientation among the paired edges.The orientation angle is define as the angle between thex-axis and the major axis that has the same secondmomentsas the region.Lettwo orientationEO(S)andEO(H)indicate a pair of the shaded region and the highlight one,respectively.If the ratio ofEO(S)andEO(H)is within a certain range,they can be matched.Given an upper threshold valueEO(max)and a lower oneEO(min), the following formula should be satisfied

Constraint 4Gray intensity restriction.With a certain intensity of sunlight,the intensity distribution between the two matching edges should be within a certain proportion scope.Defin a gray average of the input image asGμ,an upper and a lower coefficient ascmaxandcmin,respectively.When the following inequality holds,the matching is allowed.

Constraint 5Length similar in size.If coupled edges belong to the same crater,the ratio between the shaded and light edge lengths should be relatively close in size. The shadow edges length and light ones are represented byEP(S)andEP(H).Defin an upper coefficienrmaxand a lowercoefficienrmin.Thematchingis allowedwhenthe following expression is satisfied

Based ontheconstraints,we proposean effectivelyedge matching algorithm.The example of feature matching algorithmis shownin Fig.5.The fl w chart ofthis algorithm follows in Fig.6.

Fig.5 Example of result of feature matching

The input of the algorithm is the set of selected edges feature{EDAtt},and the output is matched edges table,define as{MatchList}=(ES,EH,Tag).Tagis a matched label of the shadow edge ofESwith the highlight edgeEH.Tagis defineas a matching state:Tag=0(unused),Tag=1(matched),andTag=2 (used).ESandEHare the order number in edges feature{EDAtt},together with the corresponding crater feature model CRATER.Initially,{ES}initialized in the order number corresponding to{EDAtt},both all{EH}and{Tag}are set to{0}respectively.A detailed description of the algorithm is shown as follows.

Fig.6 Flow chart of matching algorithm

Algorithm 1Feature matching based on edges

InputFeature table of edges{EDAtt},lighting vectorL,gray meanGyμ,length meanPmμ,orient angle meanAzμ.

OutputMatched edges table{MatchList}with{Tag=1}.

Step 1Initialize{MatchList}={ES,EH,Tag}. SetESrefas a currentreferenceedgefrom the firs item in{MatchList}.Let a search list{SearchList}be the remaining items in{MatchList}.Set a candidate matching table:{EHcands}=?.Set a refin table:{EHrefine}=?.

Step 2Compute the distance between centroid of the edge regionESrefand each item of the search list{SearchList}in order.Add those satisfying the distance constraints into a candidate matching table:{EHcands}.

Step 3If the candidate matching table{EHcands}/=?,lettheitemsequentiallyinEHcandsbeacurrentmatching itemEHcur,then go to Step 4;else go to Step 8.

Step 4LetESrefandEHcurbe a temporary CME feature pair.Four match features are calculated according to(4)–(7),respectively.The three ratios of the currentpair edge are computed to{Gycur,Pmcur,Azcur}.Based on (14)–(17),checktheconstraintsofthecurrentfeaturepair. If all the constraints are satisfied go to Step 5;else go to Step 6.

Step 5AddEHcurto the refin chosenset of the edge candidates{EHrefine}.Compute the similarity factor for confidenc evaluation.Save the similarity factorsmof theEHcurinto a temporary list{Listsm}.

Step 6Let the next item in{EHcands}be indicated asEHcur.Repeat Step 4 for the remaining items in{EHcands}.

Step 7Let the candidate item with the largest similarity factor from{Listsm}correspondingto{EHrefine}as a matched itemEHmatch,and combine it withESrefto a crater feature model CRATER.Mark{Tag=1}and update the current item in{MatchList}with{ESref,EHmatch,Tag=1}.Set the matched item in{MatchList}={EHmatch,EHmatch,Tag=2}.Let{EHrefine}←?.Set the temporary list{Listsm}←?.

Step 8IfESrefis not the last one with{Tag=0}, go to Step 9;else output result and end procedure.

Step 9LetESrefbe the next item with{Tag=0}in{MatchList}and update the{SearchList}with the remaining items in{MatchList}.LetEHcands←?. Repeat Step 2.

3.6Robust ellipse fittin

On the right matched edges,craters parameters are extracted by applying the ellipse fittin algorithm.Although craters are approximatedto an ellipse,eventuallythe crater edges can be fitte to a real ellipse rather than circular.The existing ellipse fittin algorithmscan be classifie into two categories:Hough transform techniques[6]and the least squares methods[20,21].The Hough transform search is an estimation of parameters by clustering in the fve dimensional parameters space due to an ellipse being spec-ifie by fve parameters.Thus,it results in heavy computational load.This paper improves the direct least square fittin [20]method with pre-normalizing edge points to increase the accuracy,and extract the main characters of craters.

On the Cartesian plane,a generic algebraic representation for a conic section is define as a set of edge points Pi=(xi,yi)satisfying the implicit equation in

where h=[a,b,c,d,e,f]Tdenotes the ellipse geometric parameter,M=[x2,xy,y2,x,y,1]with(x,y)represents a pointcoordinatesbelongingto the coniccurve,anda,b,c are not all zero.The conic polynomial could represent an ellipse,hyperbola or parabola decided by the discriminant equation 4ac?b2.For the case of an ellipse,

For general conic least square fittin approach,if a series of point set(xi,yi)(i=1,2,...,n)is acquired,the least square error method can be used to estimate the parameters of the conic section.We firs gather data in the form of a design matrix M which is of the form as

For the case of an ellipse,the inequality constraint 4ac?b2＞0 needs to be imposed.Since the inequality constraint is difficul to implement and there is the freedom coefficient in the h parameters of the general ellipse equation,we may simply incorporate the scaling into the constraint and impose the equality constraint as follows.

This quadratic constraint can be expressed as a matrix form:

The ellipse-specifi fittin approach is transformed into solving the constrained minimization problem through general eigenvalue decomposition.

then,

The exactly one elliptical solution h is the generalized eigenvector corresponding to the single positive generalized eigenvalue of MTM?λK.

Comparedwith the ordinaryleast square(OLS)method, the constraint 4ac?b2=1 can be introduced to ensure an elliptical fitting We use a linear normalization to preprocess the edges data to improve the fittin accuracy.Direct least square ellipse fittin based on normalization is proposed in this paper to increase the robustness in high noise data.By normalizingdata,the edge points are translated so that their centroid is at the origin to improve the stability of solving the scatter matrix MTM.

The following description gives the details of four steps to estimate the ellipse parameters.

(i)To normalize the points(x1,y1),(x2,y2),...,(xn, yn)with min-max normalization.Its purpose is to minimize the cutoff errors in least square fittin calculation.

(ii)To construct the new data matrix M with the row vector,where?is the Hadamard(entry wise)product operator.

(iii)To make the scatter matrix W and the constraint matrix K(see in(20)).They are both of 6×6 dimension.

(iv)To obtain the estimation of the six parameters of the fittin ellipse by calculating the eigenvalues and eigenvectors with Wh=λKh,and then fin the eigenvector h corresponding to the positive eigenvalue λ,which is the estimation of the vector h=[a,b,c,d,e,f]T.The eigenvalues and eigenvector are readily obtained by solving

4.Experimental results

We applyour proposedsaliency detection and edge matching method to Mars and lunar images as representative examples.Mars images are captured by The Mars orbitercamera(MOC)aboardtheMarsglobalsurveyorprobe.Lunar images are taken from lunar orbiter(USA).

Fig.7 presents the detection results on Mars images (R0200575).The illustrations of each step and the fina detection results are given.Observed from the figures the salient craters map after image signature processing brings together at the regions where the surface is of low contrast or complex terrain background as shown in Fig.7(b).The saliency map includes not only a few craters with sharp intensity contrast but also smaller ones with weak intensity contrast.Fig.7(c)shows thatoutstandingedgesare extracted after Canny operator,including any false detected and inconspicuous ones.Many spurious edges in Fig.7(c) disappear in Fig.7(d)by selecting with illumination direction and edges gradient rules.The matched candidate edges are shown in the Fig.7(e)in CME model and some constraints.The fina detected craters together with ellipse fittin are given in Fig.7(f).

Fig.7 Craters extracting main steps

From Fig.8,our saliency detection based method obtains efficien results on lunar images with low and higher solar elevation angles.Referring to the original images Fig.8(a),it can be found that our method is able to detect the saliency edges of lunar images precisely.The results not only contain those salient craters,but also include the inconspicuous or smaller ones,which is similar to human attention.Observed from the figures the two original images do not appear to be very similar,because the sunlight direction and the gray intensity contrast are quite different.However,the CME model based edges matching algorithm works effectively for both the two types of image.The mainreasonis thatCME has consideredthepriori knowledge about craters and the advanced algorithm.

The results of comparisons with our method against the the state of the art method[4]can be seen in Table 1 on four chosen images.From the table,the average accuracy of extractionis 90.4%and is increased by 7.9%,which is a satisfactory detection ratio,considering the differing characteristics of the selected images.

Fig.8 Experiment results of crater extraction

Table 1 Accuracy rate comparison of different methods

5.Conclusions

This paper presents a method based on salient edge information for detecting craters in planetary images.The saliency detection based on image signature provides crater extraction with robustness to different illuminations and terrains.Our methods take both human attention and specifi edges’features of crater into account.The adaptive threshold algorithm provides cater edge detection for different images.We use multi-constraints to delete some pseudo edges in effect for pairing edges of craters.The CME model describes in detail the shape and gray profil features of crater’s attribute.The edges feature matching algorithm based on the CME model shows the effectiveness of pairing edges in the same crater rim.Therefore,the illumination condition tolerance of the proposed algorithm is augmented vastly.Experiments upon the real planetary images show the effectiveness of our method,the crater detection average rate is better than 90%,and can be improvedby 7.9%compared with the state of the art method.

We also fin that the crater extraction has some difficul ties due to a big varietyin backgroundand fluctuan topography.In the future,our main focus is to introduce intelligent methods to remove the complex terrain background for detecting irregular craters.

[1]N.G.Barlow.What we know about Mars from its impact craters.Geological Society of America Bulletin,2010, 122(5/6):644–657.

[2]B.S.Shylaja.Determination of lunar surface ages from crater frequency size distribution.Journal of Earth System Science, 2005,114(6):609–612.

[3]B.Leroy,G.Medioni,E.Jonson,et al.Crater detection for autonomous landing on asteroids.Image and Vision Computing, 2001,19(11):787–792.

[4]C.Yang,A.Johnson,L.Matthies,et al.Optical landmark detection for spacecraft navigation.Proc.of the 13th Annual AAS/AiAA Space Flight Mechanic Meeting,2003:1785–1803.

[5]Y.Sawabe,T.Matsunaga,S.Rokugawa.Automated detection and classificatio of lunar craters using multiple approaches.Advances in Space Research,2006,37(1):21–27.

[6]J.Troglio,J.Le Moigne,J.Benediktsson,et al.Automatic extraction of ellipsoidal features for planetary image registration.IEEE Geoscience and Remote Sensing Letters,2012, 9(1):95–99.

[7]L.Bandeira,J.Saraiva,P.Pina.Impact crater recognition on Mars based on a probability volume created by template matching.IEEE Trans.on Geoscience and Remote Sensing, 2007,45(12):4008–4015.

[8]R.Martins,P.Pina,J.S.Marques,et al.Crater detection by a boosting approach.IEEE Geoscience and Remote Sensing Letters,2009,6(1):127–131.

[9]E.R.Urbach,T.F.Stepinski.Automatic detection of subkm craters in high resolution planetary images.Planetary and Space Science,2009,57(7):880–887.

[10]L.Banderia,W.Ding,T.F.Stepinski.Detection of subkilometer crater in high resolution planetary images using shape and texture features.Advances in Space Research,2012, 49(1):64–74.

[11]B.D.Bue,T.F.Stepinski.Machine detection of Martian impact craters from digitaltopography data.IEEE Trans.on Geoscience and Remote Sensing,2007,45(1):265–274.

[12]Y.Xie,G.Tang,S.Yan,et al.Crater detection using the morphological characteristics of Chang’E-1 digital elevation models.IEEE Geoscience and Remote Sensing Letters,2013, 10(4):885–889.

[13]L.Itti,C.Koch,E.Niebur.A model of saliency-based visual attention for rapid scene analysis.IEEE Trans.on Pattern Analysis and Machine Intelligence,1998,20(11):1254–1259.

[14]L.Gao,F.Bi,J.Yang.Visual attention based model for target detection in large-fielimages.Journal of Systems Engineering and Electronics,2011,22(1):150–156.

[15]X.Hou,J.Harel,C.Koch.Image signature:highlighting sparse salient regions.IEEE Trans.on Pattern Analysis and Machine Intelligence,2012,34(1):194–201.

[16]W.Zhang,F.Liu,L.Jiao,et al.SAR image despeckling based on edge detection and nonsubsampled second generation bandelets.Journal of Systems Engineering and Electronics,2009, 20(3):519–526.

[17]J.Li,P.Huang,X.Wang,et al.Image edge detection based on beamlet transform.Journal of Systems Engineering and Electronics,2009,20(1):1–5.

[18]R.Lai,Y.Yang,B.Wang,et al.A quantitative measure based infrared image enhancement algorithm using plateau histogram.Optics Communications,2010,283(1):4283–4288.

[19]X.Xu,S.Xu,L.Jin,et al.Characteristic analysis of Otsu threshold and its applications.Pattern Recognition Letters, 2011,32(7):956–961.

[20]A.W.Fitzgibbon,M.Pilu,R.B.Fisher.Direct least squares fit ting of ellipses.IEEE Trans.on Pattern Analysis and Machine Intelligence,1999,21(5):476–480.

[21]D.S.Barwick.Very fast best-fi circular and elliptical boundaries by chord data.IEEE Trans.on Pattern Analysis and Machine Intelligence,2009,31(6):1147–1152.

Biographies

An Liuwas born in 1968.He received his M.S. degree from Institute of Northwest Light Industry in 1995.He is now a doctor candidate of Tsinghua University.His research interests include image processing,pattern recognition and signal processing.

E-mail:liua09@mails.tsinghua.edu.cn

Donghua Zhouwas born in 1963.He received his Ph.D.degree from Shanghai Jiaotong University in 1990.He is currently a professor of the Department of Automation at Tsinghua University.He is a senior member of the IEEE.His research activities are in the areas of process identification diagnosis, and control.

E-mail:zdh@mail.tsinghua.edu.cn

Lixin Chenwas born in 1966.He received his Ph.D.degree from Beijing University of Aeronautics and Astronautics in 2005.He is currently a senior researcher of Beijing Institute of Special Electronic-mechanical Technology.His research activities are in the areas of information and system engineering.

E-mail:bjclx@126.com

Maoyin Chenwas born in 1974.He received his Ph.D.degree from Shanghai Jiaotong University in 2003.He is now an associate professor in the Department of Automation,Tsinghua University.His research interest is in the area of reliability analysis,fault prognosis,and image processing.

E-mail:myChen@mail.tsinghua.edu.cn

10.1109/JSEE.2015.00141

Manuscript received October 30,2014.

*Corresponding author.

This work was supported by the National Natural Science Foundation of China(61210012).