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.