Figure 1. A geometric model of the map-view and side view of an observed boulder in a HiRISE 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 (Eq 1) is straightforward, however determining Ha from Hm is more complex.

2.2. Manual Boulder Analysis

MBARS produces several output datasets 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 m x ~100 m) which we use as test areas for calibration. These test areas are manually selected for each image based on several qualities. 1) A lack of major non-boulder objects. 2) Low relief and lacking obvious shadow-casting topographic features. 3) Representative areas of the different boulder densities in the image (i.e., moderate and high density). 4) Areas which 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.
We have an established methodology to measure boulders, but many publications exist that present boulders measured with varying techniques (e.g. Krishna & Senthil Kumar, 2016; Watkins et al., 2019). We sought to test how using these existing datasets, 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 reflected boulder surface, using shadows only to confirm boulder presence. Three areas (1, 2, and 3, each 104 m2, Fig 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 3-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 6-point ellipse tool. Adding to these measurements, we also defined and manually measured boulders in another 104 m2area, 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 when manual analyses already exist.
In each test area, we calculate the Rock Abundance (RA) based on the manual and automated observations. Rock Abundance 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):
\begin{equation} \left(1\right)\text{\ \ \ \ \ \ \ \ \ \ }F_{k}\left(D\right)=ke^{-q\left(k\right)D}\nonumber \\ \end{equation}
Where Fk(D) is the cumulative fractional area (CFA) covered by rocks of diameter >D in m, k is the rock abundance (expressed as a fraction), and the exponential factorq(k) is given by:
\begin{equation} \left(2\right)\text{\ \ \ \ \ \ q}\left(k\right)=1.79+\frac{0.152}{k}\nonumber \\ \end{equation}
The equation for q was empirically determined in previous work (Golombek & Rapp, 1997) and has shown good agreement with martian surface rock abundances at several scales (Golombek et al., 2012). It is worth noting that values of k can exceed 1, i.e., >100% rock abundance, especially in small study areas. Such a result in these analyses would imply an excess of detectable (>1.5m) boulders relative to expected distributions. Fits between the theoretical CFA and observed CFA are limited due 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.

2.3. 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 map view and an elliptical cross section when viewed from ground level (Fig. 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 re-orients 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 determined by the boulder height, boulder width, and incidence angle (Fig. 1). The shadow length is specifically related to the height on the boulder surface where the shadow is cast, labeled in the diagram as the “casting point”. Using the shadow length and the incidence angle of the sun, the measured height (Hm ) can be determined as:
\(\left(3\right)\text{\ \ \ \ \ \ }H_{m}=\ \frac{L_{s}}{\tan\theta}\)
Where Ls is the length of the shadow and θis the incidence angle of the sun, measured from zenith. Hm will consistently be an underestimate of the actual boulder height and the ratio of Hm to actual height (Ha ) can be determined as:
\begin{equation} \left(4\right)\text{\ \ \ \ }\frac{H_{m}}{H_{a}}=\ \frac{H_{a}\tan\theta}{\sqrt{H_{a}^{2}{\tan\left(\theta\right)}^{2}+r^{2}}}\nonumber \\ \end{equation}
Where θ is the incidence angle as in Eq. 3 and r is the radius of the boulder. If the boulder is assumed to follow lunar boulders (i.e., Ha/D = 0.5 or r = Ha), this underestimate will range from Hm/Ha= 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. It 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 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. Where sufficiently high-resolution Digital Elevation Models (DEMs) area 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 CO2 for HiRISE RED filter wavelengths (~1.214*10-38 m2 at 694 nm based on Mie Scattering (Hecht, 2002)) lead to expected scattering roughly ten orders of magnitude smaller than green light on Earth. The largest likely source of penumbrae on Mars is the size of the sun in the sky (0.33° (Golombek et al., 2008)), which creates penumbrae that are on the order of 3% of the shadow length. For boulders between 1.5 and 5m 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 the only source of shadow blurring.

3. 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 Fig. 2.