# Predicting Energy Consumption and Traversal Time of Ground Robots for
Outdoor Navigation on Multiple Types of Terrain

Matthias Eder1, Gerald Steinbauer-Wagner1 This work was partially supported by
the Austrian Research Promotion Agency (FFG) with the projects RoboNav and
FutureWoodTrans.1Matthias Eder and Gerald Steinbauer-Wagner are with the
Institute of Software Technology, Graz University of Technology, Graz,
Austria. {matthias.eder<EMAIL_ADDRESS>
###### Abstract

The outdoor navigation capabilities of ground robots have improved
significantly in recent years, opening up new potential applications in a
variety of settings. Cost-based representations of the environment are
frequently used in the path planning domain to obtain an optimized path based
on various objectives, such as traversal time or energy consumption. However,
obtaining such cost representations is still cumbersome, particularly in
outdoor settings with diverse terrain types and slope angles. In this paper,
we address this problem by using a data-driven approach to develop a cost
representation for various outdoor terrain types that supports two
optimization objectives, namely energy consumption and traversal time. We
train a supervised machine learning model whose inputs consists of extracted
environment data along a path and whose outputs are the predicted energy
consumption and traversal time. The model is based on a ResNet neural network
architecture and trained using field-recorded data. The error of the proposed
method on different types of terrain is within 11% of the ground truth data.
To show that it performs and generalizes better than currently existing
approaches on various types of terrain, a comparison to a baseline method is
made.

##  I Introduction

Navigation for ground robots in outdoor environments is a highly active
research domain where several disciplines meet. One of these disciplines is
path planning which aims to obtain an optimal path between two locations
within the robot’s environment based on given objectives [1]. Typical
objectives used in path planning are the safe traversal through the
environment avoiding the risk of collision [2], the minimization of the
robot’s energy consumption [3], its traversal time [4], or a combination of
multiple objectives. Finding a suitable representation for objectives is still
challenging, particularly for complex outdoor terrains. Simple objectives to
avoid non-traversable obstacles are effective on flat surfaces on similar
terrain, but they do not translate well to other types of terrain that can
have a significant impact on the objective. Furthermore, depending on whether
a robot goes up or down a slope in an outdoor environment, the outcome of an
objective for a robot in that location may be affected. This directional
influence on a given location complicates the representation of the planning
space used by the planning algorithms, which is typically a cost-based
representation of the environment [5], also known as a cost map. As a result
of these issues, determining cost values for path planning in outdoor
environments is still a challenge. Therefore, one of the key tasks in field
robotics is to define a planning space that accurately represents given
objectives in order to plan an optimal path through an environment with a
complex terrain [6]. Ideally, path planning considers several objectives at
the same time to find an optimal path. We take a step in this direction by
using the connection of two objectives, namely energy consumption and
traversal time, to create a cost representation that can be used to represent
both, energy and time optimal paths.

To generate cost maps for minimizing the energy consumption of a path, some
work in literature already exists using physics-based models [7, 8] or data-
driven approaches [9, 3]. These methods, however, only take into account a
limited or simplified understanding of the environment and the capabilities of
the robot. Similar to that, there are objectives to find time-optimal paths
that take into account the robot’s capabilities (e.g. acceleration limits
[10]) but do not consider how the environment affects the robot’s dynamics and
its limits. However, this is a critical element in finding time-optimal paths
in complex outdoor environments since specific environmental factors have an
impact on the robot’s capabilities and must be taken into account when
planning a path through the environment.

Figure 1: Overview of the proposed approach. A neural network is trained on
environment information around a path segment $P_{i}$, which consists of the
elevation $H_{i}$, terrain class $C_{i}$, and recorded orthophoto $O_{i}$ to
predict the average power $\hat{w}_{i}$ and velocity $\hat{v}_{i}$ of the
robot. The result is then used to estimate the traversal time $T$ and energy
consumption $E$.

To address these challenges, we present a data-driven approach that generates
a cost representation for the two aforementioned path-planning objectives,
energy consumption and traversal time, taking into account multiple types of
outdoor terrain. To do so, we use machine learning to train a target function
that can predict the energy consumption and traversal time along a path using
only environmental information along that path. To learn the target function,
a path is divided into segments of unit length which are used to extract
environment information around that path segment. The energy consumption and
traversal time of a robot along the path are predicted using a ResNet-based
network architecture that has been trained on a regression task. A field-
deployed ground robot is used to collect data in a realistic outdoor setting
on different types of terrain. Information about the environment is extracted
using photogrammetry and high-resolution aerial images. An overview of the
proposed method is shown in Figure 1.

The main contributions of this paper are 1) the presentation of a machine
learning-based method to predict the energy consumption and traversal time of
a ground robot in outdoor terrains using environment information; 2) the
incorporation of divers environmental knowledge to improve the prediction
performance on different types of outdoor terrain; and 3) a comparison of our
method to a baseline to show that our network outperforms the baseline method
and also eliminates the need to retrain the model for different types of
terrain.

The remainder of this paper is structured as follows. Section II discusses
related work in the field of cost estimation for outdoor traversal w.r.t.
energy prediction and traversal time estimation. Section III formulates the
problem, which is learned in Section IV. In Section V, a performance
evaluation of the proposed method and a comparison to a baseline approach is
presented. Section VI concludes this paper.

##  II Related Work

Typically, a cost map representation is used as planning space for path
planning algorithms to efficiently compute a path through the environment that
is optimal concerning a given objective. Such a representation can either
provide binary (traversable/non-traversable) [11], continuous [12], or
discrete cost estimates [13] following a given objective. There are multiple
objectives for optimizing the traverse in an environment such as risk of
collision, smoothness, traversal time, or energy consumption [14]. Another
approach is to model costs based on a robot’s motion stability by identifying
unstable motion patterns that should be avoided [15]. While some work focuses
on the estimation of a single objective (e.g. [8]), others define the
traversability as a combination of multiple objectives (e.g. [16]). In this
work, we focus on the simultaneous estimation of the objective function for
energy consumption and traversal time.

However, the development of an energy-based cost representation is a non-
trivial task for complex outdoor terrains as the energy consumption varies on
different types of terrain and thus physics-based methods lack in
applicability in non-uniform and rough terrains [17, 18, 7]. Quann et al. aim
to solve this task by acquiring friction coefficients from a designated area
and terrains first and then applying Gaussian Process Regression to estimate
the coefficients for new areas [8]. The method also proposes the usage of a
greyscale satellite image to gain more information about unseen environments.
While this method works well on terrain with available energy measurements, it
is also shown that it does not generalize well in new environments where no
data is available. Wei and Isler proposed a data-driven approach to predict
energy consumption based on the height information around the robot using deep
learning [3]. It is shown to work well on the same type of terrain but does
not generalize well on other types of terrain without retraining the network
on new data. Our proposed approach aims to overcome this problem by
incorporating terrain information in the trained model to accurately predict
the energy consumption of a robot on multiple types of terrain. Other factors
that affect the energy consumption of a ground robot are its dynamics, e.g.
acceleration and linear/angular velocity [7]. The authors of [19] present a
method to characterize the power of a robot for a skid-steered mobile robot
and optimizes the turning costs by deriving a mathematical model of friction
for skid-steered platforms. Other works focus on the optimization of a
trajectory for a ground robot which minimizes the energy consumed during
motion [20]. In this work, we focus on the estimation of the energy
consumption for a given path segment which can be used as an objective to plan
energy-optimal paths between two locations in the environment.

Planning time-optimal paths through an environment is one of the main
objectives in path planning and is explored in a variety of settings,
primarily focusing on the robot’s dynamics in a given environment [10, 21].
Therefore, research has been conducted to estimate the traversability of a
robotic system using different types of sensors and information, such as
proprioceptive measurements [22, 23], geometry-based representations [24], or
vision sensors [25]. While such traversability estimates can be used in
connection with the robot’s theoretical dynamics to compute a time-optimal
path, they do not take the influence of the environment on the robot’s
dynamics and its limits into account. To consider environmental factors, the
authors of [4] propose a time-optimal path planning algorithm for AUVs that
integrates the continuous dynamic flow of the ocean into the planning
procedure. Unlike for AUVs, terrain rarely changes for a ground robot.
However, certain environmental characteristics still have a limiting factor on
the robot’s capabilities and need to be considered when planning a path
through the environment. In this work, we aim to improve time-optimal path
planning by considering the environmental influence on the robot’s dynamics.
Therefore, we identify the robot’s translational velocity limits in different
environmental scenarios and use this information to improve the objective for
time-optimal planning.

##  III Problem Formulation

To present the used method for predicting energy consumption and traversal
time, we first formally describe the machine learning problem with the
corresponding concepts and notations. Therefore, assume we have a ground robot
that is driving along a path $\mathcal{P}$ in environment $\mathcal{M}$ and is
able to measure its energy consumption $E$ and traversal time $T$. This work
aims to predict $E$ and $T$ for any path in the environment using a target
function $F$:

$(T,E)=F(\mathcal{M},\mathcal{P}).$ (1)

For this, we define path $\mathcal{P}=\\{p_{0},...,p_{N}\\}$ as a list of $N$
points $p_{i}=(x_{i},y_{i})^{T}$, and environment $\mathcal{M}=\\{O,H,C\\}$ as
a set of three maps that store specific information about the environment. A
map is represented as a two-dimensional array holding environment information
for a dedicated point in the environment [26]. The map $H$ represents the
height map of an environment which provides information on the elevation of
the area. $O$ is a map holding information of a gray-scale orthophoto from the
environment and $C$ stores information about the type of terrain for each
global location.

When a robot is navigated through an environment, it can measure its energy
consumption $E$ using the battery voltage $V$ and its current $I$. Therefore,
when driving along a path $\mathcal{P}$ the energy consumption between two
points $p_{s},p_{g}$ can be calculated as the integral over time $t$, between
the two time instances $t_{s},t_{g}$ of $p_{s}$ and $p_{g}$:
$E=\int_{t=t_{s}}^{t=t_{g}}V(t)I(t)dt$. Since the energy of a robot is usually
measured in discrete time intervals, the energy consumption between $p_{s}$
and $p_{g}$ can be estimated over $K$ measurements as
$E\approx\sum_{i=1}^{K}V(t_{i})I(t_{i})\triangle t_{i}$, where $t_{0}=t_{s}$
is the time at point $p_{s}$ and $t_{K}=t_{g}$ the time at $p_{g}$. $\triangle
t_{i}=t_{i}-t_{i-1}$ is the time difference between the current and last
measurement. When the measurements of a robot are conducted in a constant
interval such that $\triangle t_{i}=\triangle t_{j},\forall i,j\in[1,K]$, the
equation can be reformulated as

$E\approx\sum_{i=1}^{K}V(t_{i})I(t_{i})\triangle
t_{i}=\frac{T}{K}\sum_{i=1}^{K}V(t_{i})I(t_{i})$ (2)

where $T=t_{K}-t_{0}$ corresponds to the difference between $t_{s}$ and
$t_{g}$ and thus to the traversal time as stated in Eq. 1.

##  IV Learning Method

The problem definition as formulated in Section III is difficult to solve
using machine learning, as the prediction of a path with arbitrary length
requires a large number of training samples and faces the challenge of how to
represent the paths as input. Furthermore, the estimation of energy
consumption $E$ is directly correlated with traversal time $T$, as shown in Eq
2. Due to their interdependence, it is challenging to independently learn both
variables in a single model. To resolve these issues, we break down the
formulation from Eq. 1 into a more manageable problem that can be solved using
machine learning, and then generalize it to the entire problem as stated
above.

To do so, we divide a path $\mathcal{P}$ into path segments of unit length
$d$, resulting in $l$ path segments $P_{i}=\\{p^{i}_{0},...,p^{i}_{n}\\},1\leq
i\leq l$ with $\sqrt{(p^{i}_{0}-p^{i}_{n})^{2}}=d$. $d$ is assumed to be small
enough so that it is a good approximation of the actual path length of $P_{i}$
such that $\sum_{j=1}^{n}\sqrt{(p^{i}_{j}-p^{i}_{j-1})^{2}}\approx d$. Next,
we define $M_{i}$ to be a patch of size $d\times d$, holding environment
information $\mathcal{M}$ which was extracted along the path segment $P_{i}$.
For a description of how $M_{i}$ is extracted, see Section IV-B. Since $M_{i}$
holds the environment information along $P_{i}$ and $d$ is an approximation of
the path length, it also indirectly holds information about the path $P_{i}$.

Following Eq. 2, we define $w=1/K\sum_{j=1}^{K}V(t_{j})I(t_{j})$ to be the
average power used on a path segment $P_{i}$ with $t_{0}$ being the time at
$p^{i}_{0}$ and $t_{K}$ being the time at $p^{i}_{n}$. Moreover, if $d$ is
small enough, $T$ can be estimated as the length of the path segment $d$
divided by the average linear velocity $v$ of the robot on the path segment:
$T\approx d/v$.

The assumptions made above allow us to resolve the challenges arising from Eq.
1 by learning a prediction function in a single machine-learning model on two
independent variables $w,v$ for a path of unit length $d$. Therefore, we train
a function $\hat{f}$ which takes an environment patch $M_{i}$ as input and
predicts the average power $\hat{w}_{i}$ as well as the velocity of the robot
$\hat{v}_{i}$ for traversing the patch:

$(\hat{w}_{i},\hat{v}_{i})=\hat{f}(M_{i}).$ (3)

The function is trained by minimizing the normalized root mean squared error
(NRMSE) between the prediction $(\hat{w}_{i},\hat{v}_{i})$ and the ground
truth $(w^{*}_{i},v^{*}_{i})$, which was recorded by the robot:

$\min\sum_{i=1}^{l}\sqrt{\alpha(\hat{w}_{i}-w_{i}^{*})^{2}+\beta(\hat{v}_{i}-v_{i}^{*})^{2}}$
(4)

with $\alpha=1/\max(w^{*})$ and $\beta=1/\max(v^{*})$ being the normalization
factors for both variables, and $l$ being the number of data samples used for
training.

Using the predictions of $\hat{f}$, the target function $F$ from Eq.1 can be
solved for any path $\mathcal{P}$ in $\mathcal{M}$ using the predictions
$\hat{f}$ for Eq. 2 over all $l$ path segments in $\mathcal{P}$

$(T,E)=(\sum_{i=1}^{l}d/\hat{v}_{i},\sum_{i=1}^{l}\hat{w}_{i}*d/\hat{v}_{i})$
(5)

Figure 2: Overview of the test environment. a) orthographic overview of the
experimental environment, visualizing the driven path with the robot for data
recording. The red area was used for training/testing and the rest for
validation, b) the height map of the environment, c) the segmented environment
with seven terrain classes of which the latter four are traversable by the
robot. White areas do not hold information about the environment.

###  IV-A Data Collection Process

This section describes the process of building a representation of the
environment and the data collection process using a ground robot to record the
ground truth data for prediction.

To build the maps of the environment, a commercial UAV equipped with a high-
resolution camera was flown over the recording scene to take overlapping
pictures. Photogrammetric processing was then conducted on the recorded images
as proposed in [27] using Agisoft Metashape111https://www.agisoft.com/
(Accessed on Feb 22${}^{\text{nd}}$ 2023) to generate an orthophoto $O$ and 3D
spatial data, which was used to generate the height map $H$. $H$ holds
altitude information ranging from $574.8m$ to $661.3m$ above sea level between
the lowest and highest recorded point. The orthophoto was manually segmented
to annotate the terrain types in $C$. Following the label definitions as used
in [28], a total of seven terrain classes were used for annotation, of which
four (grass, mud, paved/unpaved road) were set to be traversable by the robot.
A visualization of the three environment maps is shown in Figure 2. The
resulting maps covered an area of around $596m\times 305m$ and have a
resolution of $0.05m$.

Figure 3: The ground robot used for data collection on two different types of
terrain, namely mud (left) and grass (right).

To generate training data and ground truth measurements, the differential
drive robot Husky from Clearpath
Robotics222https://clearpathrobotics.com/husky-unmanned-ground-vehicle-robot/
(Accessed on Feb 22${}^{\text{nd}}$ 2023) was used. Figure 3 shows the robot
in the recording environment on two different types of terrain. It is powered
by a lead-acid battery with a nominal voltage of 24V. The system itself can
read the battery voltage and current consumed by the motors for propulsion.
Moreover, the robot was equipped with a RTK GNSS from GeoKonzept to record the
robot’s 3D pose with an accuracy of $5cm$. The global positioning data was
fused with the robot’s wheel odometry, an Xsens MTi-G-710 IMU, and the visual
odometry received from a Stereolabs ZED2 using an Unscented Kalman Filter
(UKF) to improve the velocity measurements during data recording. For data
recording, the robot was guided manually through the environment to cover all
the traversable terrains with varying slopes. Figure 2a shows the driven path
through the environment which was used for training and validation. During the
motion of the robot, the pose, battery voltage, current and fused velocity
were recorded by the system at a frequency of $20Hz$. To investigate the
velocity limits in the environment, the robot was permanently controlled
through the environment at its maximum speed of $1.0m/s$. Accelerations and
in-place rotations were discarded during the training data generation process.
The slope angles of the robot’s recorded path range from $-22.5^{\circ}$ to
$22.5^{\circ}$. In this work, only traversable surfaces are considered for the
proposed method. Therefore, we do not consider scenarios where the robot gets
stuck which would result in infinite energy consumption and a velocity equal
to zero. Following the patch extraction process from Section IV-B, a total of
$10.886$ patches were extracted of which $51\%$ were on grass, $9\%$ on mud,
and $20\%$ on paved and unpaved roads each.

###  IV-B Training Sample Generation

Figure 4: Patch generation process. A patch from each environment map is
extracted along the path segment of the robot and rotated to face East. This
is used as input for the neural network to predict energy consumption and
traversal time for the given patch.

To generate training samples for the training function $\hat{f}(M_{i})$, we
need to extract the environment patches $M_{i}$. The patch representation is
based on the approach from [3] and is used as it has information about the
path segment $P_{i}$ encoded and also contains information about the immediate
surrounding of the robot. For the training process, we define the length of
the path segment to be $d=1m$ and subsequently extract a patch of size
$1m\times 1m$ from all three maps $O,H,C$ held by $\mathcal{M}$. Therefore,
for any given path segment $P_{i}=\\{p^{i}_{0},...,p^{i}_{n}\\}$, we place
$p^{i}_{0}$ on the middle point ($d/2$) on one edge of patch $M_{i}$ and
select $p^{i}_{n}$ as the middle point on the opposite edge. This way, the
patch is rotated along the robot’s path in the environment by
$\theta^{\circ}$, where $0^{\circ}$ is heading towards the East. To use the
patch as a training sample on existing network architectures, it is rotated by
$-\theta^{\circ}$ such that $M_{i}$ can be represented as a matrix of size
$s\times s\times 3$, with $s=d/r$ depending on the path length $d$ and the
resolution $r$ of the environment information $\mathcal{M}$. As described in
Section IV-A, a resolution $r=0.05m$ is used, resulting in $s=40$ data points
per row and column for each layer and thus in a patch size of $M_{i}=(40\times
40\times 3)$. A constant depth of $3$ is chosen for $M_{i}$ to represent each
map ($O,H,C$) in a separate layer. The first layer in $M_{i}$ holds the grey-
scale values origin from the orthophoto $O$, scaled into $[0,1]$. The second
layer holds information on the terrain class $C$ which is scaled from $[1,c]$
to $[0,1]$, where $c$ is the number of terrain classes. As described in
Section IV-A, a total of $c=7$ terrain classes were used within this work, of
which four are traversable by the robot. The third layer holds the relative
height information $H$ of the path segment. Therefore, each absolute height
value in $M_{i}$ is subtracted by the minimum height value $\min(H_{i})$ of
the patch. The height layer is then also scaled into $[0,1]$ using a maximum
height difference of $1m$. The height difference was chosen as it represents a
$100\%$ grade on the path segment. An exemplary visualization of the sampling
process is shown in Figure 4.

###  IV-C Network Architecture and Training

Figure 5: The ResNet18-based network architecture used for training.

To learn the $\hat{f}(M_{i})$, the common ResNet-18 network architecture [29]
was used with an adapted head for image regression on two variables (see Fig.
5), as it has already proven to work for similar tasks [3]. For the training
process, environment patches $M_{i}$ of size $40\times 40\times 3$ are used as
input and two floating point values $w^{*}_{i},v^{*}_{i}$ are used to learn
the regression function. As loss function, the NRMSE as described in Equation
4 was used for the final layer. For training, a learning rate of $10^{-4}$ was
used. The network was trained using Keras. The collected data from Section
IV-A was split into two areas of which one was used for training/testing and
one for validation regarding an unseen environment. Figure 2a shows how the
collected data was separated. Furthermore, the data for training and testing
was split randomly to contain $80\%$ training samples and $20\%$ test samples.
This results in a total of $5.142$ training samples, $1.286$ test samples, and
$4.458$ samples for validation. The predicted outputs
$\hat{w}_{i},\hat{v}_{i}$ of the network were used to compute the estimated
traversal time $T$ and the energy consumption $E$ as shown in Equation 2.

##  V Performance Evaluation

In this section, we present the results of our proposed method in different
scenarios. First, we present the results of the training environment and show
that with the proposed approach it is possible to predict the traversal time
and energy consumption of a given path segment reliably on different types of
terrain. Second, we assess the performance of the network in a previously
unseen area and confirm the applicability of the trained network in unseen
environments. Third, to investigate the effect of the environment information
on the anticipated output, we conduct an ablation study on the various input
layers. Finally, we demonstrate that our approach outperforms the baseline and
eliminates the need to retrain the network for every type of terrain by
comparing it to the baseline, presented in [3].

To evaluate the performance of the trained networks, the Root Mean Squared
Error (RMSE) and Mean Absolute Percentage Error (MAPE) were chosen as metrics.
Both metrics were used for all subsequent evaluations, including the
comparison to the baseline. To compute the relative error for a single
prediction value, its prediction value ($\hat{x}$) and the ground truth
($x^{*}$) are used to compute the absolute percentage error (APE):
$e=\text{abs}(\hat{x}-x^{*})/x^{*}$.

###  V-A Results on the Known Test Environment

To evaluate the applicability of the proposed method on the task of predicting
the average power consumption and velocity of the robot driving through a
given input patch, we investigated the performance of the trained network on
test data that was recorded in the same environment where the data for
training the network was recorded. We summarized the results in Table I, which
describes the overall network performance over all given types of terrain
(’All’) and additionally investigates the performance on each type of terrain
to see how well the predictions work on different terrain. While the first two
rows show the performance of the trained values $w$ and $v$, the last two rows
show the performance of the computed traversal time $T=d/v$ and the energy
consumption $E$ as described in Equation 2.

TABLE I: Performance Results of the proposed approach on the test data,
extracted from the same environment which was used for training.

|  | Type of Terrain
---|---|---
Variable | Error | Grass | Mud |  | Unpaved
---
Road
| Paved
---
Road
All
| RMSE ($W$) | 14.35 | 13.92 | 7.81 | 16.09 | 14.04
$w$ | MAPE ($\%$) | 9.96 | 7.51 | 5.74 | 11.53 | 9.39
| RMSE ($m/s$) | 0.04 | 0.05 | 0.02 | 0.05 | 0.04
$v$ | MAPE ($\%$) | 3.32 | 5.27 | 1.78 | 3.64 | 3.41
| RMSE ($s$) | 0.06 | 0.10 | 0.03 | 0.07 | 0.06
$T$ | MAPE ($\%$) | 3.26 | 5.50 | 1.74 | 3.52 | 3.37
| RMSE ($J$) | 18.22 | 29.26 | 10.98 | 25.55 | 21.51
$E$ | MAPE ($\%$) | 11.06 | 10.64 | 6.15 | 12.90 | 10.77

The average error (MAPE) for predicting the power $w$ over all types of
terrain is $9.39\%$, and $3.41\%$ for velocity $v$. Using both variables to
compute the traversal time and energy consumption leads to an error of
$3.37\%$ and $10.77\%$. It can be seen, that the network is trained well on
all types of terrain, as the performance on all terrain types shows a good
prediction. The lowest MAPE for $w$ and $v$ could be observed on unpaved roads
($5.74\%$ and $1.78\%$), while the highest MAPE for $w$ was observed on paved
roads and for $v$ on mud. Overall, a good prediction performance could be
observed in the training environment. Figure 6 shows the prediction
performance for energy consumption $E$ on a number of selected data samples
from the test dataset on different types of terrain. It can be seen, that the
prediction works well for any type of terrain. This is also confirmed when
inspecting the relative error $e$ for the same data samples in Figure 7.

$20$$40$$60$$80$$0$$100$$200$$300$$400$$500$$600$Grass Unpaved Road Paved Road
MudPatch numberEnergy consumption (in $J$)Ground Truth $E^{*}$Prediction
$\hat{E}$

Figure 6: Exemplary comparison of the Energy prediction $\hat{E}$ and ground
truth $E^{*}$ of data samples from the test set on different types of terrain.

$20$$40$$60$$80$$0$$10$$20$$30$$40$$50$Grass Unpaved Road Paved Road MudPatch
number$\%$APE of $\hat{E}$

Figure 7: Absolute Percentage Error (APE) of the energy prediction $\hat{E}$
for exemplary data samples from the test set on different types of terrain.

###  V-B Generalization to unseen Environments

To show that our trained network generalizes well on unseen environments, we
investigate the performance of traversal time $T$ and energy consumption $E$
on the validation dataset, which was selected as shown in Figure 2a.
Therefore, a total of $4.458$ data points were used to evaluate the
performance of our approach in said environment. Table II shows the results
for the predicted values $w,v$, as well as for the computed objective values
$T,E$.

TABLE II: Performance Results of the proposed approach on the validation data,
which was recorded in an environment that was not used during training.

|  | Type of Terrain
---|---|---
Variable | Error | Grass | Mud |  | Unpaved
---
Road
| Paved
---
Road
All
| RMSE ($W$) | 18.85 | 29.20 | 15.07 | 16.15 | 17.73
$w$ | MAPE ($\%$) | 11.29 | 19.25 | 8.42 | 11.58 | 10.74
| RMSE ($m/s$) | 0.07 | 0.06 | 0.04 | 0.07 | 0.06
$v$ | MAPE ($\%$) | 5.60 | 5.22 | 2.90 | 5.70 | 5.14
| RMSE ($s$) | 0.09 | 0.09 | 0.04 | 0.09 | 0.08
$T$ | MAPE ($\%$) | 5.74 | 5.18 | 2.91 | 5.67 | 5.04
| RMSE ($J$) | 27.50 | 49.58 | 19.24 | 22.69 | 26.63
$E$ | MAPE ($\%$) | 13.04 | 22.46 | 10.12 | 13.69 | 12.80

The predictions of the average power $w$ and velocity $v$ on all terrains
deviate from the ground truth with a RMSE of $17.73W$ and $0.06m/s$ and a MAPE
of $10.74\%$ and $5.14\%$ respectively, which is close to the network’s
performance in the training environment. Moreover, the network seems to
generalize well in nearly all types of terrain. Only the error of $w$ (and
consequently $E$) on terrain ’Mud’ increases significantly on unseen data. One
reason for the error increase could be the relatively small proportion of
samples for terrain of type ’Mud’ in the validation set ($5.6\%$) which does
not fully reflect the properties of the terrain. Other than that, the network
shows good results in the prediction of traversal time and energy consumption.
This verifies that the proposed approach generalizes well to new environments
in surrounding areas since the errors are similar to the errors in the
training environment. Fig. 8 and Fig. 9 show the prediction performance for
traversal time $T$ on a number of selected data samples from the validation
dataset on different types of terrain.

$20$$40$$60$$80$$1$$1.2$$1.4$$1.6$$1.8$$2$Grass Unpaved Road Paved Road
MudPatch numberTraversal time (in $s$)Ground Truth $T^{*}$Prediction $\hat{T}$

Figure 8: Exemplary comparison of the estimated traversal time $\hat{T}$ and
ground truth $T^{*}$ of data samples from the validation set on different
types of terrain.

$20$$40$$60$$80$$0$$10$$20$$30$$40$$50$Grass Unpaved Road Paved Road MudPatch
number$\%$APE of $\hat{T}$

Figure 9: Absolute Percentage Error (APE) of the estimated traversal time
$\hat{T}$ for exemplary data samples from the validation set on different
terrain.

###  V-C Ablation of the Input Layers

To investigate the effect of the given environment information on the trained
model, an ablation study was conducted which systematically removes input
layers and observes the effect on the predicted outcome. The ablation of
individual layers has already been proven to be useful to investigate the
stability of the trained network on missing features and also to reveal causal
relationships of the features [30].

To evaluate the stability of the proposed model and to investigate the effect
of the environment features, we conducted an ablation of the trained network
by systematically replacing environment layers from the network input with
random noise. This way, the performance of the trained model could be
evaluated with different combinations of input layers. Table III gives an
overview of the performance results when ablating an environment feature. In
the context of this work, the ablation was conducted by not removing the
feature completely but by replacing the layer using random noise. The ablation
was conducted on the test set (’Test’), as well as on the validation set
(’Val’).

TABLE III: Performance Results of the proposed approach using ablated input
layers.

|  | $w$ | $v$ | $E$ | $T$
---|---|---|---|---|---
| Input
---
Patch
Set |  | RMSE
---
($W$)
| MAPE
---
($\%$)
| RMSE
---
($m/s$)
| MAPE
---
($\%$)
| RMSE
---
($J$)
| MAPE
---
($\%$)
| RMSE
---
($s$)
| MAPE
---
($\%$)
| Test | 14.04 | 9.39 | 0.04 | 3.41 | 21.51 | 10.77 | 0.06 | 3.37
$\\{O,H,C\\}$ | Val | 17.73 | 10.74 | 0.06 | 5.14 | 26.63 | 12.8 | 0.08 | 5.04
| Test | 19.33 | 10.94 | 0.05 | 3.84 | 30.35 | 12.29 | 0.08 | 3.71
$\\{H,C\\}$ | Val | 19.00 | 12.49 | 0.06 | 5.05 | 27.59 | 14.53 | 0.08 | 5.18
| Test | 14.78 | 9.86 | 0.04 | 3.62 | 22.17 | 10.92 | 0.06 | 3.63
$\\{O,H\\}$ | Val | 22.39 | 13.82 | 0.07 | 6.25 | 34.95 | 17.38 | 0.10 | 6.54
| Test | 50.63 | 41.02 | 0.07 | 5.69 | 67.59 | 47.84 | 0.09 | 5.64
$\\{O,C\\}$ | Val | 55.63 | 45.07 | 0.08 | 7.22 | 73.75 | 52.20 | 0.11 | 7.57
| Test | 54.52 | 44.63 | 0.09 | 6.49 | 57.34 | 57.17 | 0.11 | 6.41
$\\{O\\}$ | Val | 61.77 | 50.24 | 0.09 | 8.10 | 87.39 | 60.39 | 0.13 | 8.68
| Test | 66.37 | 44.10 | 0.07 | 5.73 | 91.93 | 47.16 | 0.11 | 5.32
$\\{C\\}$ | Val | 61.49 | 37.11 | 0.07 | 5.92 | 79.83 | 40.07 | 0.10 | 5.98
| Test | 19.33 | 11.69 | 0.05 | 3.71 | 28.67 | 13.31 | 0.07 | 3.60
$\\{H\\}$ | Val | 22.52 | 15.17 | 0.07 | 5.53 | 32.09 | 17.36 | 0.09 | 5.64

First, only a slight performance decrease can be observed when removing the
ortho patch on the input layer and using $\\{H,C\\}$ for prediction,
suggesting that the network is robust against the loss of ortho information.
Similarly to that, the evaluation of $\\{O,H\\}$ shows that the network
structure can cope with the loss of terrain information in the training
environment to some extent. However, an MAPE increase in energy prediction $E$
of nearly $5\%$ on the validation set suggests that information on terrain is
already an important feature for the generalization of the predictions on new
environments. When removing the height information ($\\{O,C\\}$), a
significant loss in performance can be observed, indicating that height
information is the crucial factor to predict energy consumption and traversal
time. This is further confirmed when looking at the performance results on a
single environment feature used as input. $\\{O\\}$ is not able to provide
stable predictions for the energy consumption and traversal time. $\\{H\\}$ as
a single input layer can provide relatively precise predictions, confirming
the work of [3]. However, to gain the full performance an different types of
terrain, the use of all three features as input is still required.

###  V-D Comparison to a baseline approach

To demonstrate the contribution of our proposed approach, we compare it to
baseline approaches which are used to predict the energy consumption $E$ and
the traversal time $T$. As no approach exists so far that provides both
objectives simultaneously, each objective is compared to its own baseline. We
show that our approach performs better than the baseline approaches and also
generalizes better to different types of terrain and new environments. For the
performance comparison, the same datasets as used in Section V-A and V-B are
used for baseline comparison in the training environment and for an unseen
environment.

TABLE IV: Performance comparison to the baseline approach. ’Wei & Isler [3]’
is only trained on terrain ’Grass’, while ’Wei & Isler [3] Retrained’ consists
of multiple networks which were separately trained on each type of terrain.

|  | Ours |  | Wei & Isler
---
[3]
|  Wei & Isler [3]
---
Retrained
| Prediction
---
Value
Set |  | RMSE
---
($J$)
| MAPE
---
($\%$)
| RMSE
---
($J$)
| MAPE
---
($\%$)
| RMSE
---
($J$)
| MAPE
---
($\%$)
| Test | 21.51 | 10.77 | 51.39 | 27.25 | 34.52 | 14.64
$E$ | Val | 26.63 | 12.80 | 56.16 | 28.37 | 32.92 | 16.96

For the comparison of the energy consumption, the method proposed by Wei and
Isler [3] is used as a baseline, as it has been proven to work well in outdoor
terrains and outperforms certain physics-based models. However, it has to be
noted that the baseline approach is not able to adapt to new types of terrain
within a single trained network. Therefore, our approach is compared to two
variants of the baseline: The first variant is the network trained solely on
data samples recorded on the terrain of type ’Grass’ (as done in [3]). For the
second variant, the approach from [3] is retrained separately on every type of
terrain to show the best possible performance of the baseline method on the
given terrain class. Table IV shows the overall performance comparison of our
approach to the baseline. As can be seen, when evaluated on multiple types of
terrain, our approach outperforms the baseline approach which was trained on
grassy terrain by more than $15\%$. This proves that our approach can adapt to
different types of terrain. Moreover, our approach outperforms the baseline by
$5\%$ in both the training and validation environments when compared to the
overall performance of the retrained baseline. This demonstrates that the
proposed approach can further improve energy prediction in outdoor
environments on multiple types of terrain and also generalizes better to
unseen environments.

To investigate the performance of the traversal time estimation, our approach
is compared to the expected traversal time which is based on the theoretical
system limitations of the robot. We compare the expected traversal time to the
outcome of our approach, which takes into account the robot’s immediate
surroundings while moving through a patch, to demonstrate that the expected
traversable time also depends on environmental factors such as slope and type
of terrain. Therefore, as a baseline, the robot’s expected velocity $v_{e}$
along a data sample is divided by the length $d$ of a path segment such that
$T_{e}=d/v_{e}$. As the robot is driven through the environment at the maximum
speed of $v_{e}=1m/s$ and the distance of a path segment in the recorded
dataset is $d=1$, the expected traversal time is a constant factor of one
second. Table V shows how the consideration of environmental limits can
improve the traversal time estimation, when compared to the expected traversal
time. It can be seen that the standard deviation (RMSE) of the actual time and
expected time is within $0.17-0.22s$ causing an error of up to $17.45\%$. By
estimating the velocity limit under consideration of the environment, an
overall error (MAPE) of $~{}5\%$ can be achieved. This shows, that the
consideration of the robot’s limits in different environmental settings is an
important factor when planning a time-optimal path through the environment.

TABLE V: Performance comparison to the expected time.

|  | Ours | Expected Time
---|---|---|---
| Prediction
---
Value
Set |  | RMSE
---
($s$)
| MAPE
---
($\%$)
| RMSE
---
($s$)
| MAPE
---
($\%$)
| Test | 0.06 | 3.37 | 0.22 | 17.45
$T$ | Val | 0.08 | 5.04 | 0.17 | 13.82

##  VI Conclusion

In this paper, we presented an approach to predict the energy consumption and
traversal time of a ground robot for path planning on different types of
terrain in outdoor environments. We showed that in our trained machine-
learning model, height information has the biggest impact on predicting the
energy consumption and traversal time of a given path. However, to improve the
overall prediction performance, more information about the environment and
type of terrain has to be considered. Therefore, we showed that the
incorporation of information on the type of terrain and visual information
from an orthophoto is a valuable factor to improve the prediction performance.
We showed that our method significantly reduces the error when predicting the
energy consumption when compared to a baseline approach and also gets rid of
the necessity to retrain a model for every type of terrain. We also showed
that the expected traversal time based on simple robot definitions from a fact
sheet cannot be used to precisely predict the traversal time in complex
terrain. Our approach, however, takes into account a robot’s immediate
surroundings to improve the estimation for traversal time. To further improve
the proposed approach, future work can incorporate rotational velocities in
the model, as rotations also play a crucial role in accurately predicting the
robot’s energy consumption and traversal time. Furthermore, the model can be
used as a baseline to adapt the trained information to new circumstances, such
as a new robot system or changing weather/seasonal conditions.

## References

  * [1] C. Alexopoulos and P. Griffin, “Path planning for a mobile robot,” _IEEE Transactions on Systems, Man, and Cybernetics_ , vol. 22, no. 2, pp. 318–322, 1992\.
  * [2] J. R. Sánchez-Ibáñez, P. J. Sanchez-Cuevas, and M. Olivares-Mendez, “Optimal and risk-aware path planning considering localization uncertainty for space exploration rovers,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 4092–4099.
  * [3] M. Wei and V. Isler, “Predicting energy consumption of ground robots on uneven terrains,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 1, pp. 594–601, 2022.
  * [4] T. Lolla, M. P. Ueckermann, K. Yiğit, P. J. Haley, and P. F. J. Lermusiaux, “Path planning in time dependent flow fields using level set methods,” in _2012 IEEE International Conference on Robotics and Automation_ , 2012, pp. 166–173.
  * [5] L. Jaillet, J. Cortes, and T. Simeon, “Transition-based rrt for path planning in continuous cost spaces,” in _2008 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2008, pp. 2145–2150.
  * [6] D. Fan, K. Otsu, Y. Kubo, A. Dixit, J. Burdick, and A.-a. Agha-mohammadi, “Step: Stochastic traversability evaluation and planning for risk-aware off-road navigation,” in _Robotics: Science and Systems XVII_. Robotics: Science and Systems Foundation, 2021\.
  * [7] K. Tiwari, X. Xiao, and N. Y. Chong, “Estimating achievable range of ground robots operating on single battery discharge for operational efficacy amelioration,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 3991–3998.
  * [8] M. Quann, L. Ojeda, W. Smith, D. Rizzo, M. Castanier, and K. Barton, “Off–road ground robot path energy cost prediction through probabilistic spatial mapping,” _Journal of Field Robotics_ , vol. 37, no. 3, pp. 421–439, 2020.
  * [9] M. Wei and V. Isler, “Air to ground collaboration for energy-efficient path planning for ground robots,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 1949–1954.
  * [10] M. Lepetič, G. Klančar, I. Škrjanc, D. Matko, and B. Potočnik, “Time optimal path planning considering acceleration limits,” _Robotics and Autonomous Systems_ , vol. 45, no. 3, pp. 199–210, 2003.
  * [11] W. Walch, M. Eder, K. Mautner-Lassnig, and G. Steinbauer-Wagner, “Offroad terrain classification for mobile robots,” _Austrian Robotics Workshop 2022_ , pp. 6–11, 2022.
  * [12] Z. Zhu, N. Li, R. Sun, D. Xu, and H. Zhao, “Off-road autonomous vehicles traversability analysis and trajectory planning based on deep inverse reinforcement learning,” in _2020 IEEE Intelligent Vehicles Symposium (IV)_ , 2020, pp. 971–977.
  * [13] J. Collier, G. Broten, and J. Giesbrecht, “Traversability analysis for unmanned ground vehicles: Interpreting the environment,” in _Defense Technical Information Center_ , 2006.
  * [14] C. Sevastopoulos and S. Konstantopoulos, “A survey of traversability estimation for mobile robots,” _IEEE Access_ , vol. 10, pp. 96 331–96 347, 2022.
  * [15] M. Prágr, P. Čížek, and J. Faigl, “Traversal cost modeling based on motion characterization for multi-legged walking robots,” in _2019 European Conference on Mobile Robots (ECMR)_ , 2019, pp. 1–6.
  * [16] Z. Jian, Z. Lu, X. Zhou, B. Lan, A. Xiao, X. Wang, and B. Liang, “Putn: A plane-fitting based uneven terrain navigation framework,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 7160–7166.
  * [17] S. Liu and D. Sun, “Minimizing energy consumption of wheeled mobile robots via optimal motion planning,” _IEEE/ASME Transactions on Mechatronics_ , vol. 19, no. 2, pp. 401–411, 2014.
  * [18] Z. Sun and J. Reif, “On finding energy-minimizing paths on terrains,” _IEEE Transactions on Robotics_ , vol. 21, no. 1, pp. 102–114, 2005.
  * [19] S. Dogru and L. Marques, “Power characterization of a skid-steered mobile field robot with an application to headland turn optimization,” _Journal of Intelligent & Robotic Systems_, vol. 93, no. 3-4, pp. 601–615, 2019.
  * [20] P. Tokekar, N. Karnad, and V. Isler, “Energy-optimal trajectory planning for car-like robots,” _Autonomous Robots_ , vol. 37, no. 3, pp. 279–300, 2014\.
  * [21] A. Wolek, E. M. Cliff, and C. A. Woolsey, “Time-optimal path planning for a kinematic car with variable speed,” _Journal of Guidance, Control, and Dynamics_ , vol. 39, no. 10, pp. 2374–2390, 2016.
  * [22] G. G. Waibel, T. Löw, M. Nass, D. Howard, T. Bandyopadhyay, and P. V. K. Borges, “How rough is the path? Terrain traversability estimation for local and global path planning,” _IEEE Transactions on Intelligent Transportation Systems_ , pp. 1–12, 2022.
  * [23] A. Jacoff, E. Messina, and J. Evans, “Performance evaluation of autonomous mobile robots,” _Industrial Robot: An International Journal_ , vol. 29, no. 3, pp. 259–267, 2002.
  * [24] Z. Jian, Z. Lu, X. Zhou, B. Lan, A. Xiao, X. Wang, and B. Liang, “Putn: A plane-fitting based uneven terrain navigation framework,” 2022. [Online]. Available: https://arxiv.org/abs/2203.04541
  * [25] A. Broggi, E. Cardarelli, S. Cattani, and M. Sabbatelli, “Terrain mapping for off-road autonomous ground vehicles using rational b-spline surfaces and stereo vision,” in _Intelligent Vehicles Symposium (IV), 2013 IEEE_. IEEE, 2013, pp. 648–653.
  * [26] A. Kleiner and C. Dornhege, “Real-time localization and elevation mapping within urban search and rescue scenarios,” _Journal of Field Robotics_ , vol. 24, no. 8-9, pp. 723–745, 2007. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/rob.20208
  * [27] C. Mayer, L. Pereira, and T. Kersten, “A comprehensive workflow to process uav images for the efficient production of accurate geo-information,” in _CNCG2018 - IX National Conference on Cartography and Geodesy_ , 2018, pp. 1–8.
  * [28] A. Valada, G. Oliveira, T.Brox, and W. Burgard, “Deep multispectral semantic scene understanding of forested environments using multimodal fusion,” in _International Symposium on Experimental Robotics (ISER 2016)_ , 2016. [Online]. Available: http://lmb.informatik.uni-freiburg.de/Publications/2016/OB16c
  * [29] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 770–778.
  * [30] R. Meyes, M. Lu, C. W. de Puiseau, and T. Meisen, “Ablation studies in artificial neural networks,” _CoRR_ , vol. abs/1901.08644, 2019. [Online]. Available: http://arxiv.org/abs/1901.08644


# Deep Prototypical-Parts Ease Morphological Kidney Stone Identification and
are Competitively Robust to Photometric Perturbations

Daniel Flores-Araiza1 Francisco Lopez-Tiro1 Jonathan El-Beze2 Jacques Hubert2
Miguel Gonzalez-Mendoza1 Gilberto Ochoa-Ruiz1 Christian Daul3
1Tecnologico de Monterrey, School of Engineering and Sciences, Mexico
2CHU Nancy, Service d’urologie de Brabois, Nancy, France
3CRAN UMR 7039, Université de Lorraine and CNRS, Nancy, France

###### Abstract

Identifying the type of kidney stones can allow urologists to determine their
cause of formation, improving the prescription of appropriate treatments to
diminish future relapses. Currently, the associated ex-vivo diagnosis (known
as Morpho-constitutional Analysis, MCA) is time-consuming, expensive and
requires a great deal of experience, as it requires a visual analysis
component that is highly operator dependant. Recently, machine learning
methods have been developed for in-vivo endoscopic stone recognition. Deep
Learning (DL) based methods outperform non-DL methods in terms of accuracy but
lack explainability. Despite this trade-off, when it comes to making high-
stakes decisions, it’s important to prioritize understandable Computer-Aided
Diagnosis (CADx) that suggests a course of action based on reasonable
evidence, rather than a model prescribing a course of action. In this
proposal, we learn Prototypical Parts (PPs) per kidney stone subtype, which
are used by the DL model to generate an output classification. Using PPs in
the classification task enables case-based reasoning explanations for such
output, thus making the model interpretable. In addition, we modify global
visual characteristics to describe their relevance to the PPs and the
sensitivity of our model’s performance. With this, we provide explanations
with additional information at the sample, class and model levels in contrast
to previous works. Although our implementation’s average accuracy is lower
than state-of-the-art (SOTA) non-interpretable DL models by 1.5%, our models
perform 2.8% better on perturbed images with a lower standard deviation,
without adversarial training. Thus, Learning PPs has the potential to create
more robust DL models. Code at:
https://github.com/DanielF29/Prototipical_Parts

Figure 1:  Traditional Morpho-constitutional Analysis (MCA) involves an
endoscopic procedure to locate and extract a kidney stone (Ex-vivo), and
subject it to visual and Fourier Transform Infrared Spectroscopy (FTIR)
analysis by a specialist for classification, in the MCA section of the figure.
Recent works have proposed using Deep Learning (DL) to classify kidney stones
during an endoscopic procedure (In-vivo) automatically, represented in the XAI
section of the figure. We argue and show, eXplainable AI (XAI) is needed for
the adoption of DL as a Computer Assisted Diagnosis (CADx) tool. In the orange
box is the extension to the current DL systems this proposal enables.

(a)

(b)

Figure 2:  (a) The model classifies correctly a kidney stone patch, sub-type
”WD”, by measuring the similarity of learned Prototypical-Parts (PPs) to parts
in the input image. As seen in (b), first, the most relevant area on the input
image is indicated by a yellow bounding box and a heatmap representing levels
of similarity on the input image for each Prototypical-Part (PP). Second, the
image with the representative part from the training set that looks like the
input indicated part is also shown on a yellow box, with its corresponding
heatmap. Third, the level of saliency of visual characteristics and the
descriptors are obtained. In this example descriptors for the most relevant PP
are shown, indicating hue as the most relevant feature for the activation of
this PP. Then, PPs identified work both as our explanations and for the
classification, rendering this architecture interpretable, easing the
urologist task of kidney stone classification, as they are provided with the
reasons used by the model for the suggested classification.

##  1 Introduction

Urolithiasis is a disease in which kidney stones form somewhere in the urinary
tract [28]. In developed countries, this disease has become a public health
problem and presents a high incidence of kidney stones episodes, affecting up
to 10% of their population [18]. The formation of kidney stones is strongly
related to factors such as diet [13, 25]. However, there are other factors
that can lead to the production of a kidney stone such as age, chronic
diseases, hereditary-family history, among others [4, 6]. The timely
identification of a kidney stone aids clinical specialists in determining the
causes of its formation and to prescribe a personalized treatment, which can
prevent relapses [13].

The morpho-constitutional analysis (MCA) is the standard technique for
determining the different types of whole kidney stones removed during surgery
(ex-vivo stones) [5, 6]. Using this method it is possible to identify 21
different types of kidney stones, with both pure and mixed compositions [4].
Although the MCA efficiently establishes the type of ex-vivo kidney stones, it
is very difficult to provide a reliable diagnosis during an endoscopic
intervention (in-vivo), and the results of the MCA on extracted stones (ex-
vivo) may take several days. Also, it can be difficult to carry out the
inspection during surgery, as the extraction process can take up to one hour
[4]. Furthermore, it is very difficult to train specialists on this technique,
even considering a high incidence of kidney stone episodes.

Given the importance and difficulty of performing the visual inspection in a
repeatable and highly reliable manner, recently, several Artificial
Intelligence (AI) methods have been proposed for automating the kidney stone
classification process. Using AI to automate the classification of kidney
stones entails either replacing MCA or aiding the surgeon during the
extraction process.

Several DL-based techniques in particular have demonstrated that is possible
to replicate the results obtained by well-trained specialists [15, 10, 2, 20,
21]. However, DL models lack interpretation of the extracted features, making
these models of limited help in clinical settings. To facilitate clinical
acceptance, a DL model will need to provide the specialist with its reasoning
process, this implies providing evidence on the what, how and why behind any
suggested classification [7, 26]. Therefore, to be a useful aid for CADx
purposes, any evidence should reflect the same key distinctions used for MCA
by a specialist.

In order to pave the way for AI-based interpretable MCA using deep learning
techniques, we leverage recent strides in explainability that seek to base
image classification on case-based reasoning [3, 22, 1]. In this work, both
visual explanations and quantitative information about visual characteristics
deemed important by the network are provided [19]. It must be emphasized that
the explanations provided by our approach follow the reasoning processes of
urologists during MCA by showing the detected morphological relevant features
of each image as learned PPs.

Overall, this work presents an interpretable DL classifier for MCA. This
proposal facilitates human-machine collaboration for morphological analysis
whose explicit reasoning can be easily inspected, understood and verified by
urologists and medical specialists. In contrast to previous works, we evaluate
our interpretable models and their DL counterparts under perturbations of
global visual characteristics to describe their relevance to the PPs and the
sensitivity of our model.

The of the paper is organized as follows. In the following Section 2 the
proposed approach reasoning and justification are presented. In Section 4 the
experimental setup is described for our experiments. In Section 5 we share our
results and present a discussion on their implications, sharing one
explanation case comparing a correct and incorrect classification. Finally, in
Section 6 the conclusion and the insights obtained for further research in the
field are presented.

##  2 Motivation

Figure 3:  Visual comparison of different types of model interpretability on
an example of a kidney stone image as an example. (a) Traditional non-
interpretable deep learning output classification of a kidney stone [15], (b)
object-level attention map (e.g., class activation map [24], (c) part
attention (provided by attention-based interpretable models [30]); and (d)
part attention with similar prototypical parts [3], with additional (e)
quantitative evaluations of visual features, the descriptors [12, 19].

MCA consists of two complementary analyses of the kidney stone ex-vivo in the
laboratory. First, a specialist performs a visual inspection of the
characteristics (such as color, texture and shapes) of the surface (external)
and section (cross-sectional) views. Second, a biochemical study is performed
to determine the main components in the stone through an Fourier Transform
Infrared Spectroscopy (FTIR) study [9]. Both studies complement each other to
provide a detailed report to determine the type of kidney stone. This
traditional MCA process is represented in the MCA area in Fig. 1.

Recently proposed AI systems for kidney stone classification are presented as
automated CADx tools. This approach could enable faster and more reliable
diagnosis based solely on kidney stone images, reducing the need for invasive
procedures, on in-vivo [10] or ex-vivo [15, 2, 20] conditions, as represented
in MCA-Based DL, the grey area at Fig. 1. Also, for in-vivo scenarios, the
operations would be quicker to perform and less traumatic, due to the
classification of the kidney stones right before performing dusting, a process
intended to fragment and destroy the kidney stones inside the urinary tract.
Current DL classifiers of kidney stones are black boxes, as those have the
limitation of being unable to describe the inner workings that led to a given
prediction beyond the class label. Therefore, these types of models cannot
provide useful information to the specialist to understand how features on the
input image were used to output the classification proposed by the DL model.

The eXplainable AI (XAI) field seeks to remedy these limitations. XAI methods
are used to provide an explanation of how DL models arrived at a particular
output, typically after the classification was assigned. If DL models achieve
high performance and can provide explanations for the classification based on
relevant features for the specialists, then it would be easier to maintain
accountability for recommended treatments based on those analyses. In the area
indicated by the XAI bracket at Fig. 1, a holistic explanation is represented,
by providing the description of what and why visual features are relevant from
the input image, as quantitative evaluations of visual features, and a
saliency map highlighting where the described features are located.
Explanations that could be considered complete do not yet exist for kidney
stone identification. Nonetheless, improving explanations by indicating
precise visual characteristics would facilitate the adoption of DL CADx tools,
resulting in faster, more accountable, and repeatable diagnosis
recommendations for the treatment of urolithiasis.

Taking into account the previously mentioned desired characteristics of
explanations, in this work, we seek to improve the current methods by adding
the required explainability to facilitate the morphological analysis of kidney
stones, our proposal of how to achieve this is depicted in Fig. 2. For
context, Fig. 3 presents a visual comparison of different explanations. Then,
to improve DL-generated explanations for the classification of kidney stones
we began generating saliency maps of learned representative parts, the learned
PPs, for each class [12], after, we evaluate which global visual
characteristics PPs are relaying for such classification [19].

##  3 Proposed Approach

Having a single heatmap as a visual explanation, as most current XAI
visualization methods do, is an oversimplification of all the characteristics
expected and used to classify a kidney stone into its corresponding type. To
alleviate the unmet needs of most of the current DL models and XAI
visualization methods, this proposal leverages a case-based reasoning process.
Our proposal pipeline does this by extracting semantic features, representing
parts, from an input image with a Convolutional Neural Network (CNN). Then,
the extracted features are measured against the learned prototypical parts, in
the form of similarity scores, as indicated by equation 1, and represented in
Fig. 2(b). Finally, the output prediction of the classification task is a
weighted combination of the similarity scores. Using these same learned PPs
for the output of the model and the explanations, the method guarantees
faithfulness of the explanations to the inner workings of the model. Current
case-based reasoning methods, instead of using high-dimensional features from
the convolutional layers of a CNN, use the identification of a few
prototypical parts [3, 1, 22, 23]. This identification of a limited number of
prototypical parts enables user understanding of the learned reasoning behind
the DL model output, which we show, achieves competitive performance against
its non-interpretable counterparts. The advantage of the case-based reasoning
process is that the explanations generated follow the same reasoning process a
specialist used to classify a kidney stone image.

However, for images of kidney stones, a great level of expertise is required,
and thus, explanations generated from PPs for kidney stone patches still
depend on non-obvious characteristics from the input image for non-experts,
reason for us to quantify the sensitivity of PPs to a set of perturbations
[19], which we dub herein as “Descriptors”.

Each perturbation is directed at a particular visual characteristic
(Brightness, Saturation, Shape, Texture, Contrast, and Hue). These descriptors
allow measuring individually the level of importance of visual features to
each of the PPs, as shown in the table in Fig. 2(b).

Our model shows the underlying decision-making process for each case, based on
several explanations which correspond to the most similar cases identified and
used by the model to give its suggested classification of the image. Fig. 2
shows an example of how our proposal works: I) Our model processes an input
image through a feature extractor backbone (CNN in Fig. 2(a)). II) Similarity
of the learned PPs to the extracted features from the CNN are calculated
accordingly to equation 1 (Similarity score in Fig. 2(a)). III) From the
Prototypes layer, in Fig. 2(a), the scores of PPs calculated against the input
image are processed through a linear weighted combination, a fully connected
layer, to provide a final classification of a kidney stone image. IV) The main
PPs that determined the classification are provided as the most similar cases
to their respective parts from the input image, as seen in the blue box in
Fig. 2(b). V) Finally, the saliency of certain characteristics for each
prototypical part is evaluated, the Descriptors, listed in the table in Fig.
2(b) and provide a description of the input image features that mainly
triggered the activation of the corresponding learned Prototypical-Part (PP).

This capacity yields predictions that are easy to understand by specialists
familiarized with the features related to each kidney stone type and subtype,
rendering classifications interpretable for urologists. Thus, aiming to help
in the morphological analysis of each kidney stone image, the trained model
provides the specialist with the means to trust the AI, check its output for
plausibility, and easily overrule it when necessary. Our approach, with its
inherently interpretable reasoning process, contrasts directly with previous
work that relied on post-hoc explanation techniques to explain a trained
black-box model on particular classifications [20, 11] or with global
explanations [8].

Following the architecture proposed on [3, 22], in this study, we explore
ProtoPNet architecture with some important differences, I) the analysis
performed explores different backbones, with the intention to compare CNN with
different levels of intricate connections. II) Different number of
prototypical parts per class and III) The relevance of data augmentation for
training. A relevant note is that the training of the models was performed
without part annotations, only relying on the class label.

Table 1: Number of complete endoscopic images acquired.

Subtype | Main component | Key | Surface | Section | Mixed
---|---|---|---|---|---
Ia | Whewellite | WW | 62 | 25 | 87
IIa | Weddellite | WD | 13 | 12 | 25
IIIa | Acide Urique | AU | 58 | 50 | 108
IVc | Struvite | STR | 43 | 24 | 67
IVd | Brushite | BRU | 23 | 4 | 27
Va | Cystine | CYS | 47 | 48 | 95
|  | TOTAL | 246 | 163 | 409

Figure 4: Examples of kidney stones patches extracted from the dataset used in
[8]. This endoscopic dataset consists of six of the most common kidney stones
namely Whewellite (WW), Weddellite (WD), Acide Urique (AU), Struvite (STR),
Brushite (BRU), and Cystine (CYS).

###  3.1 Kidney stone dataset

A simulated in-vivo dataset of endoscopic images was used for the experiments
reported in this contribution (Table 1). This dataset consists of surface (246
images) and section (163 images) views. Each view contains six of the most
common kidney stone sub-types: Ia (Whewellite, WW), IIa (Weddellite, WD), IIIa
(Uric Acid anhydrous, UA), IVc (Struvite, STR), IVd (Brushite, BRU), and Va
(Cystine, CYS). Two different reusable digital flexible ureteroscopes Karl
Storz (Storz Image 1 Hub and Storz image1 S) captured images of real kidney
stones in an environment that simulates real conditions (such as illumination
or blur) during surgery. In addition, the acquired endoscopic images were
labeled manually by the urologist Jonathan El-Beze (see [8] for more details).

In order to train DL-based models, 12,000 square-patches of 256$\times$256
pixels were extracted from 409 endoscopic images of both views. This combined
data is referred to as mixed dataset. Patches of the dataset are shown in Fig.
4. For more details, patches were extracted with the method described in [20]
and [15].

###  3.2 Model architecture:

The architecture implemented uses a feature extractor $f(x)$, in this proposal
a CNN, to perform the extraction of semantic features from the input images
$x$ and learn $n$ prototypes per class, corresponding to parts in the input
image. This first step of feature extraction from the input image implies the
potential diversity and quality of possible prototypical parts to be learned
from. Hence, three CNN backbones are explored (VGG16, ResNet50 and
DenseNet201) in order to compare the different possible performances obtained
with our approach and compared against pure non-interpretable CNNs present in
the estate of the art for kidney stone identification.

Next, two layers of $1\times 1$ convolutions are added, to adjust the depth of
the feature activation maps to the depth selected of 128 for the layer of PPs,
layer $g$. The variable $l$ is used to index each of the $128\times 1\times 1$
patches $f(x)_{l}$ across the spatial dimensions. The layer of PPs $g$,
contains $n$ prototypes to be learned per class, with dimensions $128\times
1\times 1$.

Since a prototype has the same number of channels but a smaller spatial
dimension than the convolutional feature maps $f(x)$, we can interpret the
prototype as representing a prototypical activation pattern of its class and
we can visualize the prototype as a patch of the training image it appears in.
The distance between $p_{j}$ and each of the $l$ $1\times 1$ spatial patches
of the convolutional feature map $f(x)$ is measured by $d_{j,l}=\lVert
p_{j}-f(x)_{l}\rVert_{2}^{2}$, and converts these distances to similarities
$s_{j,l}$ by:

$s_{j,l}=log\frac{d_{j,l}+1}{d_{j,l}+\epsilon}$ (1)

Where $\epsilon$ is a small value to avoid division by zero. This similarity
scores $s_{j,l}$ can be arranged spatially into a similarity map, upsampled
and superimposed on the input image to generate a saliency map as the ones
presented on the blue box in Fig. 2(b). The overall similarity score of the
corresponding PPs is the highest score $s_{j}$ obtained from $s_{j,l}$. These
similarity scores are then processed with a fully connected layer, as to
combine them in a weighted combination to obtain an output per class, which is
then normalized with a softmax function. The details for the initialization
and training hyperparameters, except for the data augmentation used, follow
[3].

###  3.3 Descriptors calculation:

Importance scores for image characteristics of the PPs are analyzed by
measurement of the change in the similarity score $s_{j}$ when perturbing the
input image with certain modifications $\hat{s_{j}}$. Let
$i\in\left\\{brightness,contrast,saturation,hue,shape,texture\right\\}$ denote
the type of modification, as further detailed in [19]. Then the local
importance score $\Phi_{local}^{i,j,k}$ of characteristic $i$ for test image
$j\in S_{test}$ with class $k$ on the $j-th$ prototype is the difference in
similarity score:

$\Phi_{local}^{i,j,k}=s_{j}-\hat{s_{j}}$ (2)

And its global importance is a weighted arithmetic mean by weighing the local
scores of all images in $S_{train}$ by their respective original similarity
score per PP:

$\Phi_{global}^{i,j}=\frac{\sum_{k=1}^{|S_{train}|}\Phi_{local}^{i,j,k}\cdot
s_{j}}{\sum_{k=1}^{|S_{train}|}s_{j}}$ (3)

Hence, if unmodified image $k$ gets a low similarity score with prototype $j$,
it will provide a low weight during the global importance calculation. In
contrast, if the prototype $j$ is clearly present in image $k$, the
architecture will assign a high similarity score and hence the local
importance score for image $k$ gets a high weight. As demonstrated in equation
4, these importance scores can be used to create global explanations that
explain a prototype, and local explanations, equation 2, that explain the main
visual features responsible for a high similarity score between a given image
and a prototype. As appreciated from equation 2 and equation 4, with this
descriptors the intention is to explain the identification of PPs on the input
image from the models perspective, not from human perception, nonetheless,
since the descriptors explored are humanly interpretable actually those do
enrich the explanations obtained. The specifics for the different
perturbations used to obtain the values of the PPs descriptors follow [19].

###  3.4 Evaluation metrics

To determine the different types of kidney stones in endoscopic images, MCA
specialists perform a visual inspection of the surface (external part) and
section (internal part) view, and provide a detailed report of each view of
the kidney stone. In addition, in order to improve performance in this task,
MCA also proposes combining information from surface and section views, since
both views provide relevant information for classification. [4]

In this work, we perform different training of mixed views models (similar to
how MCA specialists perform it). The accuracy and precision metrics are mainly
used to quantitatively evaluate different implementations of DL models
proposed in the SOTA mentioned to be intended as a CADx tool [2, 10, 15, 20,
27]. Therefore, our particular implementations were measured and compared on
accuracy as well. It can be observed in Table 2 the comparison of the accuracy
metric of recent works in the SOTA for the task of kidney stone identification
against our implementations.

##  4 Experimental Setup

Three different CNN architectures are explored for our proposal: i) a VGG16 to
examine the performance of a simple deep CNN, ii) a ResNet50 to consider a
medium size CNN with residual connections and iii) a DenseNet201 for the
evaluation of dense connections of a deep model to the SOTA for kidney stone
identification. These architectures were pre-trained on ImageNet, to then
train them with the kidney stone dataset. This enables having these CNNs as
baselines, and ResNet 50 and VGG16 are also used on other implementations in
the current SOTA for kidney stone identification [17, 11, 2, 14, 16, 29].

For the prototype layer $g$, the number of PPs per class explored are
$(1,3,10)$ PPs per each of the six classes of the three configurations of our
dataset, described in section 3.1. The training of our implemented ProtoPNet
models is a cycle consisting of three stages: ($A_{1}$) Updating the weights
of the convolutional layers used as a feature extractor $f(x)$ per 10 epochs,
using the loss function presented in [3]. ($A_{2}$) Projection of the
prototype layer $g$ to their corresponding closest parts after the
convolutional layers $f(x)$ updates. Detailed description in supplementary
material section S1.2 ($A_{3}$) Update of the fully connected layer $h_{1}$
weights for 20 epochs (FC & Softmax layer in Fig. 2(a)). We repeated this
cycle three times and selected the best model on the mixed test data obtained
during one of the ($A_{3}$) phases.

Preparation for the train and test sets were 80% (1,600 images per class) and
20% (400 images per class), respectively. The patches were also “whitened”
using the mean $m_{i}$ and standard deviation $\sigma_{i}$ of the color values
$I_{i}$ in each channel $(I^{w}_{i}=(I_{i}-m_{i}\sigma_{i})$, with $i=R,G,B)$.
The different mentioned configurations for training were performed on the
mixed data, described in Section 3.1, without data augmentation and with data
augmentation. The data augmentation used consists in randomly applying one
data augmentation at a time to an input image, from a set of visual
transformations. This set of transformations consists of horizontal flip,
vertical flip, rotation between -180 and 180 degrees, a distortion scale of
0.4, translation up to 0.2, or symmetric padding of 50 pixels. After selecting
the transformation it has a 50% chance of been applied.

After training the different model’s configurations, evaluation of the
perturbations, mentioned in Section 3.3, is carried out as PPs descriptors
and, additionally, we evaluate the performance of the models under these
perturbations (graph in the supplementary material as Fig. S1). To preserve
clarity of results, supplemantary Figures S1, S3 and S2, Fig. 6 and Tables 2
and 3, only show the results on the mixed test data. Additionally, the
performance of the trained models under different levels of change on the most
impactful perturbation, the hue channel perturbation, is evaluated. The
original change in hue explored in [19] applies a color jitter modification
with the torchvision library to the input images, fixed on a change of 0.1. We
explore a range of values for the change in the hue channel to characterize
the sensitivity of the performance of our models. Fig. S2 illustrates the
performance of one model training configuration under different hue channel
settings.

Figure 5:  UMAP of the Part-Prototypes and descriptors activations on the test
patches of kidney stone images, for the particular case of a ProtoPNet using
DenseNet201 as the backbone and 3 Prototypical parts per class. Our approach
allows obtaining separate clusters of the output classes.  Figure 6: The six
descriptor activations of the Prototypical Parts (PPs) per class, for the
particular case of a ProtoPNet using DenseNet201 as the backbone and 3
Prototypical parts per class.  Table 2:  Accuracy of state-of-the-art models
and three of ours for identifying kidney stones. Models trained with data
augmentation.

Method | Surface | Section | Mixed
---|---|---|---
Martinez, et al. [17] | 56.2±23.3 | 46.6±12.8 | 52.7±18.9
Estrade, et al. [11] | 73.7±17.9 | 78.8±10.6 | 70.1±22.3
Black, et al.[2] | 73.5±19.0 | 88.8±02.8 | 80.1±13.8
Lopez, et al.[14] | 81.0±03.0 | 88.0±02.3 | 85.0±03.0
Lopez, et al.[16] | 83.2±01.2 | 90.4±04.8 | 85.6±00.1
Villalvazo-Avila, et al.[29] | 88.8±02.8 | 84.4±06.0 | 96.6±00.5
Ours | Densenet201 3pps | 80.2±00.6 | 87.6±06.2 | 85.2±01.5
Resnet50 1pps | 80.9±01.3 | 91.0±03.3 | 88.2±02.1
Vgg16 10pps | 76.7±02.2 | 82.6±04.9 | 82.1±00.9

Table 3: Performance metrics of our models and convolutional neural networks used as backbones. Evaluations are divided between tests performed with the same type of data the model was trained, considered Independent Identically Distributed (IID) data and tests performed under perturbations, which are considered Out Of Distribution (OOD) data. All models were pre-trained on ImageNet.  |  |  | IID |  | OOD
---|---|---|---|---|---
Model | Accuracy (%) | Precision (%) | F1 (%) |  | Accuracy (%) | Precision (%) | F1 (%)
DenseNet201 | 89.67±3.60 | 90.51±3.60 | 89.47±3.46 |  | 58.44±28.06 | 60.22±29.62 | 54.54±32.73
ResNet50 | 83.40±5.32 | 85.58±4.30 | 83.11±5.21 |  | 45.12±20.82 | 49.16±24.08 | 38.25±23.47
Vgg16 | 81.47±2.72 | 84.71±2.64 | 80.62±3.84 |  | 58.35±22.62 | 63.73±19.93 | 54.38±26.30
CNNs Average | 84.84±5.29 | 86.93±4.35 | 84.40±5.60 |  | 53.97±24.62 | 57.70±25.40 | 49.06±28.56
DenseNet201 | 1 PPs | 86.29±1.91 | 87.54±1.34 | 86.25±1.78 |  | 60.31±20.32 | 70.29±14.38 | 57.45±23.23
3 PPs | 85.19±1.50 | 86.08±1.27 | 85.21±1.42 |  | 59.97±19.91 | 69.21±12.59 | 57.66±22.14
10 PPs | 87.29±0.92 | 87.99±0.94 | 87.23±0.92 |  | 61.20±19.33 | 71.26±11.92 | 59.01±21.37
ResNet50 | 1 PPs | 88.21±2.07 | 88.61±1.88 | 88.21±2.05 |  | 58.26±24.23 | 66.48±18.53 | 55.16±27.50
3 PPs | 86.66±1.37 | 87.11±1.55 | 86.62±1.45 |  | 57.90±21.82 | 63.24±19.16 | 54.21±25.32
10 PPs | 85.44±1.44 | 86.16±1.15 | 85.43±1.42 |  | 56.31±22.66 | 64.81±16.55 | 53.77±25.52
Vgg16 | 1 PPs | 81.78±1.60 | 83.38±1.94 | 81.82±1.62 |  | 51.78±19.02 | 53.82±19.40 | 46.63±22.77
3 PPs | 82.21±3.33 | 82.92±3.82 | 81.84±3.75 |  | 51.23±19.40 | 54.09±19.80 | 45.88±23.20
10 PPs | 82.08±0.90 | 83.39±0.92 | 82.13±0.87 |  | 54.69±18.05 | 56.48±17.56 | 49.89±21.51
ProtoPNets Average | 85.02±2.83 | 85.91±2.66 | 84.97±2.89 |  | 56.85±20.63 | 63.30±17.90 | 53.30±23.83

##  5 Results and Discussion

###  5.1 Implementation details

Experiments were carried out to explore different backbones for the feature
extractor network $f(x)$, the numbers of PPs per class and the significance of
data augmentation. This, is to evaluate the model’s performance and the
obtained explanations. In contrast to CNNs, which are black-box classifiers,
our proposal provides explanations for input images in three different ways:
i) by showing the activation area, where each PP is identified. ii) The
corresponding representative example image to each activated PPs. iii)
Quantitative evaluations of how relevant six different visual descriptors are
for the activation of each PP.

###  5.2 Visualization of a case example of explanations

In Fig. 5, an example of our model classifications is shown. We found that
averaging the global scores of the PPs descriptors per class results in a
hierarchy of descriptor significance, as depicted in Fig. 6. The WD class
reveals a unique pattern of descriptor activation. An examination of several
instances shows that the average activation of a class descriptors can hint at
misclassifications. As shown in Fig. 5, a green diamond represents a
misclassified AU image, which is misidentified as WD due to its activations
corresponding to those of the WD class. This causes the image to appear within
the WD class cluster on the UMAP plot. We use a UMAP visualization to present
PPs activations for input images in the context of the top three
discriminative dimensions. This UMAP highlights class separability for each
output class of the ProtoPNet, as showcased in Fig. 5. The use of descriptors
mitigates the cases for visually similar PPs by providing details on the
characteristics most relevant for each PP [19]. Additionally, learning PPs
from patches of kidney stone images allows observing PPs are learned
proportionally to the architecture of the CNN used.

###  5.3 Performance results

A summary of the accuracy of the different configurations of models trained is
displayed in Table 3. Tests carried out on the original mixed test dataset,
corresponding to the training mixed dataset, are considered Independent
Identically Distributed (IID) data, tests performed under perturbations
applied to the original test dataset are considered Out Of Distribution (OOD)
data and named accordingly in Table 3. Usage of data augmentation yielded an
average 3% improvement accuracy compared to the non-augmented approach
consistently across different evaluations. To ease presentation, we report
results just on models trained with data augmentation.

The performance of our implemented ProtoPNet models is comparable with its
corresponding non-interpretable CNN backbone models as also appreciated in
Table 3. The average accuracy of our implementations keeps competitiveness
with a difference not bigger than $\leq 0.18\%$, when compared with the
average accuracy of the baseline CNN models. The mean values and standard
deviations displayed in Table 3 are computed based on the configurations of
the models outlined therein. It is worth noticing that Table 3 shows our
interpretable models to outperform traditional CNNs from the SOTA on kidney
stone classification when evaluated on perturbed input images. Specifically,
our models achieve 1.9% higher accuracy, 5.7% higher precision, and 3.2%
higher F1 scores compared to traditional models.

##  6 Conclusions and Future Directions

We showed that by adapting and fine tunning CNN models into ProtoPNets is
possible to convert considered black-box CNN models into interpretable ones,
bringing detailed and faithful explanations. The main visual features in the
input image that produce the activation of PPs in particular and in the
average cases per class and per dataset are able to be extracted from the
evaluated behavior of the network, as shown in Fig. 6. This case-based
reasoning approach allows to generate explanations and additional activation
details enabling urologists to use these models as assistance tools for the
MCA. More importantly, we showed evidence for how to use the overall behavior
of descriptors per class in combination with UMAP projection of the behavior
of the model to identify the reasons behind misclassification, as shown in
Fig. 5. Also, the UMAP visualization illustrates the blending of certain
classes that share similar color and texture characteristics. To enhance
cluster separability a metric learning approach could be explored. This to
penalizes close PPs from different classes, encourages proximity for PPs of
the same class but avoids their collapse.

Upon examining Table 3, on one hand, it is evident that despite a minor
decrease in average performance for IID evaluations in PPs models compared to
their corresponding CNNs. This could be due to restricting the classification
to a few PPs per class, in comparison CNN models use a few hundred features.
On the other hand, PPs models consistently exhibit reduced variance in
performance metrics. Additionally, these models consistently perform better in
OOD evaluations, suggesting that PPs models tend to be more robust than their
CNN counterparts. This characteristic could be attributed to a regularization
effect resulting from measuring the similarity between the learned PPs and the
image features. PPs models inherit most of the performance characteristics of
CNNs, indicating that methods to enhance CNN performance should also improve
PPs models, thus expanding potential improvements.

In future research, the model’s reasoning and reporting mechanisms could be
refined to align more closely with the structured lexicon employed by
urologists. This could be achieved by incorporating natural language
descriptions for the PPs identified and their descriptors to make explanations
more accessible to specialists and a broader range of users.

## Compliance with Ethical Standards

The images were captured in medical procedures following the ethical
principles outlined in the Helsinki Declaration of 1975, as revised in 2000,
with the consent of the patients.

## Acknowledgments

The authors extend their gratitude to the AI Hub and the Centro de Innovación
de Internet de las Cosas at Tecnológico de Monterrey for generously providing
access to an NVIDIA DGX computer, which facilitated the experiments conducted
in this paper. Additionally, we would like to express our gratitude to CONACyT
for the Ph.D. scholarships for Daniel Flores-Araiza and Francisco López-Tiro.

## References

  * [1] Alina Barnett, Fides Regina Schwartz, Chaofan Tao, Chaofan Chen, Yinhao Ren, Joseph Y. Lo, and Cynthia Rudin.  Iaia-bl: A case-based interpretable deep learning model for classification of mass lesions in digital mammography.  Nat. Mach. Intell., 3:1061–1070, 2021.
  * [2] Kristian M Black, Hei Law, Ali Aldoukhi, Jia Deng, and Khurshid R Ghani.  Deep learning computer vision algorithm for detecting kidney stone composition.  2020\.
  * [3] Chaofan Chen, Oscar Li, Chaofan Tao, Alina Jade Barnett, Jonathan Su, and Cynthia Rudin.  This Looks Like That: Deep Learning for Interpretable Image Recognition.  CoRR, jun 2018.
  * [4] Mariela Corrales, Steeve Doizi, Yazeed Barghouthy, Olivier Traxer, and Michel Daudon.  Classification of stones according to michel daudon: a narrative review.  European Urology Focus, 7(1):13–21, 2021.
  * [5] Michel Daudon, Arnaud Dessombz, Vincent Frochot, Emmanuel Letavernier, Jean-Philippe Haymann, Paul Jungers, and Dominique Bazin.  Comprehensive morpho-constitutional analysis of urinary stones improves etiological diagnosis and therapeutic strategy of nephrolithiasis.  Comptes Rendus Chimie, 19(11-12):1470–1491, 2016.
  * [6] Michel Daudon, Paul Jungers, Dominique Bazin, and James C Williams.  Recurrence rates of urinary calculi according to stone composition and morphology.  Urolithiasis, 46(5):459–470, 2018.
  * [7] Brian Edwards.  Fda guidance on clinical decision support: Peering inside the black box of algorithmic intelligence, Nov 2019.
  * [8] Jonathan El Beze, Charles Mazeaud, Christian Daul, Gilberto Ochoa-Ruiz, Michel Daudon, Pascal Eschwège, and Jacques Hubert.  Evaluation and understanding of automated urolithiasis recognition methods.  BJU International, 2022.
  * [9] Laurence Estepa and Michel Daudon.  Contribution of fourier transform infrared spectroscopy to the identification of urinary stones and kidney crystal deposits.  Biospectroscopy, 3(5):347–369, 1997.
  * [10] Vincent Estrade, Michel Daudon, Emmanuel Richard, Jean-Christophe Bernhard, Franck Bladou, Gregoire Robert, and Baudouin Denis de Senneville.  Towards automatic recognition of pure & mixed stones using intraoperative endoscopic digital images.  BJU International, abs/2105.10686, 2021.
  * [11] Vincent Estrade, Michel Daudon, Emmanuel Richard, Jean-Christophe Bernhard, Franck Bladou, Gregoire Robert, and Baudouin Denis de Senneville.  Towards automatic recognition of pure and mixed stones using intra-operative endoscopic digital images.  British Journal of Urology (BJU) International, 2021.
  * [12] Daniel Flores-Araiza, Francisco Lopez-Tiro, Elias Villalvazo-Avila, Jonathan El-Beze, Jacques Hubert, Gilberto Ochoa-Ruiz, and Cristian Daul.  Interpretable deep learning classifier by detection of prototypical parts on kidney stones images.  arXiv preprint arXiv: Arxiv-2206.00252, 2022.
  * [13] Justin I Friedlander, Jodi A Antonelli, and Margaret S Pearle.  Diet: from food to stone.  World journal of urology, 33(2):179–185, 2015.
  * [14] Francisco Lopez, Andres Varelo, Oscar Hinojosa, Mauricio Mendez, Dinh-Hoan Trinh, Yonathan Elbeze, Jacques Hubert, Vincent Estrade, Miguel Gonzalez, Gilberto Ochoa, and Christian Daul.  Assessing deep learning methods for the identification of kidney stones in endoscopic images.  In 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society, EMBC 2021, Mexico, Mexico, Nov. 2021.
  * [15] Francisco Lopez, Andres Varelo, Oscar Hinojosa, Mauricio Mendez, Dinh-Hoan Trinh, Yonathan ElBeze, Jacques Hubert, Vincent Estrade, Miguel Gonzalez, Gilberto Ochoa, et al.  Assessing deep learning methods for the identification of kidney stones in endoscopic images.  In 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pages 2778–2781. IEEE, 2021\.
  * [16] Francisco Lopez-Tiro, Juan Pablo Betancur-Rengifo, Arturo Ruiz-Sanchez, Ivan Reyes-Amezcua, Jonathan El-Beze, Jacques Hubert, Michel Daudon, Gilberto Ochoa-Ruiz, and Christian Daul.  Boosting kidney stone identification in endoscopic images using two-step transfer learning.  arXiv preprint arXiv: Arxiv-2210.13654, 2022.
  * [17] Adriana Martínez, Dinh-Hoan Trinh, Jonathan El Beze, Jacques Hubert, Pascal Eschwege, Vincent Estrade, Lina Aguilar, Christian Daul, and Gilberto Ochoa.  Towards an automated classification method for ureteroscopic kidney stone images using ensemble learning.  In 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pages 1936–1939. IEEE, 2020\.
  * [18] Anmar Nassir, Hesham Saada, Taghreed Alnajjar, Jomanah Nasser, Waed Jameel, Soha Elmorsy, and Hattan Badr.  The impact of stone composition on renal function.  Urology Annals, 10(2):215, 2018.
  * [19] Meike Nauta, Annemarie Jutte, Jesper C. Provoost, and Christin Seifert.  This looks like that, because … explaining prototypes for interpretable image recognition.  CoRR, abs/2011.02863, 2020.
  * [20] Gilberto Ochoa-Ruiz, Vincent Estrade, Francisco Lopez, Daniel Flores-Araiza, Jonathan El Beze, Dinh-Hoan Trinh, Miguel Gonzalez-Mendoza, Pascal Eschwège, Jacques Hubert, and Christian Daul.  On the in vivo recognition of kidney stones using machine learning.  arXiv preprint arXiv:2201.08865, 2022.
  * [21] Ege Gungor Onal and Hakan Tekgul.  Assessing kidney stone composition using smartphone microscopy and deep neural networks.  BJUI Compass, 3:310 – 315, 2022.
  * [22] Dawid Rymarczyk, Lukasz Struski, Jacek Tabor, and Bartosz Zielinski.  Protopshare: Prototype sharing for interpretable image classification and similarity discovery.  CoRR, abs/2011.14340, 2020.
  * [23] Dawid Rymarczyk, Łukasz Struski, Michał Górszczak, Koryna Lewandowska, Jacek Tabor, and Bartosz Zieliński.  Interpretable image classification with differentiable prototypes assignment.  arXiv preprint arXiv: Arxiv-2112.02902, 2021.
  * [24] Ramprasaath R. Selvaraju, Abhishek Das, Ramakrishna Vedantam, Michael Cogswell, Devi Parikh, and Dhruv Batra.  Grad-cam: Why did you say that? visual explanations from deep networks via gradient-based localization.  CoRR, abs/1610.02391, 2016.
  * [25] R Siener and A Hesse.  Fluid intake and epidemiology of urolithiasis.  European journal of clinical nutrition, 57(2):S47–S51, 2003.
  * [26] Shelly Soffer, Avi Ben-Cohen, Orit Shimon, Michal Marianne Amitai, Hayit Greenspan, and Eyal Klang.  Convolutional neural networks for radiologic images: A radiologist’s guide.  Epub, pages 590–606, 2019.
  * [27] Alejandro Torrell Amado.  Metric learning for kidney stone classification.  2018\.
  * [28] Adie Viljoen, Rabia Chaudhry, and John Bycroft.  Renal stones.  Annals of clinical biochemistry, 56(1):15–27, 2019.
  * [29] Elias Villalvazo-Avila, Francisco Javier Lopez-Tiro, Jonathan El Beze, Jacques Hubert, Miguel González-Mendoza, Gilberto Ochoa-Ruiz, and Christian Daul.  Improved kidney stone recognition through attention and multi-view feature fusion strategies.  ArXiv, abs/2211.02967, 2022.
  * [30] Heliang Zheng, Jianlong Fu, Tao Mei, and Jiebo Luo.  Learning multi-attention convolutional neural network for fine-grained image recognition.  In Proceedings of the IEEE International Conference on Computer Vision (ICCV), Oct 2017.

##  S1 Supplementary Material

###  S1.1 Performance under visual perturbations

Accuracy evaluated on the mixed test dataset, for the different ProtoPNet
backbones and the number of PPs learned was obtained to identify the level of
loss in performance of the models learning PPs. With this evidence for the
kidney stone classification task made available, future users will have a
reference for which visual characteristics PPs models tend to be more
susceptible to and to which degree. It is found that hue perturbation on the
input images is on average the perturbation with the most dramatic loss in
performance for models learning PPs training with the Kidney Stones datasets.
In contrast, this type of architecture tends to be robust against moderated
changes in the brightness and saturation of the input images. With this
evidence, a call for precaution when dealing with assumed small perturbation
on the input images [19] is implemented.

Also, the performance of the traditional CNN used as feature extractors were
explored under the same perturbations, these results are shown in the first
four rows under the OOD performance metrics in Table 3 for each of the CNN
backbones explored and their performance for IID test data and OOD is
summarized in Fig. 3(a). These same evaluations are shown for the different
backbones of PPs models trained in Fig. 3(b). It is observed PPs models
present a lower standard deviation when compared to their CNN counterparts for
IID evaluations, per this behavior for each of the different backbones used by
the models, as seen in Fig. 3(c).

Figure S1: Accuracy evaluated on the test dataset, averaging the different
ProtoPNet backbones and the number of Prototypical parts learned. It is found
that hue perturbation on the input images is on average the perturbation with
the most dramatic loss in performance for ProtoPNet models training with the
Kidney Stones dataset. In contrast, the ProtoPNet architecture tends to be
robust against moderated changes in the brightness and saturation of the input
images.  Figure S2: Accuracy evaluated on mixed test dataset under
modifications of the images hue channel, with values between -0.1 and 0.1.
These evaluations reflect the level of robustness ProtoPNet achieves with a
DenseNet201 CNN backbone with 3 part-prototypes per class, under the
perturbation that most affect the model.

(a)

(b)

(c)

Figure S3:  (a) Accuracy of CNN models evaluated separated on Independent and
Identically Distributed (IID) test data and Out Of Distribution (OOD) tests.
(b) Accuracy of ProtoPNet models evaluated on IDD and OOD test data. (c) Same
as in ”b” for each of the three backbones used.

###  S1.2 Prototypes projection details

The Projection of Prototypical parts (PPs) is the intermediate step $A_{2}$,
mentioned in Section 4, performed in the training process that allows
visualization of the learned PPs. In this step, each prototype $p_{j}$ is
assigned the value of its nearest latent training patch $f(x)_{l}$ from all
the images of the same class $k$ initially assigned to $p_{j}$. Therefore,
this distance is $d_{j,l,k}=\lVert p_{j,k}-f(x)_{l,k}\rVert_{2}$. In this way,
each abstract PPs learned conceptually equates to one training image patch.
Allowing the faithful visualization of the learned PPs, use to generate the
visual explanations as well as the final output classification given by the
model. Mathematically, for the projection prototype $p_{j}$ of class $k$,
i.e., $p_{j}\in P_{k}$, we perform the following update:

$p_{j,k}\xleftarrow{}\operatorname*{arg\,min}_{f(x)_{l,k}\in S_{train}}\lVert
p_{j,k}-f(x)_{l,k}\rVert_{2}$ (4)

where $S_{train}=\\{f(x)_{l,k}:f(x)_{l,k}\in f(x_{k})$ for all training images
$S_{train}$, $(x,y):y\in k\\}$.