\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{The Martian Boulder Automatic Recognition System, MBARS}



\author[1]{Don R Hood}%
\author[2]{Steven F Sholes}%
\author[3]{Suniti Karunatillake}%
\author[4]{Caleb I. Fassett}%
\author[5]{Ryan C. Ewing}%
\author[6]{Joseph Levy}%
\affil[1]{Baylor University}%
\affil[2]{University of Washington}%
\affil[3]{Unknown}%
\affil[4]{Marshall Space Flight Center (NASA)}%
\affil[5]{Texas A\&M University}%
\affil[6]{Colgate University}%


\vspace{-1em}



  
  \date{June 15, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{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 (MRO).
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 and describe an automated tool to locate and
measure boulders on the martian surface: the Martian Boulder Automatic
Recognition System (MBARS). The open-source Python-based toolkit
autonomously measures boulder diameter and height in HiRISE images
enabling rapid and accurate assessments of boulder populations. We
compare our algorithm with existing boulder-counting methodologies,
manual analyses, and objects of known size to verify accuracy and
precision. Additionally, we test MBARS quantitatively characterizing
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.%
\end{abstract}%



\sloppy


\section*{The Martian Boulder Automatic Recognition System,
MBARS}

{\label{the-martian-boulder-automatic-recognition-system-mbars}}

\textbf{Don R. Hood\textsuperscript{1,2}, S.F.
Sholes\textsuperscript{3,4}, S. Karunatillake\textsuperscript{5}, C.I.
Fassett\textsuperscript{6}, R.C. Ewing\textsuperscript{2}, J.
Levy\textsuperscript{7}}

\textsuperscript{1}Department of Geosciences, Baylor University, Waco,
Texas 76706, USA

\textsuperscript{2}Department of Geology and Geophysics, Texas A\&M
University, College Station, Texas 77843, USA

\textsuperscript{3}Department of Earth and Space Sciences, University of
Washington, Seattle, Washington 98195, USA

\textsuperscript{4}Jet Propulsion Laboratory, California Institute of
Technology, Pasadena, California 91109, USA

\textsuperscript{5}Geology and Geophysics Department, Louisiana State
University, Baton Rouge, Louisiana 70803, USA

\textsuperscript{6}NASA Marshall Space Flight Center, Huntsville, AL
35805, USA

\textsuperscript{7}Department of Earth and Environmental Geosciences,
Colgate University, Hamilton, New York 13346, USA

Corresponding author: Don Hood (Don\_Hood@baylor.edu)

\subsection*{Key Points}

{\label{key-points}}

\begin{itemize}
\tightlist
\item
  The Martian Boulder Automatic Recognition System is a new tool to
  detect and measure boulders on the martian surface.
\item
  MBARS is comparably or more accurate than prior published algorithms
  that measure boulders.
\item
  MBARS can readily reproduce results originally acquired through manual
  measurement
\end{itemize}

\subsection*{Abstract}

{\label{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 (MRO).
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 and describe an automated tool to locate and
measure boulders on the martian surface: the Martian Boulder Automatic
Recognition System (MBARS). The open-source Python-based toolkit
autonomously measures boulder diameter and height in HiRISE images
enabling rapid and accurate assessments of boulder populations. We
compare our algorithm with existing boulder-counting methodologies,
manual analyses, and objects of known size to verify accuracy and
precision. Additionally, we test MBARS quantitatively characterizing
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.

\subsection*{Plain Language Summary}

{\label{plain-language-summary}}

Large boulders (\textgreater{}1 m diameter) are widely distributed on
the martian surface. They are easily observed from orbit, making them
visible with high resolution imaging. Mapping the location, number, and
size of boulders is helpful for understanding which geological processes
bring boulders to the surface, move them around, and fragment them into
smaller rocks and soil. Here, we present and describe the Martian
Boulder Automatic Recognition System (MBARS), a set of tools that
automatically locates and measures boulders in high-resolution images of
the martian surface. We compare results generated by MBARS with results
from other automated boulder-measuring tools, as well as with results
from manual boulder measurements to ensure accuracy. We also use MBARS
to map boulders around an impact crater on Mars and compare the boulder
distribution to a similar-sized crater on the Moon.

\subsection*{Keywords}

{\label{keywords}}

Remote Sensing (5464), Surface Materials and Properties (5470),
Instruments and Techniques (5494), Impact phenomena and cratering (5420)

\subsection*{1. Introduction}

{\label{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 (Golombek et al.,
2006; de Haas et al., 2013), and glacial processes (Levy et al., 2021).
Landing site assessments 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\emph{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, the 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 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 may be insufficient for some
investigations. Here we present an advancement from the existing boulder
measuring approaches with an open-source Python-based methodology to
automatically identify, locate, 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 known size, compare MBARS results to other
algorithms and manual analysis, 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).

\subsection*{2. Methods}

{\label{methods}}

\subsubsection*{2.1. The HiRISE Dataset}

{\label{the-hirise-dataset}}

The primary data for this work are images from the HiRISE camera, a
high-resolution camera that provides up to \textasciitilde{}25 cm/pixel
images of the martian surface. Typical HiRISE footprints are
\textasciitilde{}3 km x 6 km and images are globally distributed and not
contiguous due to their small size. The Point-Spread Function (PSF) of
HiRISE is \textasciitilde{}1.5 pixels (Kirk et al., 2008; McEwen et al.,
2007), so objects [?]1 meter across are resolvable at typical resolutions.
The PSF sets a lower bound on measurements made in HiRISE images, as
precision below \selectlanguage{ngerman}±1.5 pixels is not achievable without deconvolution or
other image-enhancing methods. 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 \emph{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 \selectlanguage{english}[?]1 meter, 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 (Fig. 1) we can explore if there are incidence
angles for which we do not expect to detect boulders within any given
size range. An object three pixels across is considered a reasonable
lower limit to reliably detect an object, so we use a minimum shadow
length of 3 pixels (\textasciitilde{}0.75 m) to determine the lower
incidence angle limit. The minimum detectable boulder height in meters
can be calculated as \(h_{\min}=0.75/tan(\theta)\), where \selectlanguage{greek}θ \selectlanguage{english}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_{\min}=1.4/tan(\theta)\). Due to the sun-synchronous orbit of MRO,
the incidence angle of most HiRISE images falls between
\textasciitilde{}40\selectlanguage{ngerman}° and \textasciitilde{}75°. In these images boulders
of diameters above 0.37 m and 1.67 m respectively, will cast detectable
shadows. Therefore, boulders that are sufficiently wide to be detected
(\selectlanguage{english}[?]1m) will also cast detectable shadows in most HiRISE images, and
boulders \textgreater{} 1.5 m wide are predicted to cast detectable
shadows in nearly all HiRISE images if they are present.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\selectlanguage{greek}\textbf{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 (H\textsubscript{m}), the actual height
(H\textsubscript{a}), shadow length (L\textsubscript{s}), boulder radius
(r) and the sun incidence angle (θ). The relationship between
H\textsubscript{m} and L\textsubscript{s} (Eq 1) is straightforward,
however determining H\textsubscript{a} from H\textsubscript{m} is more
complex.}\selectlanguage{english}

\subsubsection*{2.2. Manual Boulder
Analysis}

{\label{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 \emph{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
(\textasciitilde{}100 m x \textasciitilde{}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 10\textsuperscript{4} m\textsuperscript{2}, 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 10\textsuperscript{4}
m\textsuperscript{2}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 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 (\emph{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 F\textsubscript{k}(D) is the cumulative fractional area (CFA)
covered by rocks of diameter \textgreater{}\emph{D} in m\emph{, k} is
the rock abundance (expressed as a fraction), and the exponential
factor\emph{q(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., \textgreater{}100\% rock abundance,
especially in small study areas. Such a result in these analyses would
imply an excess of detectable (\textgreater{}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.

\subsubsection*{2.3. Detection by Shadow
Segmentation}

{\label{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 (\emph{H\textsubscript{m}} ) can be determined as:

\(\left(3\right)\text{\ \ \ \ \ \ }H_{m}=\ \frac{L_{s}}{\tan\theta}\)

Where \emph{L\textsubscript{s}} is the length of the shadow and
\selectlanguage{greek}\emph{θ}\selectlanguage{english}is the incidence angle of the sun, measured from zenith\emph{.
H\textsubscript{m}} will consistently be an underestimate of the actual
boulder height and the ratio of \emph{H\textsubscript{m}} to actual
height (\emph{H\textsubscript{a}} ) 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 \selectlanguage{greek}\emph{θ} \selectlanguage{english}is the incidence angle as in Eq. 3 and \emph{r} is the
radius of the boulder. If the boulder is assumed to follow lunar
boulders (i.e., H\textsubscript{a}/D = 0.5 or r = H\textsubscript{a}),
this underestimate will range from
H\textsubscript{m}/H\textsubscript{a}= 0.5 - 0.86 for \selectlanguage{greek}θ\selectlanguage{english}=30\selectlanguage{ngerman}°- 60°, with a
greater underestimate at lower \selectlanguage{greek}θ. \selectlanguage{english}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
\textless{}10\selectlanguage{ngerman}°, 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 (\textasciitilde{}6 mbar) and scattering cross
section of CO\textsubscript{2} for HiRISE RED filter wavelengths
(\textasciitilde{}1.214*10\textsuperscript{-38} m\textsuperscript{2} 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 \textasciitilde{}10 cm, which is
smaller than the 1.5 pixel (\textasciitilde{}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.

\subsection*{3. Algorithm Description}

{\label{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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2. A flowchart of the MBARS algorithm. The blue box shows
the required input data, all of which are available on the PDS. 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 sequence. The final product, the
GIS-ready list of boulder objects and attributes, compiles 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.}

\subsubsection*{3.1. Image Preparation and Shadow Boundary
Selection}

{\label{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 broken into panels, which we
carry out with the Split Raster tool in ArcMap This provides crucial map
orientation information to MBARS (as \emph{.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-1000 pixel square image panels
(\textasciitilde{}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 non-shadows, and modifies the brightness
boundary between those two classes to maximize the inter-class entropy.
MBARS takes a different approach, predicting a shadow boundary, i.e.,
the expected maximum brightness to be considered a part of a shadow, for
each image by forward modeling (Fig. 3). To predict a shadow boundary,
three key components are used: the darkness of shadows, the brightness
of non-shadow pixels, and the PSF of the HiRISE instrument. The interior
of shadows larger than \textasciitilde{}5 pixels are not perfectly dark
but are generally dark enough to register as Digital Number (DN,
corresponds to pixel brightness in the image) =1 in HiRISE images. The
brightness of non-shadow 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 on-board 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 \selectlanguage{greek}λ \selectlanguage{english}(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
non-shadow pixels, we construct a model image with a dark (DN=1) shadow
and convolve it with a background (non-shadowed area, Fig. 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-70 (i.e., 40\textsuperscript{th} and
70\textsuperscript{th} percentile) produce MBARS results consistent with
manual observations (Table S1).

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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{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 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.}

\subsubsection*{3.2. Boulder Segmentation}

{\label{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., Fig. 4d, S2). Local minima in the DN 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 (\textless{}4 pixel
area) to reliably identify, or too large to likely be boulders
(\textgreater{}30 m, i.e., very coarse blocks or other megagravel,
(Blair and McPherson, 1999)) are disregarded. With each shadow
identified and flagged, MBARS begins the process of determining boulder
morphometry based on the shadow dimensions.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 4. Portion of TRA\_000828\_2495 with a) no annotations,
b) Manual Results, c) N-M Method results (yellow = high, purple = low),
d) 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.}

\subsubsection*{3.3. Measuring Boulder
Dimensions}

{\label{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,
i.e., 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 (Fig. 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 Fig. 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
re-creating the full ellipse. The algorithm then, for each shadow, fits
an ellipse to the combined original and reflected shadow boundaries
(Figs. 1,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
(\textless{}1\selectlanguage{ngerman}° 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 (Figs.
1,4). Fits that return diameters or heights that are too large to likely
be a boulder (i.e., \textgreater{}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 mis-measured. At this stage, a list of shadow
objects, each containing information about the shadow, shadow ellipse,
and boulder, is recorded for each image panel.

After an image panel is complete, the measured boulders are re-visited
by MBARS, looking for overlapping shadows. First, all boulder objects
with shadows that touch (i.e., have adjacent 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\ldots{}n) up to \emph{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 towards 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\emph{.shad}
files, a custom filetype created using the C\emph{Pickle}module of
\emph{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.

\subsubsection*{3.4 MBARS Output and
Calibration}

{\label{mbars-output-and-calibration}}

MBARS outputs GIS-ready data tables which describe each boulder object
in the analyzed image. To save space, only the most important parameters
of each boulder object are included in this table. The \emph{image}value
records the partition of the original image in which the boulder is
located. The \emph{flag} value (as discussed in 2.4.2) is a unique ID to
the boulder within its partition. The \emph{xloc} and \emph{yloc}values
give the x- and y- location of the boulder in meters Easting and
Northing in the projection of the HiRISE image. The \emph{bouldwid}
,\emph{bouldheight} , and \emph{shadlen} values give the boulder
diameter and boulder height in meters and the shadow length in pixels.
The last four parameter are related to the results of the
ellipse-fitting process: \emph{measured,} \emph{fitgood,} and
\emph{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\emph{fitgood} =0 are potentially problematic objects,
and\emph{fitgood} =1 suggests a confident boulder detection.
The\emph{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 (3000 pixels), or
having a non-convergent shadow fit. Boulders are exported into two
files:\emph{FILENAME\_Clean\_boulderdata}
and\emph{FILENAME\_All\_boulderdata} . The \emph{All} file contains all
boulder records, but the \emph{Clean} data table only has the boulders
for which \emph{fitgood=1.} The problematic objects are preserved to
help the user check for false negatives and adjust running parameters as
needed, as non-convergent 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
which 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 m 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 (Fig 5). Fig. 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. Fig. 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).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Figure 5. An example of the calibration of 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 (Fig. 9).}

\subsection*{4. Model Verification}

{\label{model-verification}}

Due to accompanying manual analysis, the accuracy of MARS 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 MBARS and compare it to other methods:
(1) measuring objects of known size, (2) comparing results to previous
manual analysis, and (3) comparing results to existing algorithms. In
the plots below, error bars are smaller than the symbols on the plot.

\subsubsection*{4.1. Measuring Known
Objects}

{\label{measuring-known-objects}}

Robotic landers represent one of the few objects on the martian surface
with fully known 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. 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.

HiRISE images of VL1 and VL2 are shown in Fig. 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.7m diameter of the landers (Table 1). Across the six images, MBARS
measures the VL1 lander between 2.0 m and 2.5 m, and VL2 between 2.3 m
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 RTG 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\textsubscript{a}) of the landers (See Fig
1, Eq. 4). These correction amount to an \textasciitilde{}20 cm increase
in height between H\textsubscript{m} and H\textsubscript{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 20cm of the known height of the lander, suggesting good
agreement between our spheroidal model and the landers.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{Figure 6. 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 4 other images of
landers.}

\textbf{Table 1: Example of objects of known size measured with MBARS.
Most images provide results that are within 1-pixel uncertainty of the
known diameter. The measured heights (H\textsubscript{m}) and actual
heights (H\textsubscript{a}, Fig. 1) are shown, with
H\textsubscript{a}\textasciitilde{}20 cm taller than H\textsubscript{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.}\selectlanguage{english}
\begin{longtable}[]{@{}llllllllll@{}}
\toprule
Lander & HiRISE image & Sun Incidence Angle, deg & Lander
Dimensions\textsuperscript{a} & Lander Dimensions\textsuperscript{a} &
G-H\textsuperscript{b} observations & G-H\textsuperscript{b}
observations & MBARS Dimensions & MBARS Dimensions & MBARS
Dimensions\tabularnewline
& & & D (m) & H (m) & D (m) & H (m) & D (m) & H\textsubscript{m} (m) &
H\textsubscript{a} (m)\tabularnewline
VL1 & PSP\_001521\_2025 & 48 & 2.7 & \textasciitilde{}1.5 & 2.01 & 1.67
& 2.01 & 1.65 & 1.8\tabularnewline
VL1 & PSP\_001719\_2025 & 50 & 2.7 & \textasciitilde{}1.5 & 2.18 & 2.22
& 2.5 & 1.98 & 2.2\tabularnewline
VL1 & ESP\_046170\_2025 & 50 & 2.7 & \textasciitilde{}1.5 & I & I & 2.02
& 1.75 & 1.9\tabularnewline
VL2 & PSP\_001976\_2280 & 58 & 2.7 & \textasciitilde{}1.5 & 2.54 & 1.41
& 2.3 & 1.7 & 1.8\tabularnewline
VL2 & PSP\_002055\_2280 & 57 & 2.7 & \textasciitilde{}1.5 & 2.46 & 1.73
& 2.6 & 1.8 & 2\tabularnewline
VL2 & PSP\_001501\_2280 & 51 & 2.7 & \textasciitilde{}1.5 & 2.61 & 1.64
& 2.6 & 1.9 & 2.1\tabularnewline
\bottomrule
\end{longtable}

\textsuperscript{a}: Viking 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.

\textsuperscript{b}: From Golombek et al., 2008

\subsubsection*{4.2. Comparison to Prior Manual
Results}

{\label{comparison-to-prior-manual-results}}

Manual measurements are an essential calibration step when applying
MBARS. However, the manual analyses presented here were taken with an
understanding of how MBARS detected boulders, creating some bias in how
the boulders were measured. It is worthwhile to compare MBARS results to
other manual boulder counts to determine if existing manual datasets are
useful for calibration. 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. 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 technique or fundamental challenges with the
image (Section 2.2). The comparison between automated and manual results
is shown in Fig. 7, using four areas within HiRISE image
PSP\_007718\_2350.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

\textbf{Figure 7. Comparison of 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). 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
fully protrude from the subsurface, 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), suggesting that there are not systematic issues with the prior
manual analyses. In all cases, MBARS underestimates the abundance of
smaller boulders (\textless{}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). 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 over-predict 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 (Fig. 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\selectlanguage{ngerman}° from
zenith, so 1 m boulders may be undetectable due to short shadows
(Section 2.1). However, several \textasciitilde{}2m wide bright spots
are visible that lack significant shadows (Fig 7), and a boulder of that
size and with an h/D \textasciitilde{}0.5 (Sec. 2.3) would be expected
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 \textasciitilde{}10-20\% of the manually estimated values
(Table S1) and calibration to manual analyses provides \emph{in
situ}quantification of these biases and uncertainties.

\subsubsection*{4.3. Comparison to Previous
Algorithms}

{\label{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 (Fig. 17 in Golombek et al., 2008) were manually measured
as an independent evaluation and used as the test areas for calibration
purposes.

4.3.1. Comparison to G-H Method

{\label{comparison-to-g-h-method}}

Similar to MBARS, the G-H method detects boulders based on their shadows
and uses these shadows to determine boulder morphology. Fig. 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 (Fig. 9).
These same results are shown for a portion of Test Area B in Fig. 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

\textbf{Figure 8. Comparison of MBARS results, G-H method results (solid
lines), 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 area A and B were taken from prior
publications (Golombek et al., 2008) and shown in Fig. 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.}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image12/image12}

\end{center}
\end{figure}

\textbf{Figure 9. HiRISE image TRA\_000828\_2495 annotated with analysis
Areas A and B as measured in Fig. 8. This image is an example of typical
terrain in the northern lowlands: low relief, moderate craters, and
polygonal terrain covering most of the surface.}

In area A (high boulder density) MBARS and the G-H method produce
similar results. Both predict an RA of \textasciitilde{}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 (Fig. 8). In both cases this overestimate is most
obvious at the higher 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
(\textasciitilde{}37\% RA (Fig. 8)), though the same tendency of
overestimation at higher boulder diameters is present. In this area, the
G-H method underpredicts the boulder abundance (\textasciitilde{}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\textsuperscript{2} are compared between the two results (Fig. 10). The
maps demonstrate an agreement in boulder density variations, but MBARS
detects a lower boulder density in all locations. The accompanying CNR
plot (Fig. 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
(Fig. 10).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image13/image13}

\end{center}
\end{figure}

\textbf{Figure 10. Comparison of observed rock density between MBARS
results (\emph{top}) and previous work (\emph{bottom}, Fig. 19, Golombek
et al., 2008) in HiRISE image PSP\_001391\_2465 and a comparison of
cumulative number of rocks/m\textsuperscript{2} (CNR, \emph{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 Fig. 11 and discussed
in Section 5 is marked.}

Overall, these comparisons show that MBARS performs comparably to the
G-H method in the tested images. Using our manual results as a reference
point, 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.

4.3.2. Comparison to N-M Method

{\label{comparison-to-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 \textgreater{}10m in any dimension
(Nagle-McNaughton et al., 2020). One of the central steps of the N-M
method is a 2x2 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 Figs. 8 and 4 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 this is an overestimate of the boulder population.
Visual inspection of the N-M results (Fig. 4) show 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
explains 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 soil-boulder as well as
soil-shadow contrast may be a deciding factor in which analysis method
will be most accurate in a given study area.

\subsubsection*{4.4 Biases}

{\label{biases}}

One trend that is apparent in many of the CFA analyses is the tendency
of MBARS to overestimate the CFA contribution of larger (\textgreater{}2
m) boulders compared to both manual analyses and theoretical
distributions (Fig 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 discernable.
A bias towards 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 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 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 2.4.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 helps 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 has 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 (Fig. S3). At each
point in the grid (generated via the Fishnet tool in ArcGIS), we
identify and manually measure the nearest boulder with diameter
\textgreater{}1.5 m. Then, comparing these boulders to the MBARS
results, we determine if MBARS failed to identify the boulder (false
negative), mis-measured the boulder by \textgreater{}0.5 m
(over/under-estimation), or correctly identified and measured the
boulder (accurate estimation). Plots of MBARS-measured boulder size vs.
manually-determined boulder size show significant scatter around unity,
suggesting mostly random, rather than systematic error (Fig. S3). 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 \textgreater{} 2.5m)
prohibit confident extrapolation. Absolute detection rates (both
accurate and mis-measured 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 \emph{a posteriori} .

\subsection*{5. Application}

{\label{application}}

Boulders are widespread throughout the images used above, but there are
enhanced boulder populations apparent immediately surrounding the
\textasciitilde{}0.5 km impact craters in each image (e.g., Figs. 9,
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 ejected distances. This is the kind of
work that can be rapidly enabled by MBARS, outpacing manual analysis
substantially and increasing sample size of examined craters. To compare
with this work, we isolate the boulders detected by MBARS near an
unnamed crater (65.98 N, 103.95 W) in HiRISE image PSP\_001391\_2465
(Fig.11), and calculate statistics that are comparable to those
determined in prior work (Watkins et al., 2019).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image14/image14}

\end{center}
\end{figure}

\textbf{Figure 11. The unnamed crater in HiRISE image PSP\_001391\_2465
with rim, 1R and 2R circles (black lines) and MBARS-detected boulders as
white dots. For context, the rim is marked in Fig. 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 \textasciitilde{}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 (Fig. 11), and is likely a relatively young impact into an
ice-rich substrate (Levy et al., 2018). The crater has a radius
R\textasciitilde{}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 Digital Elevation Models
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 cutoff 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 are
shown in Fig. 12, including the BRFD of only boulders \textgreater{}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 show that much larger boulders (\textgreater{}7 m)
are preserved at both the lunar craters, particularly at South Ray
Crater, compared to the unnamed crater (D\textsubscript{max}=4.3 m). The
unnamed crater replicates the trend of decreasing maximum boulder size
at greater ejected distances seen 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 \textasciitilde{}1.5-2R that has mostly
diminished by 3 R. This broad peak is also visible when only boulders
\textgreater{}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 rim, or the
more even distribution of boulders around South Ray crater.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image15/image15}

\end{center}
\end{figure}

\textbf{Figure 12. Boulder Range-Frequency Distribution
(BRFD,\emph{top}) and Boulder Size-Range Distribution (BSRD,
\emph{bottom}) for the unnamed martian crater (Fig. 11), and the two
lunar craters: Camelot, and South Ray (Watkins et al., 2019). Data for
the unnamed crater are cutoff 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 re-plotted here in radii from the crater
center.}

The comparison of BSRDs and BRFDs around these craters show several
intriguing similarities and disparities in the boulder populations. The
unnamed martian crater shows the expected radial decrease in boulder
abundance, and 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 \textasciitilde{}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 vs. 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.

\subsection*{6. Conclusions}

{\label{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 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 compare well with
prior analyses of similarly 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
is available. This methodology enables many avenues of scientific
investigation into impact, erosional, and glacial processes, which can
now be pursued more rapidly and at broader scope than is viable with
manual analysis alone.

\subsection*{Acknowledgements}

{\label{acknowledgements}}

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.

\subsection*{Data and Software Availability
Statement}

{\label{data-and-software-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''' DOI: 10.18738/T8/EFPJWY .The code to
operate MBARS, along with operating instructions and other information
is available at https://github.com/dhood14/MBARS .

\subsection*{}

{\label{section}}

\subsection*{References}

{\label{references}}

Boggs, P. T., Byrd, R. H., Rogers, J. E., \& Schnabel, R. B. (1992).
User's Reference Guide for ODRPACK Version 2 . 01 Software for Weighted
Orthogonal Distance Regression. \emph{U.S. Department of Commerce} ,
(June), 99.

Demidov, N. E., \& Basilevsky, A. T. (2014). Height-to-diameter ratios
of moon rocks from analysis of Lunokhod-1 and -2 and Apollo 11--17
panoramas and LROC NAC images. \emph{Solar System Research} ,\emph{48}
(5), 324--329. https://doi.org/10.1134/s0038094614050013

Eppes, M. C., Willis, A., Molaro, J., Abernathy, S., \& Zhou, B. (2015).
Cracks in Martian boulders exhibit preferred orientations that point to
solar-induced thermal stress. \emph{Nature Communications} , \emph{6} .
https://doi.org/10.1038/ncomms7712

Golombek, M. P., \& Rapp, D. (1997). Size-frequency distributions of
rocks on Mars and Earth analog sites Implications for future landed
missions. \emph{Journal of Geophysical Research E: Planets} ,\emph{102}
(E2), 4117--4129. https://doi.org/10.1029/96JE03319

Golombek, M. P., Grant, J. a., Crumpler, L. S., Greeley, R., Arvidson,
R. E., Bell, J. F., et al. (2006). Erosion rates at the Mars Exploration
Rover landing sites and long-term climate change on Mars. \emph{Journal
of Geophysical Research E: Planets} , \emph{111} (12), 1--14.
https://doi.org/10.1029/2006JE002754

Golombek, M. P., Huertas, A., Marlow, J., McGrane, B., Klein, C.,
Martinez, M., et al. (2008). Size-frequency distributions of rocks on
the northern plains of Mars with special reference to Phoenix landing
surfaces. \emph{Journal of Geophysical Research} , \emph{113} , E00A09.
https://doi.org/10.1029/2007JE003065

Golombek, M. P., Huertas, A., Kipp, D., \& Calef, F. (2012). Detection
and Characterization of Rocks and Rock Size-Frequency Distributions at
the Final Four Mars Science Laboratory Landing Sites. \emph{The
International Journal of Mars Science and Exploration} , \emph{2} .
https://doi.org/10.1555/mars.2012.0001

Golombek, M. P., Kipp, D., Warner, N., Daubar, I. J., Fergason, R.,
Kirk, R. L., et al. (2016). Selection of the InSight Landing
Site.\emph{Space Sci Rev} , 5--95.
https://doi.org/10.1007/s11214-016-0321-9

de Haas, T., Hauber, E., \& Kleinhans, M. G. (2013). Local late
Amazonian boulder breakdown and denudation rate on
Mars.\emph{Geophysical Research Letters} , \emph{40} (14), 3527--3531.
https://doi.org/10.1002/grl.50726

Hecht, E. (2002). \emph{Optics} (4th ed.). San Francisco, CA: Addison
Wesley.

Huck, F. O., Burcher, E. E., Jobson, D. J., \& Wall, S. D. (1976).
Prediction of Viking Lander Camera Image Quality. \emph{NASA Technical
Note, TN D-8148} .

Kirk, R. L., Rosiek, M. R., Anderson, J. A., Archinal, B. A., Becker, K.
J., Cook, D. A., et al. (2008). Ultrahigh resolution topographic mapping
of Mars with MRO HiRISE stereo images : Meter-scale slopes of candidate
Phoenix landing sites. \emph{Journal of Geophysical Research}
,\emph{113} , 1--31. https://doi.org/10.1029/2007JE003000

Kneissl, T., Van Gasselt, S., \& Neukum, G. (2011).
Map-projection-independent crater size-frequency determination in GIS
environments - New software tool for ArcGIS. \emph{Planetary and Space
Science} , \emph{59} (11--12), 1243--1254.
https://doi.org/10.1016/j.pss.2010.03.015

Krishna, N., \& Senthil Kumar, P. (2016). Impact spallation processes on
the Moon: A case study from the size and shape analysis of ejecta
boulders and secondary craters of Censorinus crater. \emph{Icarus}
,\emph{264} , 274--299. https://doi.org/10.1016/j.icarus.2015.09.033

Levy, J. S., Fassett, C. I., Rader, L. X., King, I. R., Chaffey, P. M.,
Wagoner, C. M., et al. (2018). Distribution and Characteristics of
Boulder Halos at High Latitudes on Mars: Ground Ice and Surface
Processes Drive Surface Reworking. \emph{Journal of Geophysical
Research: Planets} , \emph{123} (2), 322--334.
https://doi.org/10.1002/2017JE005470

Levy, J. S., Fassett, C. I., Holt, J. W., Parsons, R., Cipolli, W.,
Goudge, T. A., et al. (2021). Surface boulder banding indicates Martian
debris-covered glaciers formed over multiple
glaciations.\emph{Proceedings of the National Academy of Sciences of the
United States of America} , \emph{118} (4), 1--7.
https://doi.org/10.1073/pnas.2015971118

Li, Y., \& Wu, B. (2018). Analysis of Rock Abundance on Lunar Surface
From Orbital and Descent Images Using Automatic Rock
Detection.\emph{Journal of Geophysical Research: Planets} .
https://doi.org/10.1029/2017JE005496

McEwen, A. S., Eliason, E. M., Bergstrom, J. W., Bridges, N. T.,
Delamere, W. A., Grant, J. A., et al. (2007). MRO's High Resolution
Imaging Science Experiment (HiRISE). \emph{Journal of Geophysical
Research} .

Nagle-McNaughton, T. P., Williams, J. M., Gallegos, Z. E., Wilkie, H. A.
A., Martinez, D. C., \& Scuderi, L. A. (2020). Geographic information
system based detection and quantification of boulders using HiRISE
imagery: a case study in Jezero Crater. \emph{Journal of Applied Remote
Sensing} , \emph{14} (1), 14522. https://doi.org/10.1117/1.JRS.14.014522

Sholes, S. F., Mushkin, A., \& Catling, D. C. (2017). Boulder-Size
Distributions as Indicators for Deposition Processes on Earth and Mars.
In \emph{GSA Annual Meeting} . Retrieved from
http://dx.doi.org/10.1130/abs/2017AM-304073

Watkins, R. N., Jolliff, B. L., Mistick, K., Fogerty, C., Lawrence, S.
J., Singer, K. N., \& Ghent, R. R. (2019). Boulder Distributions Around
Young, Small Lunar Impact Craters and Implications for Regolith
Production Rates and Landing Site Safety. \emph{Journal of Geophysical
Research: Planets} , \emph{124} (11), 2754--2771.
https://doi.org/10.1029/2019JE005963

Wu, B., Dong, J., Wang, Y., Rao, W., Sun, Z., Li, Z., et al. (2022).
Landing Site Selection and Characterisation of Tianwen-1 (Zhurong Rover)
on Mars. \emph{Journal of Geophysical Research: Planets} , \emph{1} ,
1--17. https://doi.org/10.1029/2021je007137

\subsection*{}

{\label{section-1}}

\selectlanguage{english}
\FloatBarrier
\end{document}

