# Throwing Objects into A Moving Basket While Avoiding Obstacles

Hamidreza Kasaei1 and Mohammadreza Kasaei2 1Hamidreza Kasaei is with the
Department of Artificial Intelligence, Bernoulli Institute, Faculty of Science
and Engineering, University of Groningen, The Netherlands. Email:
<EMAIL_ADDRESS>Mohammadreza Kasaei is with the School of
Informatics, University of Edinburgh, UK. Email<EMAIL_ADDRESS>
###### Abstract

The capabilities of a robot will be increased significantly by exploiting
throwing behavior. In particular, throwing will enable robots to rapidly place
the object into the target basket, located outside its feasible kinematic
space, without traveling to the desired location. In previous approaches, the
robot often learned a parameterized throwing kernel through analytical
approaches, imitation learning, or hand-coding. There are many situations in
which such approaches do not work/generalize well due to various object
shapes, heterogeneous mass distribution, and also obstacles that might be
presented in the environment. It is obvious that a method is needed to
modulate the throwing kernel through its meta-parameters. In this paper, we
tackle object throwing problem through a deep reinforcement learning approach
that enables robots to precisely throw objects into moving baskets while there
are obstacles obstructing the path. To the best of our knowledge, we are the
first group that addresses throwing objects with obstacle avoidance. Such a
throwing skill not only increases the physical reachability of a robot arm but
also improves the execution time. In particular, the robot detects the pose of
the target object, basket, and obstacle at each time step, predicts the proper
grasp configuration for the target object, and then infers appropriate
parameters to throw the object into the basket. Due to safety constraints, we
develop a simulation environment in Gazebo to train the robot and then use the
learned policy in real-robot directly. To assess the performers of the
proposed approach, we perform extensive sets of experiments in both simulation
and real robots in three scenarios. Experimental results showed that the robot
could precisely throw a target object into the basket outside its kinematic
range and generalize well to new locations and objects without colliding with
obstacles. The video of our experiments can be found at
https://youtu.be/VmIFF__c_84

##  I Introduction

Almost all humans are familiar with the ability to throw objects, as we learn
how to throw a ball during a game (e.g., basketball) or an object into a bin
(e.g., tossing dirty clothes into the laundry basket). We throw objects either
to speed up tasks by reducing the time of pick-and-place or to place them in
an unreachable place [1]. Therefore, adding such a throwing motion to a
robotic manipulator would enhance its functionality too. In particular,
throwing object is a great way to use dynamics and increase the power of a
robot by enabling it to quickly place objects into the target locations
outside of the robot’s kinematic range. However, the act of precisely throwing
is actually far more complex than it appears and requires a lot of practice
since it depends on many factors, ranging from pre-throw conditions (e.g.
initial pose of the object inside the gripper) to the physical properties of
the object (e.g. shape, size, softness, mass, material, etc.). Many of these
elements are challenging to describe or measure analytically, hence, earlier
research has frequently been limited to assuming predefined objects (e.g.,
ball) and initial conditions (e.g., manually placing objects in a per-defined
location). When obstacles are present in the environment and the target basket
is moving, throwing becomes even more difficult. To the best of our knowledge,
we are the first group to address such a challenging object throwing problem.

To accomplish the throwing task successfully, a robot must process the visual
information to realize which objects exist in the scene (i.e., target object,
basket, and obstacles), what are the state of the objects (i.e., pose, speed,
etc.), and how to grasp the target object (grasp synthesis). The robot then
finds an obstacle-free trajectory to grasp the object. Afterward, given the
obstacles that exist in the scene and the state of the target basket, it needs
to predict throwing parameters to throw the object to the desired location
precisely (e.g., the velocity of executing the throw trajectory, time of
release, etc.). Lastly, the robot executes the throwing motion using those
parameters.

Figure 1: An example scenario of throwing an object into a moving basket
located outside of the robot’s workspace while an obstacle obstructing the
path. To accomplish this task successfully, the robot should perceive the
environment through its RGB-D camera, and then infer the proper parameters to
throw the object into the basket.

In this paper, we formulate object throwing as a RL problem to enable the
robot to generalize well across a variety of objects and react quickly to
dynamic environments (i.e., moving basket). For RL, the exploration phase is
often unsafe in the real-world. It takes a while to build up enough experience
to train the policy to function successfully in a dynamic environment with
moving targets and obstacles. Therefore, we develop a simulation in Gazebo,
very similar to our real-robot setup, and train the robot in Gazebo initially.
Afterwards, the learned policy is used in real-robot settings directly. We
extensively evaluate the performance of our approach in both simulation and
real-robot using three different tasks with ascending levels of difficulties.
Experimental results show that the proposed method produces throws that are
more accurate than baseline alternatives. In summary, our key contributions
are threefold:

  * •

To the best of our knowledge, we are the first group that addresses object
tossing while obstacles are present in the environment and the target basket
is moving.

  * •

Despite only trained using simulation data, the proposed approach can be
directly applied to real-robot. Furthermore, it shows impressive
generalization capability to new target locations and unseen objects.

  * •

Our experiments show that the trained policy could achieve above $80\%$ object
throwing accuracy for the most difficult task (i.e., throwing object into the
basket while there is an obstacle obstructing the path) in both simulation and
real robot environments.

##  II Related Work

The robotics community has long been interested in giving service robots the
ability to throw objects [2, 3, 4, 5, 6]. Throwing formulae were mostly
influenced by analytical models in the late 1990s and early 2000s [7], while
such formulations are increasingly moving toward learning approaches today [8,
4]. In the following subsections, we briefly review these approaches.

###  II-A Analytical Approaches

Earlier throwing systems relied on handcrafting or mechanical analysis and
then optimizing control parameters to execute a throw such that the projectile
(typically a ball) lands at a target location. As we previously highlighted,
precisely modeling of dynamics is difficult because it calls for knowledge
about the physical characteristics of the object, gripper and environment,
which are hard to quantify [7]. For instance, Y. Gai et al., derived an
analytical approach for throwing a ball using a manipulator with a single
flexible link through Hamilton’s principle [3]. This is an example of tuning
for a single object, a ball in this case. In another work, Jwu-Sheng Hu et
al., [2] discussed a stereo vision system for throwing a ball into a basket.
They calculated the ball-throwing transformation for a specific ball object
based on cubic polynomial. In [9], an analytical approach is used to predict
the end-effector velocity (magnitude and direction) as well as a duration
movement for underhand throwing task by a humanoid robot. Such approaches to
some extend work for specific scenario but have difficulties generalizing over
changing dynamics and various objects.

###  II-B Learning Approaches

Unlike analytical approaches for throwing, learning-based methods enable
robots to learn/optimize the main task directly through success or failure
signals. In general, learning-based throwing approaches demonstrate better
performance than analytical methods[10, 11]. In [10], a deep predictive policy
training architecture (DPPT) is presented to teach a PR2 robot object-grasping
and ball-throwing tasks. They showed DPPT is successful in both simulated and
real robots. In another work, Kober et al. [11] introduced an RL-based method
for dart throwing task based on a kernelized version of the reward-weighted
regression. In both of these works, the properties of the object (ball and
dart) are known a-priori. In contrast to both of these approaches, we do not
make assumptions about the physical properties of objects that are thrown.

In some other works, researchers tried to combine the potential of analytical
and learning approaches for robotic throwing tasks. In particular, analytical
models are used to approximate the initial control parameters, and a learning-
based model is used to estimate residual parameters to adjust the initial
parameters. Such approaches are called residual physics. For instance, [4]
proposed TossingBot, an end-to-end self-supervised learning method for
learning to throw arbitrary objects with residual physics. Similar to our
work, their approach was able to throw an object into a basket. Unlike our
approach, they used an analytical approach for estimating initial control
parameters, and then used an end-to-end formulation for learning residual
velocity for throwing motion primitives. We formulate the throwing task as an
RL problem that modulates the parameters of a kernel motion generator. In
contrast to all reviewed works, our formulation allows the robot to throw the
object into a moving basket while avoiding present obstacles, whereas, in all
reviewed works, the throwing task is considered in an obstacle-free
environment where the target is static and known in advance.

##  III Method

In this section, the preliminaries are briefly reviewed, followed by a
discussion of how we formulate object throwing as an RL problem. The
perception that represents the world model at each time step is the subject of
the last subsection.

###  III-A Preliminaries

Markov Decision Process (MDP): An MDP can be described as a tuple containing
four basic elements:
$(s_{t},a_{t},p(s_{t+1}|s_{t},a_{t}),r(s_{t+1}|s_{t},a_{t}))$, where the
$s_{t}$ and $a_{t}$ are the continuous state and action at time step $t$,
respectively. $p(s_{t+1}|s_{t},a_{t})$ shows the transition probability
function to reach to the next state $s_{t+1}$ given the current state $s_{t}$
and action $a_{t}$. The $r(s_{t+1}|s_{t},a_{t})$ denotes the immediate reward
received from the environment after the state transition.

Off-policy RL: In online RL, an agent continuously interacts with the
environment to accumulate experiences for learning the optimal
policy$\pi^{*}$. The agent seeks to maximize the expected future return
$R_{t}=\mathbb{E}[\sum_{i=t}^{\infty}\gamma^{i-t}r_{i+1}]$ with a discounted
factor $\gamma\in[0,1]$ weighting the future importance. The expected return
under a policy $\pi$ after taking action $a$ in the state $s$ is computed by a
corresponding action-value function
$Q^{\pi}(s,a)=\mathbb{E}[R_{t}|s_{t}=s,a_{t}=a]$. By following policy $\pi$,
we can compute $Q^{\pi}$ through the Bellman equation:

$\begin{split}\scalebox{0.87}{\mbox{$\displaystyle
Q^{\pi}(s_{t},a_{t})=\mathbb{E}_{s_{t+1}\sim
p}[r(s_{t},a_{t})+\gamma\mathbb{E}_{a_{t+1}\sim{A}}[Q_{\pi}(s_{t+1},a_{t+1})]],$}}\end{split}$
(1)

where $A$ states the action space. Consider $Q^{*}(s,a)$ is the optimal
action-value function. RL algorithms aim to find an optimal policy $\pi^{*}$
such that $Q^{\pi^{*}}(s,a)=Q^{*}(s,a)$ for all states and actions.

###  III-B Problem Formulation

Given the start pose $p_{s}$ (i.e., grasp synthesis of an object) and the goal
pose $p_{g}$ (i.e., pose of the basket), the throwing task is defined as the
problem of finding proper parameters to throw the object into the basket while
avoiding obstacles $o$ in the environment. A fully observable MDP can be used
to represent this task, and the off-policy reinforcement learning framework
can be used to solve it. We discuss the detailed RL formulation of throwing
task in the following subsections.

####  III-B1 States

A feature vector is used to describe the continuous state, including the
robot’s proprioception, the pose of the obstacle and the goal in the
environment, releasing time, duration of trajectory execution, and the
distance between the thrown object and the goal. In particular, we formed a
kernel trajectory for throwing an object in a straight way using trial and
error. We then let the robot learns the initial and final value for the
shoulder joint, speed of executing the trajectory, and releasing time in the
learning process.

In particular, at each step $t$, we record the initial and final shoulder
joint values ($j_{i}$ and $j_{f}$) in radians as the proprioception: ${\rm
proprio}=(j_{i},j_{f})\in\mathbb{R}^{2}$. Then, we estimate the obstacle’s
position in task space and describe it as ${\rm
obs}=(x^{o},y^{o},z^{o})\in\mathbb{R}^{3}$. We describe the position of goal
(i.e, center of the box) as a point in the task space: ${\rm
goal}=(x^{g},y^{g},z^{g},\dot{x}^{g},\dot{y}^{g},\dot{z}^{g})\in\mathbb{R}^{6}$.
We also consider the absolute distance of the thrown object relative to the
obstacle and goal, and also the distances in the $X$ and $Y$ axes, ${\rm
dist}=(d^{g},d^{g}_{x},d^{g}_{y},d^{o},d^{o}_{x},d^{o}_{y})\in\mathbb{R}^{6}$.
We also record two timing profiles, the duration of executing the throwing
trajectory $\tau$, and the time for releasing the object, $t_{r}$, where
$t_{r}<\tau$ and ${\rm time}=(t_{r},\tau)\in\mathbb{R}^{2}$. Finally, the
state feature vector can be represented as: $s=({\rm proprio},{\rm obs},{\rm
goal},{\rm dist},{\rm time})\in\mathbb{R}^{19}$.

####  III-B2 Actions

Each action is denoted by a vector $a\in([-1,1])^{4}$, which represents (i)
the initial and (ii) the final shoulder joint values, (iii) the duration of
trajectory execution, and (iv) the releasing object time.

####  III-B3 Transition function

In each training episode, we set the pose of the goal randomly and then, set
the pose of the obstacle between the robot and the goal pose with $\pm 5cm$
randomness on the x-axis. Therefore, the transition function is determined by
executing the throwing trajectory given the sampled parameters. Specifically,
the next state, $s_{i+1}$, can be computed after executing the action $f_{s}$;
i.e., $s_{i+1}=f_{s}(s_{i},a_{i})$, where $s_{i}$ and $a_{i}$ are the current
state and action respectively. It should be noted that since the transition
function is unknown, our off-policy reinforcement learning framework is also
model-free.

####  III-B4 Rewards

A success is reached if the thrown object falls into the target basket after
executing the throwing action. In particular, we calculate the absolute
distance between the object and the goal ${\rm dis}({\rm o},g)<{\rm d}$, where
${\rm d}$ denotes the radius of a cylindrical space that is fitted inside the
basket.

As a result of the current state and the action taken, if the thrown object
collides with the obstacle, severe punishment is given by a negative reward
$r=-10$. It should be noted that the collision information can only be
obtained after the action has been executed. If the next state results in
success (the thrown object falls into the basket), we encourage such behavior
by setting the reward $r=1$. In the case of throwing the object outside the
target basket is also penalized by calculating rewards based on distance
$r=-~{}{\rm dis}({\rm o},g)$. An episode is terminated after executing a
throwing action. It should be noted that successful attempts are recorded for
later use in behavioral cloning.

###  III-C Perception

Figure 2: Visualizing the output of our perception module for an example
scene: it provides world model information in terms of object pose, size, and
label. In particular, the estimated object’s pose is shown by a reference
frame, the object’s size is demonstrated by a green bounding box, and the
label of the object is highlighted on top of the object’s Z axis by red.
Moreover, the predicted grasp pose for the target object is shown by a yellow
gripper.  Figure 3: Our experimental setups in (left) simulation and (right)
real-robot settings: our dual-arm robot consists of two UR5e arms for
manipulating objects, and an Asus Xiton RGB-D camera to perceive the
environment. We have developed a simulation environment in Gazebo very similar
to our real robot to reduce the gap between the simulation and the real robot
and facilitate transfer learning. Due to safety constraints, we initially
trained the robot in the simulation environment, and then directly transferred
the learned policy to the real robot without fine-tuning.

Due to safety considerations, we trained the proposed throwing policy in
simulation first and then use the learned policy on our real-robot platform.
In the case of simulation, we developed an interface that provides all
necessary information based on Gazebo’s services, while for real-robot
experiments, we exploit our robotic setup to detect and track the pose of the
objects, recognize them, and predict stable grasp syntheses for each of the
object in 3D space [12][13][14]. In particular, our real robot uses an RGB-D
Asus Xtion camera to perceive the environment by capturing point cloud data at
$30$ Hz. Generally, a point cloud consists of a set of points,
$p_{i}:i\in\\{1,\dots,n\\}$, where each point is described by its 3D
coordinates $[x,y,z]$ and RGB-D information. To track the pose and speed of
the target object, we use a particle filter that considers shape and color
data [15]. It is worth mentioning that for this work, we force the system to
grasp the object from above near the center of mass. In particular, our
perception system provides a world model service that the agent can call at
each time step to receive the current state of the world, which includes the
unique ID, pose, speed, label, and grasp synthesis of each object. Figure 2
shows the outputs of our perception system regarding object detection (i.e.,
highlighted by green bounding boxes and reference frames), object recognition
(i.e., highlighted by red on top of each object), and pre-grasp pose (i.e.,
highlighted by yellow gripper). For more information about our perception and
grasping pipelines, please refer to our earlier works [12][14].

##  IV Experiments and Results

We performed multiple rounds of experiments in both simulation and real-robot
settings to validate our method. In this work, we evaluate the performance of
the proposed approach based on the throwing success rate, which is calculated
as the number of times a thrown object lands into the target box divided by
the number of trials. More specifically, we tried to investigate the following
questions: (1) Which of the RL approach outperforms other baselines in terms
of throwing success rate when used in the same environment and task settings?
(2) Does our method learn to safely throw an object into a target basket while
there are obstacles obstructing its path? (3) Can the policy learned in
simulation transfers well to our real-robot where noise and uncertainty exist?

###  IV-A Experimental setup and tasks settings

Our experimental setups in simulation and real-robot are depicted in Fig. 3.
In particular, we developed a simulation environment in Gazebo utilizing an
ODE physics engine, which is very similar to our real-robot setup. Our setup
consists of an Asus xtion camera, two Universal Robots (UR5e) equipped with
two-fingered Robotiq 2F-140 gripper, and a user interface to start and stop
the experiments.

Due to safety consideration, we trained the proposed throwing policy in the
simulation. After training phase, we conducted experiments in both simulation
and real-robot setups to assess the performance of the learned policy in
throwing various objects into a box located outside the robot’s reachable
area. For evaluation purposes, we designed three tasks with ascending
difficulty levels:

  * •

Task1: obstacle-free object throwing into a static basket randomly placed in
front of the robot. An example of this task in simulation environment is shown
in Fig. 3 (left);

  * •

Task2: obstacle-free object throwing into a moving basket. An example of this
task in real-robot setup is shown in Fig. 3 (right);

  * •

Task3: object throwing while an obstacle obstructs the throwing path as shown
in Fig. 1.

The dimension of the basket was $0.30\times 0.30\times 0.15$ (W $\times$ L
$\times$ H) and obstacle was $0.15\times 0.10\times 0.22$. We used $10$
simulated daily-life objects with different materials, shapes, sizes, and
weight, and $5$ real objects. In particular, five simulated objects were used
during training (i.e., milk_box, coke can, banana, bottle, apple) and the
other five simulated objects were used for testing (i.e., beer_can, peach,
soap, pringles, mustard_bottle). For the real robot experiments, we used five
household objects that are distinct in size and shape from the simulated
objects used during the training phase (i.e., ugly_toy, hello_kitty,
small_box, juice_box, hand_soap).

Figure 4: Four consecutive snapshots showing our simulated robot successfully
thrown a milk box into the basket (Task1).In this round of experiments, the
basket is not reachable by the robot, and the robot should grasp the target
object first, and then, based on the pose of the basket, infers proper
parameters to throw the object into the basket successfully.

###  IV-B Baseline Methods

We employed two sample efficient state-of-the-art off-policy RL algorithms to
train the robot: Deep Deterministic Policy Gradient (DDPG) [16, 17], and Soft
Actor-Critic (SAC) [18, 19] via stable baseline3 [20]. The architectures of
the neural networks for both SAC and DDPG consist of two hidden layers with
the size of $256$ neron per layers with the ReLU activation functions. The
hyper-parameters of SAC are listed in Table. I, and the hyper-parameters of
DDPG that are not listed in Table I

TABLE I: SAC hyper-parameters

Parameter | Value
---|---
#hidden layers (all networks) | $2$
#hidden units per layer | $256$
#samples per minibatch | $256$
optimizer | Adam
learning rate | $3\times 10^{-4}$
batch size | 256
#epochs | 50K
discount ($\gamma$) | $0.99$
replay buffer size | $10^{6}$
nonlinearity | ReLU
target update rate ($\tau$) | 0.005
target update interval | 1
gradient steps | 1

are reported in the Experiments section of our previous work [21]. We also
considered Behavior Cloning (BC), which directly learns a policy (i.e.,
similar to the actor network) by using supervised learning on observation-
action pairs from $25$k successful trials, collected during train and test
phases.

###  IV-C Results

For each of the proposed tasks, we trained the model for $50,000$ steps in
simulation using the five training objects. In the case of simulation, for
each object the robot’s throwing performance was tested for $100$ steps using
the learned policy twice: once with the five test (unseen) objects and once
with the five train (seen) objects. We also used the same learned policy for
real robot experiments and test each of the test objects for $20$ times. As
opposed to “unseen objects”, which is a mixed set of objects not seen during
training, “seen objects” is a mixed set of objects that were used during
training. The average throwing success rates for each approach is reported in
Table II.

Figure 5: A sequence of snapshots showing our simulated robot successfully
thrown a milk box into the moving basket (Task2): In this round of
experiments, the basket is not reachable by the robot, and a human user moves
the basket using an aluminum stick. The robot should first estimate the
direction and velocity of the basket based on visual information, and then
infer proper parameters to throw the milk box into the basket successfully.
Figure 6: Tossing object into the basket while avoiding the obstacle: the
robot should first detect the pose of the object, basket and obstacle, and
then infers appropriate parameters to throw the object into the basket without
colliding with the obstacle.

In the case of Task 1 (i.e., obstacle-free object throwing into a randomly
placed static basket), we observed that for the seen objects the robot with
SAC policy could throw the objects into the basket successfully with $94\%$
accuracy, while with DDPG and BC (on average) achieved $91\%$ and $86\%$
success rate, respectively. A sequence of snapshots demonstrating our robot
successfully throwing a milk box into a basket is depicted in Fig. 4.
Regarding unseen objects, the robot, on average, obtained up to $91\%$ ($91$
success out of $100$ attempts) throwing success rate in simulation and $90\%$
($18/20$) in real-robot setting. In particular, the robot with SAC method, on
average, showed the best throwing performance in real and simulation
experiments, and outperformed DDPG and BC methods with a large margin. The
robot with DDPG achieved the second-best throwing success rate (real: $85\%$
($17/20$), simulation: $81\%$), whereas using BC approach it could get $70\%$
($14/20$) and $81\%$ accuracy in real and simulation experiments respectively.
As expected, the success rate of throwing unseen objects is moderately lower
for all policies. These results showed that the learned policy performs well
both in simulation and for the real robot.

TABLE II: Object Throwing Performance (Mean %). The S/U-S/R shows the setup
configuration: the first token refers to the Seen or Unseen objects, and the
second one denotes Simulation or Real experiments.

Methods | Task 1 | Task 2 | Task 3
---|---|---|---
S-S | U-S | U-R | S-S | U-S | U-R | S-S | U-S | U-R
SAC | 94 | 91 | 90 | 91 | 86 | 85 | 86 | 83 | 80
DDPG | 92 | 81 | 85 | 89 | 86 | 75 | 85 | 77 | 65
BC | 86 | 81 | 70 | 79 | 73 | 65 | 72 | 67 | 55

In the second round of experiments (Task2) the robot should learn to throw a
target object into a moving basket. In the case of simulation experiments, we
randomly selected the direction and a linear speed for moving the basket while
in real-robot experiments, a human user moved the basket using a stick (see
Fig. 5). Similar to the previous round of experiments, the robot with the SAC
policy obtained the best results for both seen and unseen objects. More
specifically, for the seen objects, the robot obtained $91\%$ with the SAC
strategy, whereas its throwing performance with DDPG and BC dropped to $89\%$
and $79\%$, respectively. When tossing unseen objects in simulation, the robot
was $86\%$ accurate with SAC and DDPG policies, compared to $73\%$ accuracy
with BC. Intriguingly, in contrast to simulation results, the robot with DDPG
policy does not perform as competitively to SAC in the real-world.

Task 3 is much more complex than the previous tasks as the robot should learn
to infer an obstacle-free path to throw the object into the basket. An example
of such expriments is shown in Fig. 6. Similar to the previous round of
experiments, the robot with SAC policy outperformed DDPG and BC for both seen
and unseen objects. By comparing all experiments, we observed that the
difference between SAC policy and others is larger when the task is more
difficult. This becomes more visible in the case of unseen objects in real-
robot experiments. In particular, the robot with SAC policy could better model
the solution space and handle unmodeled physical parameters presented in real-
world objects (e.g., aerodynamics, materials, or stiffness). Therefore, in
simulation experiments, the robot using SAC policy obtained $86\%$ and $83\%$
throwing accuracy for seen and unseen objects respectively, while in real-
robot experiments it could achieve $80\%$ success rate. While the performance
of the robot with DDPG policy was marginally lower than SAC policy for seen
objects (i.e., 86% vs. 85%), the difference increased significantly for unseen
objects (i.e., simulation 83% vs. 77%, and in real 80% vs 65%). Experimental
results showed that our formulation maintains the flexibility needed to
express complicated dynamics system while also making learning easier through
trial and error.

###  IV-D Failure Cases

The main reason for failure for all approaches in the simulation experiments
was the inaccurate parameters prediction and throwing of the object near the
basket. In the case of real robot experiments, apart from inaccurate
parameters prediction, we observed three other type of failures: (i)
Inaccurate Tracking was one of the main reasons, where the user move the
object really fast and the tracking could not follow the pose of the object
immediately; (ii) the second primary reason was the lag in executing the
gripper commands on time. In particular, the gripper is controlled through the
robot’s controller and we can not control it directly, i.e., we need to send
the command to the robot’s controller first, and then the controller sends the
command to the gripper, which depends on the status of network and the robot.
(iii) Selecting an unstable grasp pose was the third reason for failure.
Furthermore, we found that for objects with heterogeneous shapes like
ugly_toy, hand_soap, and hello_kitty, the trajectory of the thrown object
varied depending on the grasp pose. In contrast, for homogeneous objects, the
trajectory of the thrown object was not dependent on the grasp pose.

##  V CONCLUSIONS

In this paper, we trained a policy to modulate throwing kernel parameters
through RL to enable robots to precisely throw objects into a moving basket
while there are obstacles obstructing the path. In particular, our method
learns to handle a wide range of situations and avoid colliding the thrown
object with the obstacle. With our formulation, the robot could iteratively
learn the aspects of dynamics that are difficult to model analytically. Due to
safely constraints, we trained the throwing policy in simulation and directly
applied the learned policy in real-robot. We performed extensive sets of
experiments in both simulation and real robot setups in three scenarios with
ascending level of difficulties including, tossing objects into a (i) static
basket, (ii) moving basket, and (iii) obstacle avoidance. Experimental results
showed that the proposed approach enables robot to precisely throw the object
into the basket. We also observed that the learnt policy could be directly
applied to the real robot even though it had only been trained in simulation.
It also showed outstanding generalization capabilities to new target locations
and unknown objects. In the continuation of this work, we would like to
investigate the possibility of enhancing throwing performance by taking into
account additional sensory modalities (e.g., force or tactile sensors), as
these could help the robot grasp objects steadily and adjust its throwing
parameters accordingly.

## References

  * [1] J. E. Kuhn, “Throwing, the shoulder, and human evolution.” _American Journal of Orthopedics (Belle Mead, NJ)_ , vol. 45, no. 3, pp. 110–114, 2016.
  * [2] J.-S. Hu, M.-C. Chien, Y.-J. Chang, S.-H. Su, and C.-Y. Kai, “A ball-throwing robot with visual feedback,” in _2010 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2010, pp. 2511–2512.
  * [3] Y. Gai, Y. Kobayashi, Y. Hoshino, and T. Emaru, “Motion control of a ball throwing robot with a flexible robotic arm,” _International Journal of Computer and Information Engineering_ , vol. 7, no. 7, pp. 937–945, 2013.
  * [4] A. Zeng, S. Song, J. Lee, A. Rodriguez, and T. Funkhouser, “Tossingbot: Learning to throw arbitrary objects with residual physics,” _IEEE Transactions on Robotics_ , vol. 36, no. 4, pp. 1307–1319, 2020.
  * [5] A. Takahashi, M. Sato, and A. Namiki, “Dynamic compensation in throwing motion with high-speed robot hand-arm,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 6287–6292.
  * [6] K. Ploeger and J. Peters, “Controlling the cascade: Kinematic planning for n-ball toss juggling,” _arXiv preprint arXiv:2207.01414_ , 2022.
  * [7] M. T. Mason and K. M. Lynch, “Dynamic manipulation,” in _Proceedings of 1993 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS’93)_ , vol. 1. IEEE, 1993, pp. 152–159.
  * [8] J. Kober, K. Muelling, and J. Peters, “Learning throwing and catching skills,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 5167–5168.
  * [9] D. M. Lofaro, R. Ellenberg, P. Oh, and J.-H. Oh, “Humanoid throwing: Design of collision-free trajectories with sparse reachable maps,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 1519–1524.
  * [10] A. Ghadirzadeh, A. Maki, D. Kragic, and M. Björkman, “Deep predictive policy training using reinforcement learning. in 2017 ieee,” in _RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 2351–2358.
  * [11] J. Kober, E. Oztop, and J. Peters, “Reinforcement learning to adjust robot movements to new situations,” in _Twenty-Second International Joint Conference on Artificial Intelligence_ , 2011.
  * [12] S. H. Kasaei, M. Oliveira, G. H. Lim, L. S. Lopes, and A. M. Tomé, “Towards lifelong assistive robotics: A tight coupling between object perception and manipulation,” _Neurocomputing_ , vol. 291, pp. 151–166, 2018\.
  * [13] S. H. Kasaei, N. Shafii, L. S. Lopes, and A. M. Tomé, “Interactive open-ended object, affordance and grasp learning for robotic manipulation,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 3747–3753.
  * [14] H. Kasaei, S. Luo, R. Sasso, and M. Kasaei, “Simultaneous multi-view object recognition and grasping in open-ended domains,” _arXiv preprint arXiv:2106.01866_ , 2021.
  * [15] S. Kasaei, J. Sock, L. S. Lopes, A. M. Tomé, and T.-K. Kim, “Perceiving, learning, and recognizing 3d objects: An approach to cognitive service robots,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, 2018.
  * [16] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.
  * [17] D. Silver, G. Lever, N. Heess, T. Degris, D. Wierstra, and M. Riedmiller, “Deterministic policy gradient algorithms,” in _International conference on machine learning_. PMLR, 2014, pp. 387–395.
  * [18] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _International conference on machine learning_. PMLR, 2018, pp. 1861–1870.
  * [19] T. Haarnoja, A. Zhou, K. Hartikainen, G. Tucker, S. Ha, J. Tan, V. Kumar, H. Zhu, A. Gupta, P. Abbeel _et al._ , “Soft actor-critic algorithms and applications,” _arXiv preprint arXiv:1812.05905_ , 2018.
  * [20] A. Raffin, A. Hill, A. Gleave, A. Kanervisto, M. Ernestus, and N. Dormann, “Stable-baselines3: Reliable reinforcement learning implementations,” _Journal of Machine Learning Research_ , vol. 22, no. 268, pp. 1–8, 2021\. [Online]. Available: http://jmlr.org/papers/v22/20-1364.html
  * [21] S. Luo, H. Kasaei, and L. Schomaker, “Accelerating reinforcement learning for reaching using continuous curriculum learning,” in _2020 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2020, pp. 1–8.


11institutetext: EPITA Research and Development Lab. (LRDE), EPITA, France
22institutetext: Univ. Gustave Eiffel, IGN-ENSG, LaSTIG 33institutetext:
LaDéHiS, CRH, EHESS

# Combining Deep Learning and Mathematical Morphology for Historical Map
Segmentation

Yizi Chen 1122 Edwin Carlinet 11 Joseph Chazalon 11
Clément Mallet 22 Bertrand Duménieu 33 Julien Perret 2233

###### Abstract

The digitization of historical maps enables the study of ancient, fragile,
unique, and hardly accessible information sources. Main map features can be
retrieved and tracked through the time for subsequent thematic analysis. The
goal of this work is the vectorization step, i.e., the extraction of vector
shapes of the objects of interest from raster images of maps. We are
particularly interested in closed shape detection such as buildings, building
blocks, gardens, rivers, etc. in order to monitor their temporal evolution.
Historical map images present significant pattern recognition challenges. The
extraction of closed shapes by using traditional Mathematical Morphology (MM)
is highly challenging due to the overlapping of multiple map features and
texts. Moreover, state-of-the-art Convolutional Neural Networks (CNN) are
perfectly designed for content image filtering but provide no guarantee about
closed shape detection. Also, the lack of textural and color information of
historical maps makes it hard for CNN to detect shapes that are represented by
only their boundaries. Our contribution is a pipeline that combines the
strengths of CNN (efficient edge detection and filtering) and MM (guaranteed
extraction of closed shapes) in order to achieve such a task. The evaluation
of our approach on a public dataset shows its effectiveness for extracting the
closed boundaries of objects in historical maps.

###### Keywords:

Deep Learning Convolutional Neural Networks Mathematical Morphology Historical
Map Segmentation Object Extraction

##  1 Introduction

(a) Some geographical entities typically depicted in city maps: building
blocks (orange), roads (green) and rivers (blue).

(1)(2)(3) (b) Challenges in historical maps: (1) planimetric overlap, (2) text
overlap, (3) paper folds.

Figure 1: Contents of a 1925 urban topographic map along with an overview of
their challenging properties for automatic feature extraction. ††Extra
material for this paper (full-size figures, results, code, dataset) available
at:
https://github.com/soduco/paper-dgmm2021

The massive digitization of archival collections carried out by heritage
institutions provides access to huge volumes of historical information encoded
in the available documents. Among them, maps are unfortunately still little
exploited. Yet they are a gold mine of geographic data that allows to
reconstruct and analyze the morphological and social evolution of a place over
time [21, 11]. Topographic maps, in particular, engrave many geographical
features: their distribution in space, their topological relationships and
various information encoded by the map legend or by text labels [14, 7].
Transforming such graphical representations of geographic entities into
discrete geographic data (or vector data) is a crucial step for numerous
spatial and spatio-temporal analysis purposes. Such a transformation is most
often manually retrieved by historians or with the help of crowdsourcing
tools. This is extremely time-consuming, non-reproducible, and leads to
heterogeneous data quality. Automating this tedious task is a key step towards
building large volumes of reference geo-historical data.

Unfortunately, historical maps exhibit characteristics that hinder standard
pattern recognition approaches and make them relatively inefficient at
extracting data of good quality, i.e., that do not need to be manually post-
processed. Unlike modern computer-generated maps which follow roughly the same
semiotic rules, these maps vary in terms of legend, level of generalization,
type of geographic features and text fonts [14]. They also usually lack
texture information, which creates ambiguities in the detection of objects.
For instance, building blocks and roads have very similar textures despite
being of completely different nature (Figure 1(a)). Popular semantic [15, 2,
25] and instance [6, 24, 3] image segmentation algorithms detect objects based
on textures and are prone to fail in our context. Color is not a relevant cue
either: the palette is usually highly restricted due to the technical
limitations and financial constraints of their production. Objects in maps are
often overlapping, some are thus partially hidden and hardly separable.
Occlusion happens with overlayed textual and carto-geodetic information in
particular (Figure 1(b), rectangles (1) and (2)). Last, preservation
conditions of historical maps play a role as stains, folds or holes might
cause gaps in the cartographic information. Such artifacts may lead to
incorrect object detection (Figure 1(b), rectangle (3)).

Historical MapDeep ContourDetectionEdge Probability MapFiltering
&WatershedClosed Shapes

Figure 2: Overview of the approach presented in the paper: we combine an
efficient edge detection and filtering stage using a deep network with a fast
closed shape extraction using mathematical morphology tools.

Our contributions in this paper are as follows. After reviewing the
limitations of the current approaches for segmenting maps in section 2, we
propose a simple pipeline (fig. 2) that combines deep networks and
mathematical morphology for object detection in maps. It takes benefit from
their complementary strengths, namely image filtering and strong guarantees
with respect to closed shapes. We derive edge probability maps using a multi-
scale deep network approach depicted in section 3 and then leverage
mathematical morphology tools to extract closed shapes as explained in section
4. Eventually, in section 5, the second contribution lies in a thorough
evaluation of the relevance of the mathematical morphology stage with novel
visualizations and metrics to objectively assess our approach and better
identify the strengths and weaknesses of each stage and of the workflow.

##  2 Approaches for map segmentation

We target to recover geometric structures from scans of historical maps. As
mentioned above, due to the limited texture and colour content of such data
sources, standard semantic segmentation approaches of the literature would
fail for most cases. Instead, we cast our problem as a vectorization challenge
that can be turned into a region-based contour extraction task. Such a problem
is traditionally solved through a two-step approach: the detection of edges or
local primitives (lines, corners) followed by the retrieval of structures
based on global constraints [31]. Recent works have shown the relevance of a
coupled solution [12]. They remain tractable and efficient only for a limited
number of structures. Region-based methods (e.g., based on PDEs [18]) may lead
to oversimplified results and will not be further analysed here.
The main issue of two-step solutions is the edge detection step. This low-
level task is achieved by measuring locally pixel gradients. Due to the amount
of noise (overlapping objects, map deformation), this would result in many
tiny and spurious elements that any global solution would manage connecting.
Instead, we focus on boundary detection, i.e., a middle-level image task that
separates objects at the semantic level according to different geometric
properties of images. This offers two main advantages: (i) a limited
sensitivity to noise in maps and (ii) the provision of more salient and robust
primitives for the subsequent object extraction step. We do not focus on a
primitive-based approach since shapes on maps cannot be simply assumed.

Recently, among the vast amount of literature, convolution neural networks
(CNN) have shown a high level of performance for boundary detection [30, 13].
However, they only provide probability edge maps. Without topological
constraints, image partitioning is not ensured. Conversely, watershed
segmentation techniques in mathematical morphology can directly extract closed
contours. They run fast for such a generation, but may lead to many false-
positive results. Indeed, using only low-level image features such as image
gradients, watershed techniques may not efficiently maintain useful boundary
information [3]. Consequently, we propose here to merge the CNN-based and
watershed image segmentation methods in order to benefit from the strengths of
both strategies [29]. A supervised approach is conceivable since we both have
access to reference vectorized maps and CNN architectures pre-trained with
natural image.

##  3 Deep Edge Detection

We detail how we selected the network architecture used to detect and filter
edges, with illustrations of the strengths of such approach, and describe the
training procedure we followed to use the selected network (BDCN) on our
dataset.

Network architecture Contour detection was first addressed with the design of
handcrafted features based on brightness, color, textures [16]. Then,
improvements lied in their efficient group through mono- or multi-scale
attributes retrieving micro-structures: textons are a salient example [32].
Afterwards, main methods focused on combining all available cues, such as [1].
They used a global probability boundary by learning the weights of manually
selected features (gradients and textons as features in several image scales)
in order to detect contours and form better closed boundaries to represent the
objects in images. Since CNNs have proved their relevance to extract and
combine meaningful image features, a large amount of research has focused on
detecting contours. The most famous one is the so-called Holistically edge
detector (HED) [30], which is an end-to-end multi-scale deep learning network.
The novelty consisted in using skip-connections to merge different levels of
features and learn different losses from intermediate layers of VGG-16 [27].
This allowed to recover multiscale representations of image features.
Eventually, He et al. [13] proposed a so-called Bi-Directional Cascade Network
(BDCN) by designing a scale enhancement module (SEM) on top of HED to enhance
multiscale spatial contexts in images resulting in a better performance than
humans in the BSDS500 dataset.

One advantage of BDCN is that the multiscale representatives combine
semantically meaningful features to efficiently filter out the image textures
and text information while maintaining useful contours and lines in the
images. It is particularly suited for handling noise in our maps. Another
advantage is that learnable dilated convolutions in SEM can learn fine-grained
features with larger receptive fields that are beneficial when we want to
separate the texts with object contours accurately. It is because building
contours have much longer pixel continuity than text, resulting in higher
activation after dilated convolution. After several iterations, the
probability of text pixels will vanish, leading to their removal, similarly to
texture, as shown in fig. 3. However, the BDCN network works only at the pixel
level and cannot guarantee the required topological properties in predicted
edge probability maps without additional topological constraints [8], thus the
current solution requires knowledge of the number of structures to be
retrieved.

Figure 3: BDCN produces an Edge Probability Map (_right_) with texts and
textures removed from the input (_left_).

Training Annotated historical maps are used to train a BDCN network. The final
prediction which is a probability map where each pixel in the maps contain
values in range $[0,1]$ (zero means the pixel does not belong to a contour,
one that it does). We train our network from scratch instead of using transfer
learning on the edge weights learned from BSDS500 (dataset developed for image
boundary detection and segmentation tasks): the features in natural images are
very different from our historical map images.We need to filter out most of
the texts in our maps, but the network trained on the BSDS dataset does not
provide any useful features related to geometric filtering tasks. In order to
handle data imbalance during training, we proceed as follows. We define our
input image as $x\in$ $\mathbb{R}^{H\cdot W}$ and ground truth label $y\in$
${0,1}^{H\cdot W}$. The output of predicted image is
$\hat{y}=f(x,w)\in{0,1}^{H\cdot W}$ and every element of $\hat{y}$ is
interpreted as the probability of pixel i having label 1: $\hat{y}\equiv
p(Y_{i}=1|x,w)$. Since the edge detection is a binary classification task,
binary cross entropy loss is used as loss function between predictions and
ground truths. Due to highly imbalanced edge (97.5%) and non-edge (only 2.5%)
classes, extra parameters $\alpha,\beta$ are used as weights to re-balance the
binary cross entropy loss, as $\mathcal{L}_{BCE}=-\alpha\sum_{j\in
Y_{-}}log(1-\hat{y_{j}})-\beta\sum_{j\in Y_{+}}log(\hat{y_{j}})$ where $Y_{+}$
is the set of indices of edge pixels, $Y_{-}$ is the set of indices of non-
edge pixels, $\alpha=\left(\lambda\cdot|Y_{-}|/(|Y_{+}|+|Y_{-}|)\right)$ is
the percentage of edge pixels in each batch of historical map image and
$\beta=\left(|Y_{+}|/(|Y_{+}|+|Y_{-}|))\right)$ is percentage of non-edge
pixels. An extra $\lambda=1.1$ factor is used to enhance the percentage of
edge pixels in order to give extra weights for edge responses.
We build our code based on the BDCN code repository to train our historical
map dataset from scratch with a few modifications. We evaluate the loss for
every epoch and also for choosing the best training weights. To make the
network converging faster, we replace SGD with ADAM optimizer. The initial
learning rate is set to $5\times 10^{-5}$ with $0.9$ momentum and $0.002$
weight decay.

##  4 Segmentation of the EPM

From the Edge Probability Map, we then need to extract boundaries of the
objects. In Mathematical Morphology, the Watershed Transform [17] is a de-
facto standard approach for image segmentation. It has been used in many
applications and has been widely studied in terms topological properties [10,
23], in terms of algorithms and in terms on computation speed [23, 9].

It has two well-known issues: the over-segmentation due to the high number of
minima, and the gradient leakage that merges regions. There is a third general
issue with the watershed that concerns the separation of overlapping or
touching objects, but this is not a problem in our case since the map
components do not overlap.

Solutions to the over-segmentation problem The first problem is generally
solver by filtering the minima first. In [28], the _h_ -minima characterize
the importance of each local minimum through their _dynamic_. When flooding a
basin, it actually refers to the water elevation required to merge with
another basin. Attributes filter, filters by reconstruction [26] also allow to
eliminate some minima based on their algebraic properties: size, shape,
volume… Another efficient approach consists in first ordering the way the
basins merge to create a hierarchy of partitions and then performing a cut in
the hierarchy to get a segmentation with non-meaningful basins removed [5, 20,
4].

Solutions to the early leakage problem The second problem lies in the quality
of the gradient. It has been noted [19], that (hierarchical) watersheds have
better results on non-local supervised gradient estimators. The idea of
combining the watershed with high performance contour detector dates back to
[1].

The relevance of a simple closing by area and dynamic on the edge map produced
by our deep-learning edge detector combined with the watershed for this
application lies in three points.

First, the minimum size of the components is known. Indeed, the document
represents a physical size, and regions whose area is below 100 $m^{2}$ are
not represented in the map. Thus, we have a strong _a priori_ knowledge we
want to inject in the process, the minimum size of the regions (in pixels).
This type of constrain is hard to infer in a deep-learning system and we
cannot have such guarantees from its output. Having hard guaranties about the
shapes and their size is at the foundation of the granulometries in
Mathematical Morphology. Moreover, the connected (area) filter used for
filtering the edge image ensure that we do not distort the signal at the
boundaries of the meaningful regions.

Second, the watershed segmentation method does not rely on the strength of the
gradients to select the regions. Even if the edge response is low (i.e, the
gradient is weak), the watershed is able to consider this weak response and
closes the contour of the region. We do not depend on the strength of the edge
response from BDCN which is difficult to calibrate and normalize.

Last but not least, not only the watershed outputs a segmentation but some
implementations also produce watershed lines between regions. In our
application, watershed lines are even more important than regions because we
need to extract polygons for each shape. Event if we could extract boundaries
from regions, it avoids an extra processing step. The watershed lines produced
by the algorithm is one pixel-large and are located where the edges are the
strongest, i.e., where the network has the strongest response on thick edges.
The watershed lines form closed boundaries around regions which is a guarantee
we cannot have from the output of a network.

Input EPM Ground-truth Param set A Param set B

Figure 4: Some _failures_ and some _success stories_ of the watershed
segmentation. The parameter sets are A: _h=3_ , $\lambda=250$, and B: _h=7_ ,
$\lambda=400$. The first row shows the ability to recover weak boundaries.
This sensitivity is not desirable in some cases as it leads the over-
segmentation of the 2nd row. The third row suggests that the over-segmentation
can be prevented by a stronger filtering but would also lead to a lower shape
detection.

Figure 4 shows the strength of the watershed to recover the boundaries of
objects even on weak edge responses that would be lost by thresholding the
EPM. This is especially visible in the first row where the boundaries of
“Place du Châtelet” are leaking; nevertheless they are recovered in the
segmentation. On the downside, this ability to recover weak edges is also a
bottleneck that can create false-boundaries as shown in the middle row where
the place around “Eglise Notre-Dame” is over-segmented because of some
detection noise.

The filtering parameters (dynamic $h$ and area $\lambda$) are important to
control the trade-off between the fact we want to recover small/leaking
regions (somewhat related to the _recall_) and the false-detection of
boundaries (somewhat related to the _precision_). This is illustrated with two
sets of parameters A and B where B has more restrictive filtering parameters.
The third row of Figure 4 shows that B has less over-segmentation but in the
two first rows, it misses some boundaries.

The decision to merge objects is actually very context-dependant, as it does
not depend on the size of the component, neither its volume, nor its shape.
The watershed “does its best” to create the missing boundaries and, at the
moment, we have not managed to find better rules (e.g., with extinction values
of some attributes) to filter out the basins of the watershed.

##  5 Evaluation

To assess the performance of the proposed approach, we conducted a series of
experiments on a fully manually annotated map sheet. We report here details
about this dataset we created and used, the experimental protocol as well as
the calibration procedures we followed, the metrics we designed and used, and
discuss some results.

Dataset Among the multiple map sheets of the collection of Paris atlases, our
work focuses on the particular sheet representing a central area of the city
from year 1925 [22]. We encourage the reader to refer to the supplementary
material of this paper for a full-size view of this image. Indeed, such map
sheets are large by nature and were digitized with high resolution, resulting
in a 8500$\times$6500 image for the area of interest.
We carefully annotated the original image by creating line vector information
for each edge of each object of interest in the map. It should be noted that
only a subset map strokes should be kept as many objets are not relevant for
our current study: underground lines and railways, for instance, should be
discarded. The resulting vector information was rastered to produce: i) a
reference edge map (a small dilation was applied so the resulting edges have a
thickness of 3 pixels); ii) a reference label map identifying each shape to be
detected.
We divided the image into three disjoint subsets: a training set (rows 0 to
3999); a validation set (rows 4000 to 4999); and a test set (rows 5000 to
6500). These areas were divided into 228 disjoint tiles of 500$\times$500
pixels.

Protocol Our goal in the evaluation protocol we designed was the assess the
impact of the watershed stage in our pipeline. We compared compared the
performance of a baseline system, without watershed, with our proposed
approach: the same baseline augmented by a watershed stage (see fig. 2).
The baseline (without watershed) consists in a deep edge detection stage using
the BDCN network presented in section 3. This stage produces an edge
probability map (EPM) as previously explained. The network was trained on the
training set using the validation set as control set during training. To
generate closed shapes, we simply thresholded the EPM and extracted the
connected components. We selected the best performing threshold value (9) on
the validation set for fair comparison.
The proposed approach (baseline plus watershed) consists in adding a joint
filtering on area and dynamic of the EPM followed by a watershed. This
approach produces a label map, i.e. a usable set of closed shapes, as detailed
in section 4. We selected the best performing values for area ($\lambda$) and
dynamic $h$ parameters on the validation set.
To avoid losing topological information during component labeling (baseline)
or during watershed, theses steps were performed on the full image (with
training, validation and test sets merged) but the performance indicators were
computed exclusively on the test set by masking other areas.

| CC-labeling | Watersheding
---|---|---
IoU | Precision | Recall | F-score | Precision | Recall | F-score
0.50 | 0.20 | 0.39 | 0.27 | 0.74 | 0.50 | 0.59
0.80 | 0.10 | 0.19 | 0.13 | 0.60 | 0.40 | 0.48
0.90 | 0.04 | 0.09 | 0.06 | 0.45 | 0.30 | 0.36
0.95 | 0.01 | 0.02 | 0.02 | 0.25 | 0.16 | 0.20

Figure 5: _Left_ : comparison of the evolution of the shape detection F1-score
across all possible IoU threshold with and without the watershed stage.
_Right_ : evaluation metrics with and without watershed.

Metrics While it is common in segmentation challenges to evaluate the quality
of object detection by evaluating the precision and recall of edge detection
at pixel, such an approach would only evaluate the process halfway to our
target application: closed shapes detection. To evaluate shape detection, we
need to identify pairs of matching shapes between a reference set ($R$) and a
predictions set ($P$). Because, in our particular case, shapes are disjoint
among $R$ and also among $P$ (by construction), we can leverage the following
property: as soon as the intersection over union (_IoU_) between $r_{i}\in R$
and $p_{j}\in P$ is strictly superior to $0.5$, then we know that no other
element $r_{k}\in R,i\neq k$ can have a higher IoU with $p_{j}\in R$ than
$r_{i}\in R$, and reciprocally.
As soon as a pair of shapes $(r_{i},p_{j})\in R\times P$ exhibits
$\mathrm{IoU}(r_{i},p_{j})=\mathrm{area}(\frac{r_{i}\cap p_{j}}{r_{i}\cup
p_{j}})=T>0.5$ then we count a successful match under the threshold constraint
$T$. We introduce this threshold value so as to consider all possible values
between 0.5 (exclusive) and 1 (inclusive) and create a global indicator of the
system under all potential quality requirements. This allows us to count the
number of correctly detected shapes (_true positives_ or $\mathrm{TP}$),
missed shapes ((_false negatives_ or $\mathrm{FN}$)), and wrongly predicted
shapes ((_false positives_ or $\mathrm{FP}$)) for every operating
characteristics. We derive from this set of measures two analysis tools.
First a precision ($\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$), a recall
($\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$) and a F1 score
($\frac{2\mathrm{TP}}{2\mathrm{TP}+\mathrm{FN}+\mathrm{FP}}$) curves for all
possible threshold values. They offer a condensed view of the behavior of a
system under all possible operating characteristics. We suggest to use the
area under the F1 score curve to compute a single performance indicator. The
values of this indicator ranges from 0 (worst) to 0.5 (best).
The second tool is a pair of visualization maps: a precision map which
associates for each predicted shape $p_{j}\in P$ the maximal IoU value
$b_{pj}$ such as $b_{pj}=\mathrm{argmax}_{r_{i}\in
R}(\mathrm{IoU}(r_{i},p_{j}))$, and a recall map which associates for each
expected shape $r_{i}\in R$ the maximal IoU value $b_{ri}$ such as
$b_{ri}=\mathrm{argmax}_{p_{j}\in P}(\mathrm{IoU}(r_{i},p_{j}))$. Each pixel
of each shape is then assigned a color indicating the value of the maximal
IoU: red to yellow for values between 0 and 0.5, and yellow to green for
values between 0.5 and 1. The darker the green the better the match in both
maps. The darker the red, the more serious the false positive (resp. negative)
in precision (resp. recall) map.

Results and Discussion We report here the results for the best calibrated
variant of each of the two systems (baseline+connected component labeling vs
baseline+watershed) under test. Figure 5 (left) compares the evolution of the
F1 score indicator for both systems under each possible IoU threshold. Figure
5 (right) details the different indicators for several key values of IoU
thresholds. We can see from those results that the watershed post-processing
consistently and significantly improves the quality of the results. The
precision and recall maps presented in fig. 6 illustrate the benefits that the
watershed post-processing bring to the deep edge segmentation: it adjusts the
border of the shapes (improves precision and recall); it also removes small
noise (improves precision); and it also efficiently recovers some weak
boundaries (improves recall).

|  |
---|---|---
Precision map before watershed |  | Precision map after watershed
|  |
Recall map before watershed |  | Recall map after watershed
Figure 6: Precision and recall maps without and with watershed.

##  6 Conclusion

In this paper, we propose an efficient combination of convolutional neural
networks and mathematical morphology to address the problem of closed shapes
extraction in historical maps. Convolutional neural networks (BDCN) allow us
to efficiently detect edges while filtering unwanted features (text for
instance). Mathematical morphology is applied to the edge probability map
created by BDCN to create closed shapes reliably. The efficiency of our
approach is shown by testing it on an open dataset. We believe such a method
will make the digitization process of historical maps faster and more
reliable.

Acknowledgements
This work was partially funded by the French National Research Agency (ANR):
Project SoDuCo, grant ANR-18-CE38-0013.

## References

  * [1] Arbelaez, P., Maire, M., Fowlkes, C., Malik, J.: Contour detection and hierarchical image segmentation. IEEE Trans. on Pattern Analysis and Machine Intelligence 33(5), 898–916 (2010)
  * [2] Badrinarayanan, V., Kendall, A., Cipolla, R.: Segnet: A deep convolutional encoder-decoder architecture for image segmentation. IEEE Trans. on Pattern Analysis and Machine Intelligence 39(12), 2481–2495 (2017)
  * [3] Bai, M., Urtasun, R.: Deep watershed transform for instance segmentation. In: Proc. of Conf. on Computer Vision and Pattern Recognition. pp. 5221–5229 (2017)
  * [4] Barcelos, I.B., da Fonseca, G.B., Najman, L., Kenmochi, Y., Perret, B., Cousty, J., do Patrocínio, Z.K., Guimarães, S.J.F.: Exploring hierarchy simplification for non-significant region removal. In: SIBGRAPI Conference on Graphics, Patterns and Images. pp. 100–107 (2019)
  * [5] Beucher, S.: Watershed, hierarchical segmentation and waterfall algorithm. In: Proc. of Intl. Symp. on Mathematical Morphology (ISMM), pp. 69–76 (1994)
  * [6] Chen, K., Pang, J., Wang, J., Xiong, Y., Li, X., Sun, S., Feng, W., Liu, Z., Shi, J., Ouyang, W., et al.: Hybrid task cascade for instance segmentation. In: Proc. of Conf. on Computer Vision and Pattern Recognition. pp. 4974–4983 (2019)
  * [7] Chiang, Y.Y., Leyk, S., Knoblock, C.A.: Efficient and robust graphics recognition from historical maps. In: Intl. Workshop on Graphics Recognition. pp. 25–35. Springer (2011)
  * [8] Clough, J.R., Oksuz, I., Byrne, N., Schnabel, J.A., King, A.P.: Explicit topological priors for deep-learning based image segmentation using persistent homology. In: Intl. Conf. on Information Processing in Medical Imaging. pp. 16–28. Springer (2019)
  * [9] Couprie, M., Najman, L., Bertrand, G.: Quasi-linear algorithms for the topological watershed. J. of Mathematical Imaging and Vison 22(2-3), 231–249 (2005)
  * [10] Cousty, J., Bertrand, G., Najman, L., Couprie, M.: Watershed cuts: Thinnings, shortest path forests, and topological watersheds. IEEE Trans. on Pattern Analysis and Machine Intelligence 32(5), 925–939 (2009)
  * [11] Dietzel, C., Herold, M., Hemphill, J.J., Clarke, K.C.: Spatio-temporal dynamics in california’s central valley: Empirical links to urban theory. Intl. J. of Geographical Information Science (IJGIS) 19(2), 175–195 (2005)
  * [12] Favreau, J., Lafarge, F., Bousseau, A., Auvolat, A.: Extracting geometric structures in images with delaunay point processes. IEEE Trans. on Pattern Analysis and Machine Intelligence 42(4), 837–850 (2020)
  * [13] He, J., Zhang, S., Yang, M., Shan, Y., Huang, T.: Bdcn: Bi-directional cascade network for perceptual edge detection. IEEE Trans. on Pattern Analysis and Machine Intelligence (2020)
  * [14] Leyk, S., Boesch, R., Weibel, R.: Saliency and semantic processing: Extracting forest cover from historical topographic maps. Pattern recognition 39(5), 953–968 (2006)
  * [15] Long, J., Shelhamer, E., Darrell, T.: Fully convolutional networks for semantic segmentation. In: Proc. of Conf. on Computer Vision and Pattern Recognition. pp. 3431–3440 (2015)
  * [16] Martin, D., Fowlkes, C., Tal, D., Malik, J.: A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. In: Proc. of Intl. Conf. of Computer Vision (ICCV). vol. 2, pp. 416–423 (2001)
  * [17] Meyer, F.: Topographic distance and watershed lines. Signal processing 38(1), 113–125 (1994)
  * [18] Orzan, A., Bousseau, A., Winnemöller, H., Barla, P., Thollot, J., Salesin, D.: Diffusion curves: A vector representation for smooth-shaded images. ACM Transactions on Graphics 27(3), 1–8 (2008)
  * [19] Perret, B., Cousty, J., Guimaraes, S.J.F., Maia, D.S.: Evaluation of hierarchical watersheds. IEEE Trans. on Image Processing 27(4), 1676–1688 (2017)
  * [20] Perret, B., Cousty, J., Guimarães, S.J.F., Kenmochi, Y., Najman, L.: Removing non-significant regions in hierarchical clustering and segmentation. Pattern Recognition Letters 128, 433–439 (2019)
  * [21] Perret, J., Gribaudi, M., Barthelemy, M.: Roads and cities of 18th century france. Scientific data 2(1), 1–7 (2015)
  * [22] Préfecture de la Seine, service du Plan: Atlas des vingt arrondissements de Paris [1 vol. (3 pl., 16 pl. doubles), 68 cm]. Paris. L. Wuhrer. ARK: 73873/pf0000935524 (1925), Bibliothèque de l’Hôtel de Ville, Paris.
  * [23] Roerdink, J.B., Meijster, A.: The watershed transform: Definitions, algorithms and parallelization strategies. Fundamenta informaticae 41(1, 2), 187–228 (2000)
  * [24] Romera-Paredes, B., Torr, P.H.S.: Recurrent instance segmentation. In: Proc. of Eur. Conf. on Computer Vision (ECCV). pp. 312–329. Springer (2016)
  * [25] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Proc. of Medical Image Computing and Computer Assisted Intervention (MICCAI). pp. 234–241 (2015)
  * [26] Salembier, P., Serra, J.: Flat zones filtering, connected operators, and filters by reconstruction. IEEE Trans. on Image Processing 4(8), 1153–1160 (1995)
  * [27] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition (2015)
  * [28] Soille, P.: Morphological image analysis: principles and applications. Springer Science & Business Media (2013)
  * [29] Xie, L., Qi, J., Pan, L., Wali, S.: Integrating deep convolutional neural networks with marker-controlled watershed for overlapping nuclei segmentation in histopathology images. Neurocomputing 376, 166–179 (2020)
  * [30] Xie, S., Tu, Z.: Holistically-nested edge detection. In: Proc. of Conf. on Computer Vision and Pattern Recognition. pp. 1395–1403 (2015)
  * [31] Zhang, Z., Fidler, S., Waggoner, J., Cao, Y., Dickinson, S., Siskind, J., Wang, S.: Superedge grouping for object localization by combining appearance and shape informations. In: Proc. of Conf. on Computer Vision and Pattern Recognition. pp. 3266–3273 (2012)
  * [32] Zhu, S.C., Guo, C.E., Wang, Y., Xu, Z.: What are textons? Intl. J. of Computer Vision 62(1-2), 121–143 (2005)


##  1 METHODOLOGY

While the sub-arcsecond angular separation of two images of a GRB lensed by
the mass stated before will not be resolved by Fermi/LAT and the GRB will be
recorded as one event, the light curve recorded by Fermi/GBM will be a
superposition of the two light curves corresponding to the two lensed images
and this can be resolved when the time delay is less than the GRB recording
time and larger than burst duration. [Ji2018StrongLO] used the autocorrelation
function to probe this scenario and extract the time delay in the order of 1
ms in the candidate GRB light curves due to lenses with the mass range of
$10^{2}M_{\odot}\leq M_{L}\leq 10^{3}M_{\odot}$ This approach is optimized for
Fermi/GBM so it is sensitive to lensed GRBs with time delays longer than burst
duration and shorter than 300 seconds. Here the typical time delay is in order
of 10 seconds which corresponds to $M_{L}=10^{6}M_{\odot}$.

###  Simulating Lensing Data

We simulated a lensed light curve by superimposing the count rate $I(t)$ of a
real GRB with itself as shown in Fig.1, by imposing a time delay $\Delta t$
and a magnification ratio $R$. The resulting lensed light curve count rate is

$I^{\prime}(t)=\frac{R}{R+1}I(t)+\frac{1}{R+1}I(t+\Delta t).$ (1)

The crucial feature of this equation is that the fainter peak of the light
curve, which corresponds to the image with the lower magnification, will
arrive later than the brighter one, which corresponds to the image with the
higher magnification (i.e. $t^{+}\leq t^{-}$) as derived from Eqs.
(LABEL:eq:delay) and (LABEL:eq:magnification).

For a specific lensing configuration with normalized impact parameter
$\beta/\theta_{E}=0.7$, lens mass of $5\times 10^{5}M_{\odot}$ and
$D_{s}=2D_{L}$ (mentioned in section LABEL:sec:GLensing), and using equations
(LABEL:eq:delay) and (LABEL:eq:magnification), we find the corresponding time
delay and magnification for the simulated GRB in Fig.1 to be $\Delta t=28s$
and $R=4$.

[width=0.7]Methodology/LCoriginal_mock.png

Figure :  Top: The light curve of an original GRB (GRB 200714774). Bottom:
Simulated lensed light curve of the same GRB due to a point lens with a mass
of $5\times 10^{5}M_{\odot}$ at a cosmological redshift of half that of GRB
($z_{L}=1$ distance $\sim 1Gpc$), and with a normalized impact parameter
$\beta/\theta_{E}=0.7$. This assumed lens makes two lensed images with a time
delay of 28 seconds and a magnification ratio of 4.

###  Autocorrelation Function Technique

We apply the autocorrelation function as a technique to identify two peaks in
the GRB light curve due to the two lensed images and to obtain the time delay.
The autocorrelation function of count rate $I(t)$ with a time lag of $\tau$ is
defined by,

$C(\tau)=\sum_{t}^{n(\tau)=T-\tau}I(t)I(t+\tau).$ (2)

Applying the autocorrelation function on a simulated light curve that contains
two local peaks due to two lensed images with an artificial time delay $\Delta
t$, will yield a maximum or local spike in $C(\tau)$ at a time lag $\tau$ that
corresponds to the delay time $\Delta t$. For the simulated light curve in
Fig.1(bottom panel), a spike is identified at $\tau=28s$ which represents the
time delay between the lensed images as shown in Fig.1. While [Ji2018StrongLO]
used a normalized autocorrelation function to determine the time delay between
the gravitationally lensed images, which is less sensitive or yields
suppressed spikes, we use a non-normalized autocorrelation function which
gives better results and higher sensitivity. We note in Appendix
LABEL:section:appendixA that the error in estimation of the unnormalized
($C(\tau)$) and the normalized ($C_{N}(\tau)$) autocorrelation function are
related together as $\Delta
C_{N}(\tau)\geq\sqrt{C(\tau)^{2}+(\frac{\Delta\Sigma_{1}}{\Sigma_{1}})^{2}+(\frac{\Delta\Sigma_{2}}{\Sigma_{2}})^{2}}$,
where $\Sigma_{1}=\sum I(t)^{2}\Delta t$ and $\Sigma_{2}=\sum
I(t+\tau)^{2}\Delta t$. This means that unnormalized autocorrelation has
smaller error with respect to the normalized ones.

The standard error (SE) of the autocorrelation function in each lag $\tau$ is
defined as $SE=\frac{\sigma(\tau)}{\sqrt{n(\tau)}}$, where $\sigma(\tau)$ is
the standard deviation of the population in each lag and $n(\tau)$ is the
number of data points in the summation of Eq.(2). The reported upper and lower
limits of the autocorrelation function in each lag are given by the 68%
confidence level as the ${C(\tau)}+({SE}\times 1.69),$ $C(\tau)-({SE}\times
1.69)$, respectively. To estimate the error in $\Delta t$, i.e. $Err(\Delta
t)$, where the autocorrelation function has a maximum at $\Delta t=\tau^{*}$
(the time delay between the two lensed images), we consider the error
corresponding to $\tau^{*}$ where $\tau^{*}=\Delta t=\tau^{*}\pm Err(\Delta
t)$. This time intervals assigns the error in the time delay between the two
lensed images and if the lag error is less than the time resolution, the time
bin is reported as the error in the time delay.

[width=0.7]Methodology/spike_ACF_mock.png

Figure : Autocorrelation function of the mock lensed light curve of Fig.1
(down) that showing a local peak at the imposed time delay of $28$s.

For the burst is shown in Fig.1 (top panel), we generate mock lensed GRB with
an imposed magnification ratio, then decrease time delay between two
gravitational lensed images until the autocorrelation technique can not detect
the time delay between these two images. In Fig.1 the excluded part of figure
(shaded) is the time delay below the sensitivity of autocorrelation technique
in each magnification ratio. In this example, we see that a decrease in the
magnification ratio of the two light curve local peaks makes it more sensitive
to detect shorter time delays. The minimum detectable time delay for a large
magnification ratio is about the width of the peak. As a result, this method
is not sensitive to a delay that is less than the width of each peak in the
light curve (i.e. peak duration).

[width=0.7]Methodology/mintimedelay.jpeg

Figure : An example of the time delay detection sensitivity of the
autocorrelation technique as a function of the magnification ratio, applied
for the simulated lensed light curve of the burst shown in Fig. 1 (top panel)

###  Magnification Ratio

To find the magnification ratio of real lensed GRB light curves, we subtract
the background count rate. The magnification ratio can be sensitive to the
background model. Background intervals were selected before and after the
burst at about the same duration of the burst interval and with little overlap
with the burst interval, then the lowest order polynomial is chosen to fit the
background.

For the simulated light curve, we set the polynomial order to be zero (a
constant). We first remove the background from the simulated count rate and
then the superimposed light curves according to equation (1) (Fig.1 bottom
panel). Retrieved magnification ratio calculated by the ratio of the two local
peaks of this simulated lensed GRB is $R=3.442\pm 0.512$, which approximately
equals the imposed magnification ratio (R = 4).

###  Instrumental Error

Due to purely random distribution of count rate error (which is reported in
ASCII output of RMFIT tool), for each GRB candidate, we generated one hundred
Gaussian realization of its count rate. Then we used the autocorrelation
technique of each realization to find the time delay and the lag error in the
time delay. On the other hand, we computed the magnification ratio and
propagated the error in the magnification ratio for the individual
realizations. Finally, we calculate the mean of these parameters and report
them as the time delay and magnification ratio with the corresponding
$1-\sigma$ uncertainty.

###  Efficiency and Accuracy of the Detection Technique

We performed a Monte-Carlo simulation of mock lensed GRBs to calculate the
efficiency of the Autocorrelation function technique in finding lensed GRBs.
In this simulation, we first selected 81 long GRBs from the Fermi/GBM catalog
which just have one peak in their light curves, so they have no local peaks in
their natural autocorrelation function. Because of using TTE data files of
GRBs light curves which maximally records 300 seconds after the burst, and due
to the set time bin of 64 milliseconds, gravitational time delay less than 64
milliseconds and greater than 300 seconds that correspond to the lens mass of
$M_{L}\leq 10^{2}M_{\odot}$ and $M_{L}\geq 10^{9}M_{\odot}$ are excluded from
our simulation. So we simulated lensed GRBs with a lens mass in the range of
$10^{3}-10^{8}M_{\odot}$. We assumed a point mass lens and defined mass bin of
$3\leq i=\log(M_{L}/M_{\odot})\leq 8$ then in each mass bin of $i$ we
generated 810 mock lensed data ($N_{Generated}=810$) with random parameters of
normalized impact factor $y=\beta/\theta_{E}$ and source and lens cosmological
redshifts ($z_{s}\leq 7$ and $z_{L}\leq z_{s}$) then we applied the
autocorrelation method for these simulated GRBs to evaluate efficiency in
detecting of the imposed time delay by autocorrelation technique. Also we
checked that the retrieved magnification ratio of the simulated light curve is
within the confidence interval of $2\sigma_{R}$ ($\bar{R}-2\sigma_{R}\leq
R\leq\bar{R}+2\sigma_{R}$). If for a simulated lensed GRB we find the
conditions concerning the detection of the imposed time delay and
magnification ratio are fulfilled, we consider it a detected event. We defined
the efficiency of our detection method in $i^{th}$ mass bin as

$\epsilon_{i}=\frac{N_{i}}{N_{Generated}},$ (3)

where $N_{i}$ is detected lensed GRB in each mass bin of $i$. The result of
the simulation is depicted in Fig.1 that shows the autocorrelation technique
has the best performance in the mass about of $10^{6}M_{\odot}$ for the
gravitational lenses.

Efficiency parameter must be used in the calculation of optical depth by
normalizing the number of lensed candidates

$N_{Detected}=\sum_{i}\frac{N_{i}}{\epsilon_{i}}.$ (4)

So the optical depth in Eq.(LABEL:eq:tauTotal) should be rewritten as

$\tau_{obs}=\frac{1}{N_{Total}}\sum_{i}\frac{N_{i}}{\epsilon_{i}},$ (5)

As we will shown in the section LABEL:sec:Results, we find one candidate of
lensed GRB with mass in order of $10^{6}M_{\odot}$, therefore the numerical
value of $\tau_{obs}$ should be
$\tau_{obs}=\tau_{Average}/\epsilon_{6}=(1.340\pm 0.213)f$. We use this result
to find the contribution of gravitational lenses with mass in order of
$10^{6}M_{\odot}$ to the matter content of Universe in the section
LABEL:sec:Results.

To find the accuracy of our method, for detected event in our simulation, we
retrieve the magnification ratio and the time delay from the lensed GRB light
curve, then we use equations (LABEL:eq:delay) and (LABEL:eq:magnification) to
retrieve the mass lens. The histogram of the difference between the determined
mass from the simulated light curve and the imposed lens mass to the imposed
lens mass($\Delta M_{L}/M_{L}$) is depicted in Fig.1 that shows 94.54% of
detected simulated lensed GRBs with our method have $\Delta M_{L}/M_{L}<1$. We
note that uncertainty in mass implies an uncertainty to the equation 5 (_i.e_
$\epsilon_{i}(M_{L})$). Since the efficiency function in Fig. (1) varies in
the logarithmic scale of lens mass, we can ignore our mentioned uncertainty of
mass in the optical depth calculation.

[width=0.7]Efficiency/efficiency.jpeg

Figure : The efficiency of the autocorrelation technique versus the lens mass
in a logarithmic scale as obtained from the Monte-Carlo simulation of
gravitational lensed GRBs

[width=0.7]Efficiency/hist.jpeg

Figure : Histogram of $\Delta M_{L}/M_{L}$ for detected events from the
simulated lensed GRBs. In our method of detection, 94.54% of the detected
simulated lensed GRBs have $\Delta M_{L}/M_{L}<1$