The Martian Boulder Automatic Recognition System, MBARS

Abstract

Boulder‐sized clasts are common on the surface of Mars, and many are sufficiently large to be resolved by the high resolution imaging science experiment (HiRISE) camera aboard the Mars reconnaissance orbiter. The size, number, and location of boulders on the surface and their spatial distribution can reveal the processes that have operated on the surface, including boulder erosion, burial, impact excavation, and other mechanisms of boulder transport and generation. However, quantitative analysis of statistically significant boulder populations, which could inform these processes, entails prohibitively laborious manual segmentation, granulometry, and morphometry measurements over large areas. Here, we develop, describe, and validate an automated tool to locate and measure boulders on the Martian surface: the Martian Boulder Automatic Recognition System (MBARS). Our open‐source Python‐based toolkit automatically measures boulder diameter and height in HiRISE images enabling rapid and accurate assessments of boulder populations. We compare our algorithm with existing boulder‐counting methods, manual analyses, and objects of known size to verify accuracy and precision. Additionally, we test how MBARS quantitatively characterizes boulders around an impact crater in the Martian northern lowlands. We compare this to previous work on rock excavation during impact cratering using manually counted boulders around lunar craters.

Introduction

Images taken by the high resolution imaging science experiment (HiRISE) camera aboard the Mars reconnaissance orbiter (MRO) show that meter-scale boulders, blocks, and other megaregolith observed by landers and rovers are common across the entire surface of Mars (Golombek et al., 2008, 2012). Observations of megaregolith on Mars and other bodies have been used to examine a wide variety of surface processes, including impact cratering (Krishna & Senthil Kumar, 2016; Levy et al., 2018; Watkins et al., 2019), bedrock degradation (Nagle-McNaughton et al., 2020), thermal cracking (Eppes et al., 2015), erosion (de Haas et al., 2013; Golombek et al., 2006), and glacial processes (Levy et al., 2021). Landing site assessment is also an essential step in all missions, for which boulders are a major landing hazard (Golombek et al., 2008, 2012; Wu et al., 2022). However, the fundamentally time-consuming and difficult task of manually identifying, locating, and measuring boulders, blocks, and other megaregolith (henceforth simply boulders) is a large burden on these investigations. Tools to automatically locate and measure boulders have been previously developed for Martian landing-site analyses (Golombek et al., 2008, 2012, 2016) as well as more general applications (Nagle-McNaughton et al., 2020) and lunar studies (Li & Wu, 2018) to ease this burden. For Mars, two existing algorithms have substantial challenges to their application. The algorithm developed for landing site analysis (Golombek et al., 2008), henceforth referred to as the Golombek-Huertas (G-H) method, typically requires mission team expertise for accurate use, posing challenges to widespread replicable adoption. The more recent method developed for general use (Nagle-McNaughton et al., 2020), henceforth referred to as the Nagle-McNaughton (N-M) method, does not assess boulder height, and results from that method are given as upper and lower bounds on population morphometry, which are insufficient for some investigations. Here, we advance the existing boulder-measuring approaches with an open-source Python-based methodology to automatically identify and measure boulders in high-resolution satellite images of the Martian surface. The core of this methodology is the newly developed Mars boulder automatic recognition system (MBARS), which detects boulders via shadow identification, an established technique (Golombek et al., 2008), and uses common Python libraries and standard geographic information system (GIS) files and formats.

In this paper, we describe the MBARS algorithm as well as the complete methodology to determine boulder morphometry from HiRISE images. We then test MBARS against objects of the known size, compare MBARS results to other algorithms and manual analyses, and discuss uncertainties and errors within the methodology. As a demonstration of potential applications of MBARS, we examine the boulder population surrounding an unnamed crater in one of the test images, characterizing the ejecta and comparing it to previously analyzed lunar craters of similar size (Watkins et al., 2019).

Methods The HiRISE Data Set

The primary data for this work are images from the HiRISE camera, a high-resolution camera that provides up to ∼25 cm/pixel images of the Martian surface. Typical HiRISE footprints are ∼3 × 6 km and images are globally distributed and not contiguous due to their small size. The point-spread function (PSF) of HiRISE is ∼1.5 pixels (Kirk et al., 2008; McEwen et al., 2007), so objects ≳1 m across are resolvable at typical resolutions. The PSF sets a lower bound on measurements made in HiRISE images, as precision below ±1.5 pixels is not achievable without deconvolution or other image-enhancing methods. Because these image-enhancing methods are not applied here, we assume a ± 37.5 cm uncertainty on all estimated boulder dimensions. Unless otherwise specified, we use the map-projected black-and-white JP2 data products available on the planetary data system (PDS). All images were spatially corrected, when necessary, using the fix_jp2 protocol available on the PDS (https://pdsimage2.wr.usgs.gov/pub/pigpen/).

The PSF and resolution of HiRISE images limit the diameter of resolvable boulders to ≳1 m, but the incidence angle of the Sun in each observation places further limitations. Lower incidence angles (Sun closer to zenith) will cause shadows to shorten, potentially below detection limits, rendering boulders undetectable via their shadows. Using our boulder model (Figure 1), we can explore if there are incidence angles for which we do not expect to detect boulders within any given size range. A diameter of three pixels across is considered a reasonable lower limit to reliably detect an object, so we use a minimum shadow length of three pixels (∼0.75 m) to determine the lower incidence angle limit. The minimum detectable boulder height (from the casting point, H_m) in meters can be calculated as $${h}_{\mathrm{min}}=0.75/\mathrm{tan}(\theta )$$, where θ is the solar incidence angle measured from zenith. On the Moon, boulder height to diameter ratio (h/D) averaged 0.54 across several study sites (Demidov & Basilevsky, 2014) and Martian h/D appears comparable or smaller (Golombek et al., 2012). Taking this ratio as a starting assumption for Mars, the minimum boulder diameter that will cast a sufficiently long shadow is $${D}_{\mathrm{min}}=1.4/\mathrm{tan}(\theta )$$. Due to the Sun-synchronous orbit of MRO, the incidence angle of most HiRISE images falls between ∼40° and ∼75°. In these images, boulders of diameters above 0.37 and 1.67 m, respectively, will cast detectable shadows. Therefore, boulders that are sufficiently wide to be detected (≳1 m) will also cast detectable shadows in most HiRISE images, and boulders >1.5 m wide are predicted to cast detectable shadows in nearly all HiRISE images if they are present.

Figure 1. A geometric model of the map-view and side view of an observed boulder in a high resolution imaging science experiment image. The boulder is modeled as a spheroid oriented vertically (i.e., with one axis oriented normal to the surface) on level ground, which casts a shadow expressed as an ellipse with a horizontal axis equivalent to the boulder width and a vertical axis equivalent to twice the length of the shadow. The sunward-most extent of the boulder shadow is used as the lateral center of the boulder, and the shadow boundary (black stars) is geometrically mirrored about this axis to create the mirrored points (red stars). The ellipse defined by the black and red points defines the ideal shadow cast by a model (spheroidal) boulder. Equations 3 and 4 describe the relationships between the measured height (Hm), the actual height (Ha), shadow length (Ls), boulder radius (r), and the Sun incidence angle (θ). The relationship between Hm and Ls (Equation 1) is straightforward; however, determining Ha from Hm is more complex.

Manual Boulder Analysis

MBARS produces several output data sets based on changing parameters (Section 3.1), and it is difficult to determine which of these running parameters is correct a priori. The simplest way to determine the correct value for this shadow boundary is to compare the automated output to manually measured boulders. To perform this calibration, we manually analyze several areas in each image, characterizing the boulder population within the area, and choose the MBARS boundary parameter to maximize the match between automated results and manual results.

In each image, we manually outline the boulders in two small areas (∼100 × ∼100 m), which we use as test areas for calibration. These test areas are manually selected for each image based on several qualities. (a) A lack of major nonboulder objects, (b) low relief and lacking obvious shadow-casting topographic features, (c) representative areas of the different boulder densities in the image (i.e., moderate and high density), and (d) areas that capture the range of surface albedo within the image as the contrast between shadowed areas and soil is an important factor in detecting boulders via their shadows. In cases where an image displays strong dichotomies in soil albedo (with substantial boulder populations in each), more test areas could be added to accommodate this diversity. We outline boulders within the test areas in Esri's ArcGIS using the create polygons feature. We use circular polygons to represent each boulder to match better with MBARS, which assumes a spheroidal boulder shape that is circular in map view. In the event that a calibration area with sufficient (≳10) abundance of detectable boulders can be found (this was the case in images in the vicinity of VL1, Section 4.1), random boulders sampled from the whole image can be manually measured, and MBARS results can be chosen to best match the random sample.

Boulders are measured with varying techniques in different works (e.g., Krishna & Senthil Kumar, 2016; Watkins et al., 2019). We sought to test how using these existing data sets and therefore a different boulder measurement technique to calibrate MBARS impacted the results. Prior work examining boulders in Acidalia Planitia (Sholes et al., 2017) is one such example of an existing work, which produced manual boulder measurements using a different methodology. In contrast to MBARS, which measures boulder shadows directly, boulders in image PSP_007718_2350 were measured using the illuminated boulder face, using shadows only to confirm boulder presence. Three areas (1, 2, and 3, each 10⁴ m², Figure S1 in Supporting Information S1) were randomly selected from the image, have no large topographic features or craters, and are dominated by their boulder populations. All identifiable boulders were manually measured in the three subsets using the CraterHelperTools plugin for ArcGIS (Kneissl et al., 2011). For small boulders, the three-point circle mode was used to create boulder polygons with the three points selected where the best image contrast was present. Larger irregular boulders were measured using the six-point ellipse tool. Adding to these measurements, we also defined and manually measured boulders in another 10⁴ m² area, Area A, using our circular boulder methodology. We compare MBARS results to all four test areas to determine whether the two manual methods produce compatible results, and whether new manual analyses are required where boulder surveys already exist.

In each test area, we calculate the Rock Abundance (RA) based on the manual and automated observations. RA refers to the total surface area that is covered by rocks of all sizes (k), including rocks smaller than those detectable in HiRISE images. We determine RA in a given area by a fit to the equation (Golombek et al., 2012): [Image Omitted. See PDF]where F_k(D) is the cumulative fractional area (CFA) covered by rocks of diameter > D in m, k is the RA (expressed as a fraction), and the exponential factor q(k) is given by [Image Omitted. See PDF]

The equation for q was empirically determined in previous work (Golombek & Rapp, 1997) and has shown good agreement with Martian surface RAs at several scales (Golombek et al., 2012). It is worth noting that values of k can exceed 1, that is, >100% RA, especially in small study areas, occurring due to an excess of detectable (>1.5 m) boulders relative to expected distributions. Fits between the theoretical CFA and observed CFA are limited due to two competing factors: the inability to resolve boulders ≲1.5 m and a lack of larger boulders in the small test area. These factors cause apparent deviations between observed and theoretical CFA distributions at low and high boulder diameters. However, both these effects are related to the data rather than the technique, so misfits between theoretical and observed CFA are expected to be similar in both manual and automated analyses.

Detection by Shadow Segmentation

Determination of boulder shape from shadows relies upon an assumed basic shape of the boulder. Our idealized boulder model is that of a half-buried spheroid on level ground with a circular shape in plan view and an elliptical cross section when viewed from the ground level (Figure 1). Throughout this work, directions and axes in the images are referenced in relation to the Sun: sunward, anti-sunward, and sun-perpendicular. The algorithm reorients images based on Sun orientation, so typical geographic directions and geometric terminology (minor and major axes) are less useful. The cast shadow of a vertical, spheroidal boulder is a portion of an ellipse with its sun-perpendicular axis equal to the boulder width, and the sunward axis length is determined by the boulder height, boulder width, and incidence angle (Figure 1). Using this shadow model, as opposed to measuring the sun-perpendicular and sunward shadow dimensions directly, provides shape constraints on the shadows and reduces the occurrences of misidentification of low ridges or data errors as boulders. The shadow length is specifically related to the height on the boulder surface where the shadow is cast, labeled in Figure 1 as the “casting point.” Using the shadow length and the incidence angle of the Sun, the measured height (H_m) can be determined as [Image Omitted. See PDF]where L_s is the length of the shadow and θ is the incidence angle of the Sun, measured from zenith. Prior automated methods (Golombek et al., 2008) used a cylindrical boulder model, but the difference between these two is somewhat semantic, as their boulder width was also equal to shadow width, and boulder height measured from shadow length. One advantage of assuming a spheroidal shape is that we can infer the actual height (H_a) of the boulder above the casting point (Figure 1). H_m will consistently be an underestimate of the actual boulder height and the ratio of H_m to H_a can be determined as [Image Omitted. See PDF]where θ is the incidence angle as in Equation 3 and r is the radius of the boulder. If the boulder is assumed to follow lunar boulders (i.e., H_a/D = 0.5 or r = H_a), this underestimate will range from H_m/H_a = 0.5–0.86 for θ = 30°–60° with a greater underestimate at lower θ. Correcting this systematic uncertainty is computationally costly and makes strong assumptions about the boulder shape. The correction can be used as an optional modification to the data, as described below for the Viking lander site estimations, or the MBARS boulder sizes can be treated as minimum heights. Surface slope, and its impact on shadow lengths, is also not accounted for in our approach in part due to the relatively sparse HiRISE stereo coverage. The impact of surface slope on the interpreted height is variable as it depends on both the incidence angle as well as the surface aspect. For surface slopes <10°, shadow lengths are extended (or truncated) by less than 20% of their length, and the effect is most pronounced when shadows are cast downslope and at high incidence angles. The test images used here have slopes less than 2° at the lateral scale of Mars orbiter laser altimeter data (∼500 m), so slope corrections are significantly smaller than PSF-induced uncertainty. Where sufficiently high-resolution digital elevation models (DEMs) are available, a correction to the interpreted heights based on incidence angle, azimuth angle, surface aspect, and surface slope could be derived independent of MBARS processing steps.

A further cause of uncertainty in detecting boulder shadows are shadow penumbrae, the transition from shadowed areas to fully lit areas. Atmospheric scattering is a major cause of penumbrae on Earth, but the thin Martian atmosphere (∼6 mbar) and scattering cross section of CO₂ for HiRISE RED filter wavelengths (∼1.214*10⁻³⁸ m² at 694 nm based on Mie Scattering (Hecht, 2002)) lead to expected scattering roughly 10 orders of magnitude smaller than the green light on Earth. The largest likely source of penumbrae on Mars is the 0.33° angular size of the Sun in the sky (Golombek et al., 2008), which creates penumbrae that are on the order of 3% of the shadow length. For boulders between 1.5 and 5 m wide, this penumbra is ∼10 cm, which is smaller than the 1.5 pixel (∼35 cm) PSF. For our purposes, we assume that the penumbra caused by scattering and Sun size are not detectable within the assumed 1.5 pixel PSF of the HiRISE instrument and consider the PSF as the only source of shadow blurring.

Algorithm Description

The workflow to apply MBARS involves several steps, including preparation, manual analysis, image analysis, and steps, to determine the best-fit solution. These steps are described in detail here as well as in the flowchart in Figure 2.

Figure 2. A flowchart of the Martian Boulder Automatic Recognition System (MBARS) algorithm. The blue box and arrows show the required input data, all of which are available on the planetary data system. The image is specified manually, but the image metadata are automatically fetched by MBARS locally. The green box shows the processing steps of MBARS, described in detail in Section 3. The algorithm generally works with a single image panel at a time and carries out analyses with different shadow boundaries in a sequence. The final outputs (orange box and arrows), the geographic information system (GIS)-ready list of boulder objects and attributes, compile all image panels into a single file for import to a GIS software. The user compares the results of the different boundary parameter settings to the manually counted boulders within the test areas (purple box) choosing the best-fit solution as the final MBARS output.

Image Preparation and Shadow Boundary Selection

Prior to MBARS processing, several steps are taken to prepare the image for analysis. First, the HiRISE image is partitioned into panels, which we carry out with the Split Raster tool in ArcMap. This provides crucial map orientation information to MBARS (as .pgw world files, which are required for later steps) and splits the HiRISE image into manageable sizes for processing. The size of each panel can be controlled by the user though 500–1,000 pixel square image panels (∼125–250 m) are used in this work. The individual panel size is only limited by computer hardware and MBARS can receive panels of any size.

Previous shadow segmentation methods have relied on maximum entropy thresholding (Golombek et al., 2008) or range filtering (Nagle-McNaughton et al., 2020) to define shadow boundaries. The maximum entropy thresholding approach used in the G-H method creates two classes within the image, shadows and nonshadows, and modifies the brightness boundary between those two classes to maximize the interclass entropy. MBARS takes a different approach, predicting a shadow boundary, that is, the expected maximum brightness to be considered a part of a shadow, for each image by forward modeling (Figure 3). To predict a shadow boundary, three key components are used: the darkness of shadows, the brightness of nonshadow pixels, and the PSF of the HiRISE instrument. The interior of shadows larger than ∼5 pixels is not perfectly dark but is generally dark enough to register as digital number (DN, corresponds to pixel brightness in the image) = 1 in HiRISE images. (Note: this is the minimum DN value in 16-bit HiRISE images, as 0 is reserved for NoData.) The brightness of nonshadow pixels varies among images due to changing surface properties (e.g., albedo) and photometric conditions. Instead of predicting these changes, the brightness distribution in the target HiRISE image is statistically sampled and used in the shadow model. Note that this calculation is done on the entire image, not on individual image panels, making the shadow boundary calculation consistent across image panels. Finally, the HiRISE PSF is well-quantified from in-flight imaging of onboard targets and stars (McEwen et al., 2007), though other factors (spacecraft jitter, atmospheric conditions, etc.) are more difficult to constrain. Following previous HiRISE work (Kirk et al., 2008), a Lorentzian function with λ (half-width at half-maximum) = 0.77 is used here for the HiRISE PSF. Other factors that may influence the effective PSF are assumed to be accounted for within this PSF. To predict how a dark shadow will be blurred with the surrounding, brighter nonshadow pixels, we construct a model image with a dark (DN = 1) shadow and convolve it with a background (nonshadowed area, Figure 3) which is randomly sampled from the image pixel brightness distribution. After convolution, the shadow interior becomes brighter due to blurring with the nearby background. The DN of pixels within the constructed shadow are recorded, and this process is repeated 100 times for each HiRISE image. From the 100 results, the average of a user-chosen percentile of each shadow is taken as the shadow boundary. The choice of boundary parameter is the only point of user influence. User selected boundary parameters between 40 and 70 (i.e., 40th and 70th percentiles) produce MBARS results consistent with manual observations (Table S1 in Supporting Information S1).

Figure 3. An example of the shadow convolution method applied to test data. The shadow area is fixed at a value of 1 and the rest of the image is randomly sampled from the image population. Convolution with the high resolution imaging science experiment (HiRISE) Point Spread Function (PSF) yields the theoretical observation by HiRISE. After convolution, the original shadow area, marked in black, is examined to understand the post-convolution brightness within the shadow.

MBARS segments each image panel based on the shadow boundary DN value, below which the pixels are considered to be part of a shadowed area. During this step, MBARS also retrieves relevant metadata (sub-solar latitude/longitude, resolution, incidence angle, etc.) from the RDRCUMINDEX file provided by the PDS. The image is also rotated according to the Sun direction calculated from the sub-solar and sub-spacecraft coordinates provided in the HiRISE image metadata. This first collection of functions results in one primary product: the original image rotated and filtered such that pixel intensities above the shadow boundary are set to a fixed value. This segmented image is passed onto the next major function, Boulder Segmentation.

Boulder Segmentation

The second major step is to select and define individual shadows in the segmented image created in the first step, Image Preparation and Boundary Selection. We use a watershed algorithm to identify individual shadows, treating the dark shadows as basins and the relatively bright exteriors as plateaus and separating shadows that touch one another (e.g., Figure 4d; Figure S2 in Supporting Information S1). Local minima in the DN within each boulder shadow are used as seeds (i.e., starting points) for the watershed algorithm. This process readily identifies individual shadows and divides multiple shadows that form a merged shadow region. This process can create some false splitting in shadows with multiple minima which is corrected later in the process. The watershed process also creates unique IDs (called flags) for each shadow. In this step, shadows that are too small (<4 pixel area) to reliably identify, or too large to likely be boulders (>30 m, i.e., very coarse blocks or other megagravel, (Blair & McPherson, 1999)) are disregarded. With each shadow identified and flagged, MBARS begins the process of determining boulder morphometry based on the shadow dimensions.

Measuring Boulder Dimensions

Once a dark area is identified and defined as an individual shadow, MBARS constructs a boulder object (an internally defined Python data structure) and begins the morphometric analysis portion. The G-H shadow delineation method (Golombek et al., 2008) uses an area-preserving fit, that is, the resultant shadow ellipse had the same area as the initial shadow. We find that the area-preserving constraint did not improve fits and often led to less convergent solutions in our method. Instead of using the area preserving fit, we developed a new method based on the expected shape of the shadow according to our boulder model (Figure 1). As described in Section 2.3, a model boulder shadow is a portion of an ellipse extending from the base of the boulder, but this ellipse portion can be reflected to create a full ellipse as shown in Figure 1. This can be done on shadow observations in the same way, reflecting the anti-sunward shadow boundary across the sunward end of the shadow recreating the full ellipse. The algorithm then, for each shadow, fits an ellipse to the combined original and reflected shadow boundaries (Figures 1 and 4) using orthogonal distance regression (ODR) based on the Levenberg-Marquardt procedure (Boggs et al., 1992). In practice, this technique was more effective than fitting an ellipse to a shadow directly or attempting to fit a half-ellipse to the shadow. The fit is not fixed to align with the sunward and sun-perpendicular axes, as tests undertaken enforcing these constraints caused poorer fits. However, without constraining the orientation, ellipses are closely aligned (<1° difference) with those axes. With the shadow ellipse described, the boulder diameter and boulder height are measured from the sun-perpendicular and sun-parallel axes of the shadow ellipse (Figures 1 and 4). Fits that return diameters or heights that are too large to likely be a boulder (i.e., >30 m) are not discarded but are marked as bad results. Bad results typically come from shadow casting topography (e.g., craters, knobs) or very small objects that produce shadows that are not well described by an ellipse. These boulder objects are retained in order to help distinguish when an object was not detected and when it was mismeasured. At this stage, a list of shadow objects, each containing information about the shadow, shadow ellipse, and boulder, is recorded for each image panel.

Figure 4. Portion of TRA_000828_2495 with (a) no annotations, (b) Manual Results, (c) Nagle-McNaughton (N-M) Method results (yellow = high, purple = low), and (d) Martian Boulder Automatic Recognition System (MBARS) Results. This area exemplifies the anomalous results from the N-M method, which was more sensitive to boulder shadows and not boulders. An example of adjacent boulders that were correctly identified as separate boulders by MBARS is visible in the lower left of the image.

After an image panel is complete, the measured boulders are revisited by MBARS, looking for overlapping shadows. First, all boulder objects with shadows that touch (i.e., have edge-sharing pixels) are considered as possible candidates for merging and are compiled into sets of boulder objects with touching shadows. These groups are typically chains or clusters of adjacent boulders or single boulders that have been incorrectly split by the watershed algorithm. MBARS first merges shadows of all objects in the group and recalculates the boulder object, following all the steps described above. If the shadow was fragmented incorrectly, this new, merged shadow should provide the best solution. Then, the algorithm uses k-means clustering to break the merged shadow into smaller shadows (k = 2, 3,…,n) up to n, the number of originally detected boulder objects in the group. If the group of shadows represents some number of clustered boulders, we expect the k-means solution where k is equal to the number of boulders to provide the best solution. For each potential solution, the fit error (provided by the ODR) for each boulder in the group is calculated and summed. We use the sum, as opposed to the average, of the fit error to bias the method toward fewer boulders. Without this bias, over splitting via k-means would likely lead to overfit shadow solutions and most shadow objects being broken down into the smallest possible components. The solution with the lowest total error is taken as final and saved.

Lists of all boulder objects in an image panel are stored as .shad files, a custom file type created using the CPickle module of Python that serve as the initial record of boulders identified in the image. A second program, MBARS_ANALYSIS, has several built-in tools to analyze the boulder morphometry results from these original files. At the end of each run, MBARS automatically compiles the boulders detected in each image panel into a single table, containing the essential boulder data and coordinates.

MBARS Output and Calibration

MBARS outputs GIS-ready data tables, which describe each boulder object in the analyzed image. To streamline data analysis, only the most important parameters of each boulder object are included in this table. The image value records the partition of the original image in which the boulder is located. The flag value (as discussed in Section 2.4.2) is a unique ID to the boulder within its partition. The xloc and yloc values give the x and y locations of the boulder in meters Easting and Northing in the projection of the HiRISE image. The bouldwid, bouldheight, and shadlen values give the boulder diameter and boulder height in meters and the shadow length in pixels. The last four parameters are related to the results of the ellipse-fitting process: measured, fitgood, and fiterr record whether the boulder was measured by the algorithm, the quality of the fit of the ellipse to the shadow boundary, and the error on the ODR fit, respectively.

By design, MBARS does not strongly filter the results except where choosing not to filter causes computational concerns (e.g., large boulders crashing the fitting algorithm). As a result, many shadow-casting objects in the final results need to be removed for accurate statistical assessment of boulder populations. Objects with fitgood =0 are potentially problematic objects, and fitgood =1 suggests a confident boulder detection. The fitgood parameter can be flagged by several different conditions within MBARS, such as exceeding the maximum defined boulder width or height (30 m), exceeding maximum shadow area (3,000 pixels), or having a nonconvergent shadow fit. Boulders are exported into two files: FILENAME_Clean_boulderdata and FILENAME_All_boulderdata. The All file contains all boulder records, but the Clean data table only has the boulders for which fitgood =1. The problematic objects are preserved to help the user check for false negatives and adjust running parameters as needed, as nonconvergent fits can indicate poorly set running parameters. After importing the GIS-ready table to ArcGIS, the boulder results can be projected onto the image (using the “Display XY Data” function in ArcGIS) and analyzed further.

For each test area (Section 2.2), we calculate the CFA and best fit RA for both the manual data and the MBARS results, and the MBARS result that best matches the manually derived RA is chosen as final. Previous work on similar algorithms found that they operated best at boulder diameters between 1.5 and 2.25 m (Golombek et al., 2012). Following these findings, we calculate the best-fit RA for boulders in the test areas only between 1.5 and 2.5 m (Figure 5). Figure 5a demonstrates the comparison of several MBARS outputs and a set of manual results within test areas in HiRISE image TRA_000828_2495. In general, as the boundary parameter is decreased, shadows become smaller due to more restrictive filtering, and interpreted RA decreases. Figure 5b shows how this change in running parameters is expressed on individual boulders, showing that the most conservative and most generous interpretations of a boulder are often within 1–2 pixels of one another. In some cases, in order to match manual observations, boulders were further filtered according to their height-to-diameter ratio, removing objects with ratios above 0.75 (tall objects) or less than 0.25 (flat objects) in accordance with the standard deviation of h/D observed on the Moon (Demidov & Basilevsky, 2014).

Figure 5. An example of the calibration of Martian Boulder Automatic Recognition System (MBARS) with manual analysis (a) and how the MBARS-interpreted boulder changes with carrying running conditions (b). MBARS results are denoted “##_B,” where the ## represents the boundary parameter chosen for that run (Section 3.1) clipped to test area B (Figure 9).

Model Verification

Due to accompanying manual analysis, the accuracy of MBARS results is quantified for each image. However, further testing is warranted to ensure the accuracy of MBARS outside of the test areas and compare MBARS to other methods of boulder quantification. Three avenues are pursued here to demonstrate the accuracy of MBARS and compare it to other methods: (a) measuring objects of known size, (b) comparing results to previous manual analysis, and (c) comparing results to existing algorithms. In the plots below, error bars are smaller than the symbols on the plot.

Measuring Known Objects

Robotic landers represent one of the few objects on the Martian surface with physically quantified shape and size, making them the ideal target to ground truth MBARS. Here, these act as well-defined objects that can be measured with MBARS to test the accuracy of shadow-based measures as well as the consistency of MBARS measuring the same objects in different HiRISE images. Six images of the Viking 1 (VL1) and 2 (VL2) landers (Table 1) are used to test MBARS and allow comparison with the G-H method, which was also used to measure these landers (Golombek et al., 2008). Manual calibration areas were chosen in each image far from the vicinity of the landers to further demonstrate the generality of the solution generated in each image. In the case of the VL1 images, a suitably rocky area could not be identified for calibration, so a selection of widespread, individual boulders was used to calibrate MBARS instead using the same techniques.

Table 1 Example of Objects of Known Size Measured With MBARS

Lander	HiRISE image	Sun incidence angle, deg	Lander dimensionsa		G-Hb observations		MBARS dimensions
			D (m)	H (m)	D (m)	H (m)	D (m)	Hm (m)	Ha (m)
VL1	PSP_001521_2025	48	2.7	∼1.5	2.01	1.67	2.01	1.65	1.8
VL1	PSP_001719_2025	50	2.7	∼1.5	2.18	2.22	2.5	1.98	2.2
VL1	ESP_046170_2025	50	2.7	∼1.5	I	I	2.02	1.75	1.9
VL2	PSP_001976_2280	58	2.7	∼1.5	2.54	1.41	2.3	1.7	1.8
VL2	PSP_002055_2280	57	2.7	∼1.5	2.46	1.73	2.6	1.8	2
VL2	PSP_001501_2280	51	2.7	∼1.5	2.61	1.64	2.6	1.9	2.1

Note. HiRISE, high resolution imaging science experiment; MBARS, Martian Boulder Automatic Recognition System. Most images provide results that are within 1-pixel uncertainty of the known diameter. The measured heights (H_m) and actual heights (H_a, Figure 1) are shown with H_a∼20 cm taller than H_m in most cases. Images marked I (for Inapplicable) were either not measured in previous works or not available at the time of publication. MBARS should generally underestimate the lander width as the full width could only be measured if the landers were perfectly perpendicular to the Sun incidence direction.

^aViking lander heights are taken as from the top of the cameras (Huck et al., 1976), which we expect to cast a more consistent shadow than the full 2 m height from the top of the antenna.

^bFrom Golombek et al. (2008).

HiRISE images of VL1 and VL2 are shown in Figure 6 with the two bright spots showing the two wind covers on either side of these landers. In these images, the landers cast a sharp, dark shadow, and MBARS always succeeded in detecting and measuring the shadow. In general, the rock diameters given by MBARS for the Viking landers are smaller than the 2.7 m diameter of the landers (Table 1). Across the six images, MBARS measures the VL1 lander between 2.0 and 2.5 m and VL2 between 2.3 and 2.6 m. Published analyses with the G-H method yielded similar diameters as MBARS and similarly measured slightly higher dimensions for VL2 than VL1 (Table 1). While there may be some role played by the local topography or photometric properties of the two locales, the most likely cause of difference between the VL1 and VL2 measurements is lander orientation relative to the Sun. The full 2.7 m width of the lander is measured at its widest point, along the axis of the two radioisotope thermoelectric generator wind screens, but this widest dimension is not necessarily captured by its shadow. The orientation of VL1 and VL2 is well-constrained thanks to correlation of rock populations from orbital and ground imagery (Golombek et al., 2008). VL 1 is facing E-SE with its long axis pointing roughly NE-SW, and VL2 is facing NE with its long axis pointing NW-SE. In the VL1 images, the Sun is oriented closer to the long-axis, making the apparent diameter smaller, and in the VL2 images, the Sun is oriented more perpendicular to the long axis of the lander. In short, their orientations relative to the Sun reasonably explain the slightly higher estimations for VL2, though neither are consistently measured at their full width. For these observations, an extra step was taken to calculate the actual height (H_a) of the landers (see Figure 1, Equation 4). These corrections amount to an ∼20 cm increase in height between H_m and H_a for all observations, though the precise difference depends on the Sun incidence angle and boulder width. However, with this correction, all observations are within 20 cm (<15% error) of the known height of the lander, suggesting good agreement between our spheroidal model and the landers.

Figure 6. Martian Boulder Automatic Recognition System (MBARS) results for the Viking landers: (a) VL1, PSP_001719_2025 and (b) VL2, PSP_001501_2280, with the Sun orientation noted in the corner. The Viking landers are two of few objects of known size on the Martian surface and are therefore a useful test case for MBARS. Table 1 gives results for four other images of landers.

Comparison to Prior Manual Results

Manual measurements are an essential calibration step for MBARS. However, the manual analyses presented here were taken with an understanding of how MBARS detected boulders, creating some biases in how the boulders were measured. It is worthwhile to compare MBARS results to other manual boulder counts to determine if existing manual data sets are useful for calibration. We use the results of a pilot study (Sholes et al., 2017) that manually measured boulders in Acidalia Planitia using a different technique to determine the impact of manual boulder measurement technique on MBARS. Some discrepancies are expected between the two manual methods, since the user was directly outlining the bright boulder faces, rather than using shadows to measure the boulders, but this difference may not be significant enough to impact the MBARS calibration. We measured a separate area within the same images to determine whether discrepancies between MBARS and prior manual boulder measurements were related to a difference in measurement techniques or fundamental challenges with the image (Section 2.2). The comparison between automated and manual results is shown in Figure 7, using four areas within HiRISE image PSP_007718_2350.

Figure 7. Comparison of Martian Boulder Automatic Recognition System (MBARS) results (circles) to manual measurements (solid lines) in subsets of PSP_007718_2350 with a 50-m-wide portion of Test Area 3 shown on the right. Three images (10, 20, and 30) counted manually in prior work, and area A was manually counted as part of this work. Both sets of manual analyses had the greatest match with the same MBARS boundary parameter (Table S1 in Supporting Information S1). Inset of Area 3 is also included with MBARS-identified boulders (green outlines) over the image. As the inset shows, several bright objects lack apparent shadows in this image. Those may be mostly buried boulders that do not protrude sufficiently to cast shadows or simply boulders that deviate significantly from the expected shape. Either scenario likely caused the discrepancy between the methods as the manual methods would have counted these as boulders while they are functionally invisible to MBARS.

Both the prior manual measurements and our manual measurements produce boulder size-distributions that fit well to expected CFA curves (Table S1 in Supporting Information S1), suggesting that there are not systematic issues with the prior manual analyses. In all cases, MBARS underestimates the abundance of smaller boulders (<2 m) and overestimates the population of larger boulders, especially in Area 1. This overestimate of larger boulders appears to be a systematic issue (Section 4.4), but it is particularly pronounced in this image. Significantly, when using either the prior manual measurements, or our new manual measurements, the best-fit boundary parameter of MBARS is the same (Table S1 in Supporting Information S1). The change in manual measurement strategy did not significantly impact the choice of MBARS boundary parameter, suggesting that using existing manual measurements to calibrate MBARS is as effective as generating new manual measurements.

In these images, the tendency of MBARS to overpredict the CFA of larger boulders is particularly pronounced. In the three priorly measured areas, bright spots that lack shadows are abundant and visible in the image subset shown (Figure 7). In the prior manual analysis (Sholes et al., 2017), bright spots were counted by users as boulders, but the cast shadows are often either small or missing entirely. Boulders that lack shadows or have small shadows would interfere with the measuring methodology of MBARS. The solar incidence angle in the image is 48° from zenith, so 1 m boulders may be undetectable due to short shadows (Section 2.1). However, several ∼2 m wide bright spots are visible that lack significant shadows (Figure 7), and a boulder of that size and with an h/D ∼0.5 (Section 2.3) would be expected to cast a shadow at this incidence angle. In this image, there are then two factors challenging MBARS: A low incidence angle and an abundance of short, wide boulders. Despite these challenges, MBARS can recreate RA values within ∼10%–20% of the manually estimated values (Table S1 in Supporting Information S1) and calibration to manual analyses provides in situ quantification of these biases and uncertainties.

Comparison to Previous Algorithms

As a third comparison, we compare MBARS to both the G-H method and N-M methods of boulder detection (Golombek et al., 2008; Nagle-McNaughton et al., 2020). The three subsets of TRA_000828_2495 (A,B,C) isolated in previous work (Figure 17 in Golombek et al., 2008) were manually measured as an independent evaluation and used as the test areas for calibration purposes.

Comparison to the G-H Method

Similar to MBARS, the G-H method detects boulders based on their shadows and uses these shadows to determine boulder morphology. Figure 8 shows results for the test image from MBARS, the G-H method, the N-M method, and two sets of manual analyses within test areas A and B (Figure 9). These same results are shown for a portion of Test Area B in Figure 4. A third area (Area C) was previously analyzed (Golombek et al., 2008), though both our manual and automated analyses of area C found fewer than 5 boulders above 1.5 m in diameter, making the comparison of statistics in this area unreliable.

Figure 8. Comparison of Martian Boulder Automatic Recognition System (MBARS) results, Golombek-Huertas (G-H) method results (solid lines), Nagle-McNaughton (N-M) method conservative results (squares), and manual results (dashed lines) in test areas A (yellow colors) and B (blue colors). Best-fit RA lines are included for areas A (220%) and B (37%). G-H method results and definitions of areas A and B were taken from prior publications (Golombek et al., 2008) and shown in Figure 9. The MBARS results presented are the best fit to the manual results in areas A and B. The N-M method overestimates the boulder population in both areas even in the conservative (low) results, so the liberal (higher) results are excluded. MBARS and the G-H method are more comparable in both areas. Area B, which has a less-dense boulder population than A, shows closer agreement to MBARS results.

Figure 9. High resolution imaging science experiment image TRA_000828_2495 annotated with analysis Areas A and B as measured in Figure 8. This image is an example of typical terrain in the northern lowlands: low relief, moderate craters, and polygonal terrain covering most of the surfaces.

In area A (high boulder density), MBARS and the G-H method produce similar results. Both predict an RA of ∼220% suggesting an overabundance of boulders within the 1–2.5 m boulder range sampled here. However, both methods also overpredict the CFA compared to manual results in area A (Figure 8). In both cases, this overestimate is most obvious at the larger boulder diameters. In area B (medium boulder density), there is greater disparity between the MBARS and G-H results. MBARS results are much closer to the manual results (∼37% RA (Figure 8)), though the same tendency of overestimation at higher boulder diameters is present. In this area, the G-H method underpredicts the boulder abundance (∼20% RA) compared to manual measurements.

We also use a second image, PSP_001391_2465 to compare the performance of MBARS and the G-H method. In this image, the derived boulder densities and the cumulative number of rocks (CNR) per m² are compared between the two results (Figure 10). The maps demonstrate an agreement in boulder density variations, but MBARS detects a lower boulder density in all locations. The accompanying CNR plot (Figure 10) shows that G-H results underpredict boulder abundances at diameters above 1.5 m compared to manual and MBARS results despite having a higher boulder density. Because boulder density is totaled across all boulder sizes, the higher boulder density is likely caused by a higher abundance of low-diameter boulders predicted by the G-H method (Figure 10).

Figure 10. Comparison of observed rock density between Martian Boulder Automatic Recognition System (MBARS) results (top) and previous work (bottom, Figure 19, Golombek et al., 2008) in high resolution imaging science experiment image PSP_001391_2465 and a comparison of cumulative number of rocks/m2 (CNR, right) in representative areas. The maps show the same pattern, of high and low boulder density areas, but MBARS results are generally lower. The CNR plot shows that the G-H results consistently underpredict the boulder abundance compared to our manual analysis, while MBARS predicts the CNR more accurately. The unnamed crater rim shown in Figure 11 and discussed in Section 5 is marked.

Overall, these comparisons show that MBARS is either an improvement on, or performs comparably to, the G-H method in the tested images. At higher boulder densities, both algorithms tend to overpredict the CFA of boulders, but at more moderate densities, MBARS produces results that more closely match our manual observations. Both algorithms tend to overestimate boulder CFA at higher boulder dimensions, which may be due to a shared systematic error related to shadow-based boulder detection (Section 4.4). Both techniques yield compatible results, showing the same broad trends in RA across individual images. However, manual analyses in each image provide a clear reference point for each MBARS result and help quantify any uncertainties or errors associated with MBARS.

Comparison to the N-M Method

The N-M method of boulder size estimation uses various tools within ArcGIS to estimate the size and location of boulders based primarily on brightness contrasts (Nagle-McNaughton et al., 2020). We applied the N-M method in HiRISE image TRA_000828_2495 in order to compare among all three automated techniques. The N-M method was applied as described in the original publication, including the filtering steps to remove 1-pixel-sized boulders and boulders >10 m in any dimension (Nagle-McNaughton et al., 2020). One of the central steps of the N-M method is a 2 × 2 range filter, where areas with a large local range in DN are assumed to be boulders. Following the original methods, we define conservative (C) and liberal (L) thresholds (C = 190, L = 140), which are then used to define the high (Liberal) and low (Conservative) bounds on the boulder population. In the original work, the image was partitioned into different brightness ranges, and different C and L thresholds were defined for each. In the test image used here, there is no large-scale shadow-casting topography or major variation in surface albedo, so we felt that this step was not necessary and only used one brightness range.

Results of the N-M method are included in Figures 4 and 8 compared to MBARS and G-H results. In both areas A and B, the N-M method provides the highest estimate of boulder abundance with even the conservative (plotted as N-M Low) threshold yielding higher estimates of boulder abundance than the other methods. Most importantly, the N-M method estimates higher boulder abundances than the manual analysis in both areas, suggesting that this is an overestimate of the boulder population. Visual inspection of the N-M results (Figure 4) shows that boulder shadows were more clearly outlined in most cases than the boulders themselves. In these images, the boulder shadows contrast sharply with the surrounding soil, more so than the boulders themselves, which likely explain this result. In the original application, the N-M method was used near the landing site of the Perseverance Rover in Jezero Crater, specifically using HiRISE image PSP_002387_1985. In a visual comparison of these images, the boulder-soil brightness contrast is much stronger in Jezero than in the test image used here, where shadows are a stronger contrast to the bright surface. Of the three automated methods, the N-M method provides solutions that are least compatible with manual analyses in this image. However, photometric properties of the surface play a large role in all of these methods, and the soil-boulder as well as the soil-shadow contrast may dictate which method will be most accurate in a given study area.

Biases

One trend that is apparent in many of the CFA analyses is the tendency of MBARS to overestimate the CFA contribution of larger (>2 m) boulders compared to both manual analyses and theoretical distributions (Figure 8). One potential source of bias could be due to shortcomings in the techniques to correctly split and merge shadows that touch and overlap. This would be most prominent at larger boulder sizes where the range of irregular boulder morphology may be more discernible. A bias toward merging boulders that should be distinct would create an overestimate of large boulder abundances, an underestimate of small boulder abundances, and an overestimate of CFA at a low boulder diameter. Generally, MBARS does not overestimate the CFA at low boulder diameters, and any underestimate of smaller boulders could be ascribed instead to fundamental limits on an image resolution rather than false merges. While there are certainly cases where the merging and splitting steps fail to accurately divide or recombine the shadowed areas, it does not seem that this is a major source of bias.

Another possible source is related to the selection of shadow boundary, and how it may differently impact small and large boulders. Our shadow modeling (Section 3.1) determines a single value for DN, below which a pixel is considered part of a shadow. This value is based on the convolution of dark shadow pixels and brighter soil pixels according to the HiRISE PSF. However, as the boulders (and their shadows) increase in size, we expect less brightening of pixels near the shadow edge due to less convolution with bright background and greater convolution with the dark shadow interior. Therefore, as the shadow boundary value is increased, the model may become more accurate for smaller boulders but could overestimate the size of large boulders. Because the calibration method is used to match the manual analysis specifically in the 1–2.5 m range, this may lead to an overestimation of the size of larger boulders due to the higher DN chosen as a shadow boundary. However, using only the 1–2.5 m boulders to determine the RA and excluding the larger boulders help to mitigate this bias. Furthermore, calibrating MBARS results to manual analyses on an image-by-image basis informs whether or not this bias is present as well as how significant it may be in a given image or area within an image. In future work, users can choose how to calibrate their results depending on their scientific goals and what range of boulders is most significant for their analysis.

In addition to the above statistical examinations of MBARS, it is worth looking into the boulder-by-boulder accuracy of MBARS results. To further examine this accuracy, we use a 100-point grid-search approach to estimate false negatives, overestimation, underestimation, and accurate estimation rates in two of the test images (Figure S3 in Supporting Information S1). At each point in the grid (generated via the Fishnet tool in ArcGIS), we identify and manually measure the nearest boulder with diameter >1.5 m. Then, comparing these boulders to the MBARS results, we determine if MBARS failed to identify the boulder (false negative) mismeasured the boulder by >0.5 m (over/underestimation), or correctly identified and measured the boulder (accurate estimation). Plots of MBARS-measured boulder size versus manually determined boulder size show significant scatter around unity, suggesting mostly random, rather than the systematic error (Figure S3 in Supporting Information S1). Linear trends in the data are variable and do not support an overestimation of boulder sizes above 2 m in all images, but the relatively small sample size of this investigation (especially of boulders >2.5 m) prohibits confident extrapolation. Absolute detection rates (both accurate and mismeasured boulders) are 90% and 76% for each of the test images TRA_000828_2495 and PSP_001391_2465, respectively. This detection rate is comparable to modern Convolutional Neural Network-based techniques, such as the one used to select a landing site for the recent Tianwen-1 mission in Utopia Planitia (Wu et al., 2022). This suggests that MBARS is comparable in accuracy to more complex machine-learning-based approaches, the mechanics of which are often not deterministic and harder to understand and replicate a posteriori.

Application

Boulders are widespread throughout the images used above, but there are enhanced boulder populations apparent immediately surrounding the ∼0.5 km impact craters in each image (e.g., Figures 9 and 10). This suggests that these boulders likely originated from the impact events that formed these degraded depressions. Prior work on the Moon examined boulder distributions around relatively fresh, lunar impact craters to determine how both crater size and crater age influence the boulder distribution within the impact ejecta (Watkins et al., 2019). These observations were used to make inferences about both impact and boulder degradation processes on the Moon as well as test theoretical predictions of boulder sizes and ejection distances. This can be rapidly enabled by MBARS, outpacing the manual analysis substantially and increasing the sample size of examined craters. In order to demonstrate that MBARS can be used to rapidly and effectively replicate manual boulder surveys such as these, we isolate the boulders detected by MBARS near an unnamed crater (65.98 N, 103.95 W) in HiRISE image PSP_001391_2465 (Figure 11) and calculate statistics that are comparable to those determined in prior work (Watkins et al., 2019).

Figure 11. The unnamed crater in high resolution imaging science experiment image PSP_001391_2465 with no annotations (a) and rim, 1 and 2R circles (black lines) and Martian Boulder Automatic Recognition System-detected boulders as white dots (b). For context, the rim is marked in Figure 10. The rim was selected visually, aligning the rim just interior to the areas of increased boulder density to the South and North. The crater is estimated to have a rim diameter of ∼600 m and is close in size to both South Ray and Camelot crater (Watkins et al., 2019).

The unnamed crater has low relief, contains abundant thermal contraction polygons (Figure 11), and is likely a relatively young impact into an ice-rich substrate (Levy et al., 2018). The crater has a radius R ∼ 300 m, as determined visually by a ring of more abundant, large boulders and an apparent dearth of large boulders within the ring. Given the qualitative nature of this determination, the radius must be taken as an approximation, as no HiRISE DEMs are available in the area to provide stereo-derived topographic evidence of a rim. After defining the rim, we clipped boulders from the MBARS results out to 3R from the crater center. There are several other impact craters of comparable size and degradation state within the HiRISE image, so we cut off the survey at 3R to attribute the examined boulders more confidently to the impact event of interest. We compare results of the unnamed Martian crater to the results of Camelot (R = 303 m) and South Ray (R = 350 m) craters (Watkins et al., 2019) as these are closest in size to the unnamed crater. We examine the boulder size range distributions (BSRD) and boulder range frequency distributions (BRFD) as the size-frequency distributions are not expected to behave similarly (Golombek et al., 2012; Watkins et al., 2019).

Comparison of results among the BSRD and BRFD of the three craters is shown in Figure 12, including the BRFD of only boulders >1.5 m for the unnamed crater to compensate for the higher image resolution of HiRISE compared to the lunar images used in the original work. Comparison of BRSDs shows that much larger boulders (>7 m) are preserved at both the lunar craters, particularly at South Ray Crater, compared to the unnamed crater (D_max = 4.3 m). The unnamed crater replicates the trend of decreasing maximum boulder size at greater ejected distances seen at South Ray and Camelot Craters. However, the unnamed crater seems to show a steeper drop-off in large boulders, though this may be due to the lack of boulders ≳4.0 m. Comparing the BRFDs, the unnamed crater has a wide area of increased boulder abundance from ∼1.5 to 2R that has mostly diminished by 3R. This broad peak is also visible when only boulders >1.5 m are counted, showing that this is likely not an effect of increased sensitivity to small boulders. This distribution is unlike the sharp spike in boulder frequency at the Camelot rimor the more even distribution of boulders around South Ray crater.

Figure 12. Boulder range-frequency distribution (BRFD, top) and boulder size-range distribution (BSRD, bottom) for the unnamed Martian crater (Figure 11) and the two lunar craters: Camelot and South Ray (Watkins et al., 2019). Data for the unnamed crater are cut-off at 3R to minimize contamination from nearby craters. The data for Camelot and South Ray crater were originally expressed in crater radii from the crater rim (Watkins et al., 2019), and they are replotted here in radii from the crater center.

The comparison of BSRDs and BRFDs around these craters shows several intriguing similarities and disparities in the boulder populations. The unnamed Martian crater shows the expected radial decrease in boulder abundance and the maximum boulder size for impact ejecta. However, compared to craters of similar size on the Moon, there is a substantial lack of boulders with diameters above ∼4 m. In addition, the ejected boulders lack a prominent peak in frequency near the rim (as seen in Camelot) and diminish in abundance near 3 R (unlike South Ray). Further investigation would be required to determine if these differences are attributable to the differing erosional processes, regolith thicknesses, or target lithologies, as well as how these results can inform the relative rates of boulder degradation versus crater relaxation inferred for these boulder halo craters (Levy et al., 2018). However, such an investigation would be readily enabled by MBARS and could greatly expand the number of craters with fully characterized ejecta populations.

Conclusions

MBARS is an accessible, open-source tool for planetary scientists to locate and measure boulders in HiRISE images accurately and rapidly. We demonstrate that this tool, when applied in conjunction with supplemental manual analyses, provides a reliable technique to automatically assess large boulder populations with minimal manual effort. We show this in four key ways:

1. MBARS can accurately measure the height and width of the Viking landers, demonstrating the accuracy of these methods on objects of known size.

2. Manual boulder measurements taken for other purposes and surveys can be used to calibrate the automated results, reducing the need for new, redundant analyses.

3. MBARS is capable of measuring boulders (and landers) to the same or greater degree of accuracy with previously published methods of automated boulder measurements.

4. Comparison of MBARS analyses around a small crater compares well with prior analyses of similar-sized lunar craters.

Variations in boulder distributions can be distinctive of specific geologic processes, including impacts, mass wasting, and glacial processes. We hope that MBARS can be used as an essential tool going forward in many surface investigations of these processes. It is built to function on high-resolution images of the Martian surface but can be adapted to work on other bodies where sufficiently high-resolution data are available. This methodology enables many avenues of scientific investigation into impact, erosional, and glacial processes, which can now be pursued more rapidly and at a broader scope than is viable with the manual analysis alone.

Acknowledgments

Current support for this research comes from NASA's MDAP program, grant #80NSSC21K1093D. D. Hood has been funded throughout the project by the Louisiana Space Grant Graduate Student Research Assistantship awarded to D. Hood (NNX15AH82H) and NASA's Mars Data Analysis Program to S. Karunatillake and D. Hood (80NSSC18K1375-MDAP). The authors also thank Elizabeth McKinnie and J. Patrick Brothers for their boulder delineations used as part of this work. The authors thank NASA, JPL, and the University of Arizona for supporting access to HiRISE images.

Data Availability Statement

Data in support of this publication are located on the Texas Data Repository “Supporting Data for: ‘The Martian Boulder Automatic Recognition System, MBARS’” https://doi.org/10.18738/T8/EFPJWY (Hood et al., 2022b).The code to operate MBARS along with operating instructions and other information is available at https://github.com/dhood14/MBARS (Hood et al., 2022a). The version of MBARS used in this work is also available under the https://doi.org/10.5281/zenodo.6902208.