An Iterative Black Top Hat Transform Algorithm for the Volume Estimation of Lunar Impact Craters

Volume estimation is a fundamental problem in the morphometric study of impact craters. The Top Hat Transform function (TH), a gray-level image processing technique has already been applied to gray-level Digital Elevation Model (DEM) to extract peaks and pits in a nonuniform background. In this study, an updated Black Top Hat Transform function (BTH) was applied to quantify the volume of lunar impact craters on the Moon. We proposed an iterative BTH (IBTH) where the window size and slope factor were linearly increased to extract craters of different sizes, along with a novel application of automatically adjusted threshold to remove noise. Volume was calculated as the sum of the crater depth multiplied by the cell area. When tested against the simulated dataset, IBTH achieved an overall relative accuracy of 95%, in comparison with only 65% for BTH. When applied to the Chang’E DEM and LOLA DEM, IBTH not only minimized the relative error of the total volume estimates, but also revealed the detailed spatial distribution of the crater depth. Therefore, the highly automated IBTH algorithm with few input parameters is ideally suited for estimating the volume of craters on the Moon on a global scale, which is important for understanding the early processes of impact erosion.


Introduction
The morphologic transition of impact craters with increasing size from simple to complex is the result of evolutions involving the interaction of depositional and erosional processes.The crater volume, one of the most important morphometric parameters of a crater, is a critical piece of information to better understand these processes.However, estimates of cavity volume for lunar craters are difficult to derive due to the lack of a wealth of good data, few direct measurements, and the presence of complex crater shapes.The appearance of central peaks, terraces, an unknown true rim height above the pre-impact surface, the role of structural uplift, and the degree of crater infilling by fallback ejecta all make it difficult to calculate crater volume [1][2][3][4][5][6][7].Use of lunar morphometric relationships for crater depth, rim height, and floor width versus diameter permit approximations of the crater volume.Schröter [8] measured volumes of large complex lunar craters and found that, on average, the volume of the crater rim above the original ground surface V e equaled the volume of the crater void below the original surface V i .Later, several studies found a broad range in the ratio of V e /V i for fresh craters and suggested a significantly higher ratio for large complex craters based on the fact that the accuracy of the data used was not sufficient to yield explicit results [7,9,10].Croft [11] found that craters smaller than ~12 km in diameter were morphologically simple, and their volume increased with the cube of the diameter, while craters larger than ~20 km were complex and their volume increased at a significantly lower rate.A number of models also have been developed to describe the rate and depth of the regolith gardening process that combines probabilistic calculations with laboratory cratering experiments in fine-grained unconsolidated targets.This data can be used to convert projectile kinetic energy into realistic estimates of crater volumes, diameters, and depth [12][13][14].O'Keefe and Ahrens [15] calculated the shock pressures for representative impact conditions that formed the basis for estimating the fractional crater volumes subjected to specific shock stresses.Sundararajan [16] proposed an empirical equation relating the volume of the crater (formed during high velocity oblique impact tests) to the velocity and angle of impact, and to the target material hardness.However, these limited efforts could not accurately characterize the volume of the crater accurately, since indirect calculation and poor data quality were not robust enough to eliminate subjective errors.
With the increasing availability of digital elevation model (DEM) data for the Moon, we are now in a better position to update the estimated volume of excavation by impact craters.This paper demonstrates the application of an improved Top Hat Transform method (TH) into a high-precision DEM as a relevant tool for estimating the cavity volume of lunar craters, which are dim image regions on the lunar surface.The Top Hat Transform function is a powerful image analysis technique developed from mathematical morphology, which is a theory for the analysis of the shape and form of spatial structures initiated in the late sixties by Matheron and Serra [17,18].TH allows the extraction of peaks and pits in a 1D signal and 2D image and has been successfully applied to DEM data to extract valley depth at the pixel level, thus offering more accurate estimates of volumes [19,20].TH works by differentiating current DEM and pre-carving DEM to measure the crater depth and derive crater volume.As one of the two basic functions of TH, Black Top Hat Transform function (BTH) extracts craters by comparing a closing function result to the current DEM.The closing function performs a dilation operation followed by an erosion operation to join clear zones and separate dark zones of the image (Figure 1).If the difference in elevation between the original image and the closed image is above threshold (t), it is flagged as a crater part.The threshold value (t) is used to specifically extract the crater without extracting the artificial values, which guarantees that only real craters are retained.The depth at each pixel multiplied by the pixel area is the volume of that pixel [21,22].The volume of the crater is then computed as the sum of each pixel's volume within the crater according to the following equation: where a i is the area of pixel i, and h i is the depth of pixel i [19,23].
However, traditional BTH methods that only use one size-fixed window and one corresponding threshold value to extract objects often miss the chance of catching craters not included in certain shape and size ranges.To improve the image processing ability of the classical BTH method, an updated BTH was proposed.This method was defined by importing another structuring element of slope to differentiate the region of interest and the surrounding regions.Additionally, we employed an iterative window size and slope that increased slowly and linearly per iteration up to the maximum value that the user specified to improve the estimation accuracy.
In this paper, we name the proposed method as the Iterative Black Top Hat Transform (IBTH) and processed both window sizes and slope as iterative variables, which greatly improved the performance of BTH while suppressing the effect of the noises and utilizing the slope difference between the target region and the surrounding region well.In the following section, we describe how the algorithm works.

Materials
The United States Geological Survey (USGS) lunar geologic mapping program divides the Moon into 30 quadrangles [24].To test our algorithm, we selected a test area 6.09 × 10 11 km 2 located in LQ-27, which includes various types of impact craters (Figure 2).The DEM data from Chang'E-1 was determined from digital photogrammetric image data from a stereo charge-coupled device (CCD) camera on Chang'E-1.In the Mercator projection, each of these pixels is 869 m in size at the equator [25].The lunar orbiter laser altimeter (LOLA) DEM was constructed using the lunar orbiter laser altimeter, an instrument mounted on the NASA lunar reconnaissance orbiter spacecraft.More than 6.5 billion LOLA measurements were converted into a DEM using generic mapping tool software at a resolution of 256 pixels per degree [26].In the Mercator projection, each of these pixels is 118 m in size at the equator.

Methods
This Iterative Black Top Hat (IBTH) algorithm consisted of three notionally specific stages shown in a schematic diagram of overall organization in Figure 3.The first was the basic dilation and erosion operation with window radius and slope calculation based on the closing operation result.The second was the iterative BTH algorithm processing with window size increment ∆r and slope increment ∆s.The third was the post-processing of IBTH to remove false extraction and small isolated regions to outline real craters and calculate their volumes.
The IBTH algorithm requires four parameters: the cell size of the current DEM grid, a percent slope value that governs the grid cell classification at each step, a vector of window radius that controls the BTH operation at each iteration, and a progressive threshold value that governs the ultimate classification of the DEM cell as a crater.The slope (s) was calculated in percent rise and the window radius (r) measured in number of cells, which was based on the cell size and the largest expected crater to be recognized.The initial window radius (r0) and slope (s0) corresponded to the minimum crater size and their increments were iterated to maximum crater size with a linear increment.The IBTH algorithm was defined as follows: (1) executing a dilation operation to find the maximum elevation from the original DEM within the moving window; (2) executing an erosion operation to find the minimum elevation from the dilation result; and (3) executing a closing function to subtract the erosion result from the original DEM.Next, the threshold (t) for removing fake values was calculated as t = r*s* cell size.The window radius was increased by an increment ∆r, as was the

Materials
The United States Geological Survey (USGS) lunar geologic mapping program divides the Moon into 30 quadrangles [24].To test our algorithm, we selected a test area 6.09 × 10 11 km 2 located in LQ-27, which includes various types of impact craters (Figure 2).The DEM data from Chang'E-1 was determined from digital photogrammetric image data from a stereo charge-coupled device (CCD) camera on Chang'E-1.In the Mercator projection, each of these pixels is 869 m in size at the equator [25].The lunar orbiter laser altimeter (LOLA) DEM was constructed using the lunar orbiter laser altimeter, an instrument mounted on the NASA lunar reconnaissance orbiter spacecraft.More than 6.5 billion LOLA measurements were converted into a DEM using generic mapping tool software at a resolution of 256 pixels per degree [26].In the Mercator projection, each of these pixels is 118 m in size at the equator.

Methods
This Iterative Black Top Hat (IBTH) algorithm consisted of three notionally specific stages shown in a schematic diagram of overall organization in Figure 3.The first was the basic dilation and erosion operation with window radius and slope calculation based on the closing operation result.The second was the iterative BTH algorithm processing with window size increment ∆r and slope increment ∆s.The third was the post-processing of IBTH to remove false extraction and small isolated regions to outline real craters and calculate their volumes.
The IBTH algorithm requires four parameters: the cell size of the current DEM grid, a percent slope value that governs the grid cell classification at each step, a vector of window radius that controls the BTH operation at each iteration, and a progressive threshold value that governs the ultimate classification of the DEM cell as a crater.The slope (s) was calculated in percent rise and the window radius (r) measured in number of cells, which was based on the cell size and the largest expected crater to be recognized.The initial window radius (r 0 ) and slope (s 0 ) corresponded to the minimum crater size and their increments were iterated to maximum crater size with a linear increment.The IBTH algorithm was defined as follows: (1) executing a dilation operation to find the maximum elevation from the original DEM within the moving window; (2) executing an erosion operation to find the minimum elevation from the dilation result; and (3) executing a closing function to subtract the erosion result from the original DEM.Next, the threshold (t) for removing fake values was calculated as t = r*s* cell size.The window radius was increased by an increment ∆r, as was the slope value with increment ∆s, until they reached the maximum (r max ) and (s max ).Based on mutative (t), we defined (bthmask*) as a mask grid received value of 1 when BTH value was greater than (t) and (ibth*) as the value grid retained its own value when its value was greater than (t).We merged (bthmask*) as (ibthmask) and selected the mean (ibth*) as a candidate to calculate crater volume.On the merged (ibthmask), the isolated and small islands that were not real crater parts were filtered automatically by using the mapped crater rim polygon [27] as the threshold layer to identify crater areas as (finibthmask).The final ibth was the mean (ibth*) selected by the (finibthmask) to calculate the crater volume according to Equation ( 1) by summing up all values of the cell volume.The above procedure was carried out with Python in Arcpy in ArcGIS Desktop 10.
Remote Sens. 2017, 9, 952 4 of 10 slope value with increment ∆s, until they reached the maximum (rmax) and (smax).Based on mutative (t), we defined (bthmask*) as a mask grid received value of 1 when BTH value was greater than (t) and (ibth*) as the value grid retained its own value when its value was greater than (t).We merged (bthmask*) as (ibthmask) and selected the mean (ibth*) as a candidate to calculate crater volume.On the merged (ibthmask), the isolated and small islands that were not real crater parts were filtered automatically by using the mapped crater rim polygon [27] as the threshold layer to identify crater areas as (finibthmask).The final ibth was the mean (ibth*) selected by the (finibthmask) to calculate the crater volume according to Equation ( 1) by summing up all values of the cell volume.The above procedure was carried out with Python in Arcpy in ArcGIS Desktop 10.

Results
Algorithms are typically tested against computer-simulated datasets for which the true ground is known, or are tested against available results [28].In this case, we first tested the BTH and IBTH slope value with increment ∆s, until they reached the maximum (rmax) and (smax).Based on mutative (t), we defined (bthmask*) as a mask grid received value of 1 when BTH value was greater than (t) and (ibth*) as the value grid retained its own value when its value was greater than (t).We merged (bthmask*) as (ibthmask) and selected the mean (ibth*) as a candidate to calculate crater volume.On the merged (ibthmask), the isolated and small islands that were not real crater parts were filtered automatically by using the mapped crater rim polygon [27] as the threshold layer to identify crater areas as (finibthmask).The final ibth was the mean (ibth*) selected by the (finibthmask) to calculate the crater volume according to Equation ( 1) by summing up all values of the cell volume.The above procedure was carried out with Python in Arcpy in ArcGIS Desktop 10.

Results
Algorithms are typically tested against computer-simulated datasets for which the true ground is known, or are tested against available results [28].In this case, we first tested the BTH and IBTH

Results
Algorithms are typically tested against computer-simulated datasets for which the true ground is known, or are tested against available results [28].In this case, we first tested the BTH and IBTH algorithms to a simulated surface of the Moon, which was the result of a simulated cratered landscape (Figure 4a), and then applied them to the Chang'E DEM and LOLA DEM.A cratered landscape can be modeled geometrically utilizing the approach used by Sugita and Matsui, Gaskell, and Hartmann and Gaskell [29][30][31], which is based on the power law relationship, where a broad range of crater diameters follow.
The simulation domain is 810 × 997 cells with each cell 10 m and the slope distribution generally matched that of the real Moon landform.The true depths of the craters on this simulated surface can be easily derived by calculating the difference between the initial surface and the final cratered surface (Figure 4b).The volume based on these true depths was 1.28 × 10 10 m 3 (Figure 4c) and was used as the "truth" to measure the accuracy of the volume estimation from the BTH and IBTH algorithms.For the BTH method, we used the window size r = 100 cells and slope factor s = 0.55.For IBTH, window size varied from 5 to 100, slope ranged from 0.35 to 0.55 at 0.05 increments, and thresholds varied from 17.5 to 550 m.In Figure 5, the volume generally increased rapidly with the enlargement of the window size under a certain slope, whereas volume decreased slowly if the slope was enlarged under a certain window size.The largest volume appeared at a window size of 70, after which the volume did not change too much even if the window size increased.The volume from BTH was 1.97 × 10 10 m 3 and that from IBTH was 1.35 × 10 10 m 3 , which were all similar to the true volume.However, IBTH had a higher overall relative accuracy of 95%, which was much closer to the "truth" than that for BTH at 65%.The spatial distributions of the BTH and IBTH algorithms applied to the simulated surface are shown in Figure 4d,f, respectively.The spatial correlation coefficient [32] between the estimated depth grid and the "true" grid for IBTH (0.73) was better than that for BTH (0.62).The spatial distribution of difference between the BTH estimate and the "truth" at pixel level (Figure 4e) ranged between 239 and 806 m, with a mean of 532 m and a standard deviation of 132 m.Likewise, for IBTH (Figure 4g), however, the spatial distribution of difference ranged between −311 and 256 m, with a mean of −11 m and a standard deviation of 127 m, where the less extreme values and better statistics than those from BTH results suggest that the IBTH estimate has a more precise depth value (at pixel level) than BTH.The ratio of the number of pixels that are in the mapped craters to the number of pixels that are not in the IBTH extracted craters (false negative, FN) was 1.2%, and the ratio of the number of pixels in the IBTH extracted craters to the number of pixels that are not in the mapped craters (false positive, FP) was 3.1%, while the BTH FN = 1.7%, and FP = 3.3%.The smaller values of IBTH FN and FP indicated that IBTH was better at crater recognition (Table 1).IBTH also achieved a higher precision in estimating the depth of overlapped craters.In comparison to BTH, the IBTH method improved the estimating accuracy of 45.9%, which demonstrates the superiority of the IBTH method on the simulated lunar surface.
Remote Sens. 2017, 9, 952 5 of 10 algorithms to a simulated surface of the Moon, which was the result of a simulated cratered landscape (Figure 4a), and then applied them to the Chang'E DEM and LOLA DEM.A cratered landscape can be modeled geometrically utilizing the approach used by Sugita and Matsui, Gaskell, and Hartmann and Gaskell [29][30][31], which is based on the power law relationship, where a broad range of crater diameters follow.The simulation domain is 810 × 997 cells with each cell 10 m and the slope distribution generally matched that of the real Moon landform.The true depths of the craters on this simulated surface can be easily derived by calculating the difference between the initial surface and the final cratered surface (Figure 4b).The volume based on these true depths was 1.28 × 10 10 m 3 (Figure 4c) and was used as the "truth" to measure the accuracy of the volume estimation from the BTH and IBTH algorithms.For the BTH method, we used the window size r = 100 cells and slope factor s = 0.55.For IBTH, window size varied from 5 to 100, slope ranged from 0.35 to 0.55 at 0.05 increments, and thresholds varied from 17.5 to 550 m.In Figure 5, the volume generally increased rapidly with the enlargement of the window size under a certain slope, whereas volume decreased slowly if the slope was enlarged under a certain window size.The largest volume appeared at a window size of 70, after which the volume did not change too much even if the window size increased.The volume from BTH was 1.97 × 10 10 m 3 and that from IBTH was 1.35 × 10 10 m 3 , which were all similar to the true volume.However, IBTH had a higher overall relative accuracy of 95%, which was much closer to the "truth" than that for BTH at 65%.The spatial distributions of the BTH and IBTH algorithms applied to the simulated surface are shown in Figure 4d,f, respectively.The spatial correlation coefficient [32] between the estimated depth grid and the "true" grid for IBTH (0.73) was better than that for BTH (0.62).The spatial distribution of difference between the BTH estimate and the "truth" at pixel level (Figure 4e) ranged between 239 and 806 m, with a mean of 532 m and a standard deviation of 132 m.Likewise, for IBTH (Figure 4g), however, the spatial distribution of difference ranged between -311 and 256 m, with a mean of -11 m and a standard deviation of 127 m, where the less extreme values and better statistics than those from BTH results suggest that the IBTH estimate has a more precise depth value (at pixel level) than BTH.The ratio of the number of pixels that are in the mapped craters to the number of pixels that are not in the IBTH extracted craters (false negative, FN) was 1.2%, and the ratio of the number of pixels in the IBTH extracted craters to the number of pixels that are not in the mapped craters (false positive, FP) was 3.1%, while the BTH FN = 1.7%, and FP = 3.3%.The smaller values of IBTH FN and FP indicated that IBTH was better at crater recognition (Table 1).IBTH also achieved a higher precision in estimating the depth of overlapped craters.In comparison to BTH, the IBTH method improved the estimating accuracy of 45.9%, which demonstrates the superiority of the IBTH method on the simulated lunar surface.After testing the simulated surface, we put both BTH and IBTH methods into use on the Chang'E DEM (with a resolution of 869 m) and LOLA DEM (with a resolution of 118 m) in the southeastern part of the LQ-27 region on the Moon to further delineate the feasibility of the IBTH method.The window radius used for BTH for the Chang'E data was r = 225 cells, and that for LOLA data was r = 1700 cells, which corresponded to the radiuses of the largest craters.The window sizes for the IBTH for Chang'E data were r = 100, 115, 130, . . ., 235 cells, the slope values were s = 0.01, 0.015, 0.02, 0.025, and 0.03 and that for the LOLA data were r = 800, 850, 900, . . ., 1700 cells, s = 0.06, 0.065, 0.07, 0.075, and 0.08.In Table 1, the smaller FN and FP values clearly showed that IBTH detected more pixels as part of the craters than BTH did on both the Chang'E and LOLA data.The pixels that IBTH extracted but BTH did not were mainly distributed in the small craters and the rims of heavily eroded craters (Figure 6).Hence, the crater volume from Chang'E data from BTH was 5.19 × 10 14 m 3 and that from IBTH was 5.37 × 10 14 m 3 , and the LOLA data showed the same calculation results (BTH results were smaller than IBTH results).IBTH also achieved a higher precision in estimating the depth of craters on a pixel level.In comparison to BTH, the IBTH method applied to both datasets improved the estimating accuracy of 29.5%, which demonstrates the superiority of the IBTH method on the real lunar surface (Table 1).After testing the simulated surface, we put both BTH and IBTH methods into use on the Chang'E DEM (with a resolution of 869 m) and LOLA DEM (with a resolution of 118 m) in the southeastern part of the LQ-27 region on the Moon to further delineate the feasibility of the IBTH method.The window radius used for BTH for the Chang'E data was r = 225 cells, and that for LOLA data was r = 1700 cells, which corresponded to the radiuses of the largest craters.The window sizes for the IBTH for Chang'E data were r = 100, 115, 130, …, 235 cells, the slope values were s = 0.01, 0.015, 0.02, 0.025, and 0.03 and that for the LOLA data were r = 800, 850, 900, …, 1700 cells, s = 0.06, 0.065, 0.07, 0.075, and 0.08.In Table 1, the smaller FN and FP values clearly showed that IBTH detected more pixels as part of the craters than BTH did on both the Chang'E and LOLA data.The pixels that IBTH extracted but BTH did not were mainly distributed in the small craters and the rims of heavily eroded craters (Figure 6).Hence, the crater volume from Chang'E data from BTH was 5.19 × 10 14 m 3 and that from IBTH was 5.37 × 10 14 m 3 , and the LOLA data showed the same calculation results (BTH results were smaller than IBTH results).IBTH also achieved a higher precision in estimating the depth of craters on a pixel level.In comparison to BTH, the IBTH method applied to both datasets improved the estimating accuracy of 29.5%, which demonstrates the superiority of the IBTH method on the real lunar surface (Table 1).

Discussion
The IBTH allowed the extraction of craters of a given size depending upon the chosen window size.The maximum window radius corresponded to the largest crater radius and the minimum window radius corresponded to the smallest crater radius.The slope factor corresponded to the background slope of the topography.Furthermore, the application of the threshold operation allowed us to extract craters from the background topography.
The advantage of IBTH was that it estimated crater volume with a higher accuracy by adopting iterative window sizes and slope factors: a series of window sizes covered the different size of craters, which guaranteed that the craters were extracted in full shape.Meanwhile, the IBTH obtained a more precise pixel depth than did BTH since the small window size measured the depth more accurately and the large window size made up for its disadvantage of missing most of the crater.Combining these series of window sizes, IBTH considered accuracy at both pixel level and object (crater) level and obtained much better estimating results than BTH.As for the various slope factors, its key ability was to distinguish crater pixel from the background landscape and eliminate mis-extracted crater pixels.In particular, when the window size increases linearly, the slope factor should also increase correspondingly (Figure 5), and the threshold t = r*s* cell size adjust automatically as the window size and slope factor increases.Besides, the rest of the process in the IBTH method was automated after the window size and slope factor values were determined.With the help of multi-scale window sizes and multiple slope factors, the error resulting from a fixed window size and slope factor was minimized significantly, which leads to a more precise estimation of the total volume.
On the simulated surface, the smaller FN and FP values suggested that the estimate error mainly came from window size and not slope factor since the simulated craters were in a perfect bowl shape and were easily extracted.The BTH method overestimated the total volume by using the largest window size to make sure that the largest crater could be recognized, while ignoring the large window size considered too many neighborhood pixels to estimate the local depth.Thus, the BTH volume result was bigger than the IBTH volume.On the real lunar surface, the bigger FN and FP values implied that the slope factor plays a slightly more important role than window size given that some eroded craters with blurred rims need slope to distinguish them.The IBTH method extracted a greater number of smaller craters and crater rims that the BTH missed, meaning that the BTH method obtained a smaller volume estimate than the IBTH method.The increase in IBTH volume estimation was significant since it was related to the transition from simple to complex crater morphology [11].The higher-resolution LOLA data also provided a higher estimate when compared to the Chang'E data as there were more pixels in the same area with a smaller cell size.The relative error was the difference between the BTH and IBTH divided by the IBTH.The positive relative error indicated that the IBTH worked 29.5% better than the BTH and the same relative error from different datasets also proved the robustness and accuracy of the IBTH (Table 1).Based on the high overall relative accuracy achieved when applying the IBTH method to the simulated data, the IBTH result from the real data not only offers a more precise depth estimate, but also shows the practiced approach of volume calculation.
The IBTH was also easy to use, as it did not require many parameters that other analytical equations or empirical equations do to predict crater volume.Ratner and Styller [33] obtained an analytical expression using impact velocity to predict crater volume but did not fit the experimental data well.Sundararajan [16] had proposed an empirical equation that requires the impact velocity, the angle of impact, and the target material hardness to estimate crater volume.The datasets for these three parameters were not simple to obtain, and there were also certain impact angle limits for applying the above empirical equation.The IBTH, however, only requires window sizes and slope factors as two simple parameters to estimate crater volume; where DEMs were readily available in various resolutions, window sizes and slope factors were easy to determine on the basis of the state of the craters.

Conclusions
The Iterative Black Top Hat Transform algorithm was developed to solve two problems.First, it was arranged to calculate the volume of the crater without knowing the environment of the initial surface, particularly with regard to the surface of planetary bodies.It was successful in that it achieved a high overall accuracy and a lower standard deviation of spatial distribution of difference on the simulated surface.The IBTH algorithm was successful not only when applied to a simulated surface, but even when used on a real lunar surface, suggesting that novice users could achieve good results on bowl-shaped objects on planetary bodies.
The second contribution of IBTH was that it established a baseline performance for a progressive volume estimation implemented in its iterative form.The essence of the IBTH algorithm solved the problem of overestimating the crater depth that BTH achieved with the assistance of two parameters: the window radius (that recognizes the different sizes of craters) and the slope parameter (that governs the cell-based crater rim/ground flagging at each iteration).With these two changeable parameters, the central subroutine of the IBTH considered the characterization of the rim and depth of different crater sizes, instead of using one oversized window and slope to determine the largest crater, while amplifying the pixel value of the smaller craters.The real contribution of IBTH was that it implemented a theoretically and computationally simple basis to achieve the calculation of the crater volume.
By applying the simple and effective Iterative Black Top Hat Transform function with few input parameters to high-precision DEMs, we developed a new method for estimating crater volume on the Moon inspired by a similar method used in valley volume.Changeable window size and slope factor were two significant parameters that control the precision of the volume estimate in this new method.The application of this algorithm to a simulated cratered surface landform achieved an overall relative accuracy of 95% and the crater spatial pattern was generally consistent with the truth, with a spatial correlation coefficient reaching 0.73.The application of the algorithm to the tested area on the Moon also produced robust and accurate results, where higher resolution LOLA data generated greater volume estimates than that from the Chang'E data.
While the IBTH performed well overall, future work will address the application of this method to the entire Moon to estimate the global volume of the crater with the aim of shedding new insights into the processes of impact erosion.The enhanced efficiency of the IBTH function and increasing advent of high-resolution DEM data could also assist in other cases of analysis of mechanical erosion on the Moon.In addition, it will be interesting to apply IBTH to the estimation of crater volume on other planets including craters on the Earth, Mercury, etc.It would also be worthwhile to conduct more research on the uncertainty associated with error in various resolution digital elevation models (DEMs) and the propagation of error in the process of estimating volumes.

Figure 1 .
Figure 1.1D profiles of craters after extraction by the Black Top Hat Transform function (BTH).

Figure 1 .
Figure 1.1D profiles of craters after extraction by the Black Top Hat Transform function (BTH).

Figure 2 .
Figure 2. (Left) Localization of the study area.(Right) Chang'E Digital Elevation Model (DEM) of the study area.Also shown are craters named by the International Astronomical Union (IAU).

Figure 2 .
Figure 2. (Left) Localization of the study area.(Right) Chang'E Digital Elevation Model (DEM) of the study area.Also shown are craters named by the International Astronomical Union (IAU).

Figure 2 .
Figure 2. (Left) Localization of the study area.(Right) Chang'E Digital Elevation Model (DEM) of the study area.Also shown are craters named by the International Astronomical Union (IAU).

Figure 4 .
Figure 4. Simulated test surface and depth estimates.(a) Simulated cratered surface; (b) difference between initial and cratered surfaces; (c) difference surface thresholded at 0 m to extract crater depth; (d) crater depth extracted by BTH; (e) difference between (d) and (c), the arrows point out that the more red the color, the greater the difference, and the less accurate the estimation; (f) crater depth extracted by IBTH; (g) difference between (f) and (c), the arrows pointed that the darker the color, the greater the difference, and the less accurate the estimation.

Figure 4 .
Figure 4. Simulated test surface and depth estimates.(a) Simulated cratered surface; (b) difference between initial and cratered surfaces; (c) difference surface thresholded at 0 m to extract crater depth; (d) crater depth extracted by BTH; (e) difference between (d) and (c), the arrows point out that the more red the color, the greater the difference, and the less accurate the estimation; (f) crater depth extracted by IBTH; (g) difference between (f) and (c), the arrows pointed that the darker the color, the greater the difference, and the less accurate the estimation.

Figure 5 .
Figure 5.The iterative calculating process of IBTH with the changeable window size and slope value.Different color lines with points represent the volume change under different window sizes with a certain slope factor.The purple line is the true volume of all craters on the simulated surface and the orange line is the IBTH volume of simulated craters.

Figure 5 .
Figure 5.The iterative calculating process of IBTH with the changeable window size and slope value.Different color lines with points represent the volume change under different window sizes with a certain slope factor.The purple line is the true volume of all craters on the simulated surface and the orange line is the IBTH volume of simulated craters.

Figure 6 .
Figure 6.Crater depth estimates for the tested area located in LQ-27on the Moon.(a) Crater depth extracted by BTH based on Chang'E DEM.In comparison to the IBTH method, the red boundary of craters indicates the missing crater rim that BTH did not detect and the blue rectangles point out the missing craters that BTH did not recognize; (b) crater depth extracted by IBTH based on Chang'E DEM; (c) crater depth extracted by BTH based on LOLA DEM.In comparison to the IBTH method, the red craters boundary indicates the missing crater rim that BTH did not detect and the blue rectangles point out the missing craters that BTH did not recognize; (d) crater depth extracted by IBTH based on LOLA DEM

Figure 6 .
Figure 6.Crater depth estimates for the tested area located in LQ-27on the Moon.(a) Crater depth extracted by BTH based on Chang'E DEM.In comparison to the IBTH method, the red boundary of craters indicates the missing crater rim that BTH did not detect and the blue rectangles point out the missing craters that BTH did not recognize; (b) crater depth extracted by IBTH based on Chang'E DEM; (c) crater depth extracted by BTH based on LOLA DEM.In comparison to the IBTH method, the red craters boundary indicates the missing crater rim that BTH did not detect and the blue rectangles point out the missing craters that BTH did not recognize; (d) crater depth extracted by IBTH based on LOLA DEM

Table 1 .
Comparison of volume estimate by BTH and IBTH algorithms