A Riemannian approach for free-space extraction and path planning using catadioptric omnidirectional vision

Image and Vision Computing 95 (2020) 103872

Contents lists available at ScienceDirect

Image and Vision Computing journal homepage: www.elsevier.com/locate/imavis

A Riemannian approach for free-space extraction and path planning using catadioptric omnidirectional vision夽 Fatima Aziz a, b , Ouiddad Labbani-Igbida a, * , Amina Radgui b , Ahmed Tamtaoui b a b

XLIM Institute, UMR CNRS 7252, Université de Limoges, Parc Ester Technopole, 16 rue Atlantis, Limoges 87068, France STRS Laboratory, Institut National des Postes et Télécommunications, 2 Avenue Allal El Fassi, 10100 Rabat, Morocco

A R T I C L E

I N F O

Article history: Received 1 August 2019 Received in revised form 20 December 2019 Accepted 29 December 2019 Available online 19 January 2020 Keywords: Catadioptric vision Riemannian metric Geodesic distance Free-space segmentation Path planning

A B S T R A C T This paper presents a Riemannian approach for free-space extraction and path planning using color catadioptric vision. The problem is formulated considering color catadioptric images as Riemannian manifolds and solved using the Riemannian Eikonal equation with an anisotropic fast marching numerical scheme. This formulation allows the integration of adapted color and spatial metrics in an incremental process. First, the traversable ground (namely free-space) is delimited using a color structure tensor built on the multidimensional components of the catadioptric image. Then, the Eikonal equation is solved in the image plane incorporating a generic metric tensor for central catadioptric systems. This built Riemannian metric copes with the geometric distortions in the catadioptric image plane introduced by the curved mirror in order to compute the geodesic distance map and the shortest path between image points. We present comparative results using Euclidean and Riemannian distance transforms and show the effectiveness of the Riemannian approach to produce safest path planning. © 2020 Elsevier B.V. All rights reserved.

1. Introduction Catadioptric sensors are a class of wide angle imaging systems that have recently attracted high research interest. They offer images with large field of view and are obtained by combining a conventional camera and a convex mirror. The catadioptric system is of great interest because it provides an omnidirectional image in a single shot, and can be used in varied environments and applications such as surveillance and robot navigation. Processing perspective images is based on mathematical concepts and tools like the measurement of lengths, paths, and partial differential operators. Such measurements are usually borrowed from the Euclidean geometry. This is a major issue for catadioptric imaging that contains high geometric distortions and non-uniform resolution. The latter are due to the image formation process, where the 3D scene is projected first on a curved mirror before being captured by a conventional camera (Fig. 1). Setting a catadioptric image as a Riemannian manifold allows us to define a local Riemannian metric that encodes the geometric

夽 This paper has been recommended for acceptance by Sinisa Todorovic. * Corresponding author. E-mail addresses: [email protected] (F. Aziz), ouiddad.labbani-igbida@unilim. fr (O. Labbani-Igbida), [email protected] (A. Radgui), [email protected] (A. Tamtaoui).

https://doi.org/10.1016/j.imavis.2020.103872 0262-8856/© 2020 Elsevier B.V. All rights reserved.

information [3]. This metric should be taken into account in the processing, like the computation of geodesic distances and partial differential operators. Also, setting a color catadioptric image as a surface in a higher dimensional space allows to define a color metric that includes the image multi-components. This paper introduces different Riemannian metrics for central catadioptric omnidirectional imaging systems, each adapted to the desired task. We consider the problem of free-space detection and path planning in the catadioptric image plane and show how it is relevant to use these metrics in computational geometry for robotic applications.

2. Related work Catadioptric vision is still an active research domain for many applications [9,21]. A recent study on imaging models and unwrapping algorithms of catadioptric omnidirectional vision systems are discussed in [12]. While almost omni-image processing use classical Euclidean tools, a few work emphasize the necessity for adapted processing because of the high anamorphosis of these images. The authors in [8] have demonstrated the need to consider the mirror geometry when processing catadioptric images in smoothing and geodesic active contours experiments. They proposed to include the deformation of the mirror by computing an adapted Riemannian metric. More recently, based on the sphere camera model [15], an

2

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

Fig. 1. (a) Catadioptric image. (b) Omnidirectional catadioptric sensor.

approach for computing a generic metric has been proposed in [33]. The main benefit of these approaches is to allow processing directly in the image plane while considering the correction of the geometrical distorsions. In [14], the authors used a catadioptric panoramic system for feature tracking and visual path following tasks. The problem of geometric distortions is solved using a ground unwrapped representation, namely bird’s-eye view of the ground plane. More recently, [27] used an adapted metric, restricted to ground points in omnidirectional images, to produce a skeleton of the local navigable space. The latter is extracted using a catadioptric sensor, based on the algorithm [28] which computes a real-time propagation of active contours to detect free-space boundaries. The authors in [38] have considered the distortion in unwrapping panoramic annular images from a catadioptric imaging system, and proposed a correction method based on the tangential and radial distortion center for panoramic annular images. The method proposed in [22] predicts in real-time the trajectory of multiple moving vehicles at road intersections using a catadioptric omnidirectional camera equipped with an equiangular mirror. For image matching, an important step in the 3D reconstruction, the authors in [17] adapted the matching method in order to deal with non-uniform resolution and geometric distortions in catadioptric images. The approach is based on a geodesic distance used for a new definition of similarity measures. Regarding the color aspect, interesting contributions [32,37] have been proposed for color perspective image processing. However, catadioptric omnidirectional vision is still mainly restricted to grayscale images. A recent work [2] considered this aspect for catadioptric color images by tackling Riemannian Embedding and proposed a generalized hybrid structure tensor. It proved its interest in a smoothing process with an iterative Laplace-Beltrami algorithm. In this paper, the color information and the geometric deformation are both involved in a novel framework to carry out the extraction of the navigable free-space and the path planning task, using a catadioptric sensor. Based on Riemannian geometry tools, mainly the Riemannian Eikonal equation, the problem is solved using an anisotropic fast marching scheme. We propose a three-step process, operating directly in the catadioptric image plane, where each resolution of the Riemannian Eikonal equation requires the choice of an adequate metric.

– First, we solve the problem of free-space detection using a color structure tensor to compute a color distance map. Thanks to this color metric, we delimit the boundaries of the free-space.

– Then, we compute the distance map from the free-space boundaries. The metric used for this purpose is designed in order to correct the distortions due to the mirror of the catadioptric camera during the image formation process. It leads to derive the medial axis of the free-space form, which provides the support for safe path planning. – Thereafter, the path planning requires the computation of the distance map from the goal point. To be the safest one, we exploit at this stage the previously obtained distance map from free-space boundaries as a potential field. We use a gradient descent method to iteratively update the path approximating the geodesic between any given pair of image points in the free-space. The rest of the paper is organized as follows. In Section 3, we present a generic Riemannian metric for catadioptric images using tools of Riemannian geometry. In Section 4, we introduce the distance transform on Riemannian manifolds and its computation based on Hamilton-Jacobi equation, which involves a Riemannian metric. In Section 5, we apply these metrics to solve the free navigable space detection and path planning problems. To demonstrate the efficiency of the proposed approach, comparative results based on Euclidean and adapted Riemannian metrics, using indoor and outdoor catadioptric omnidirectional images, are discussed before concluding in Section 6. 3. Generic Riemannian metric for catadioptric systems After a brief introduction of the unified spherical model used in central catadioptric imaging and some basic tools in Riemannian geometry, we will define adapted metrics for catadioptric processing. 3.1. The equivalent sphere model The sphere camera projection [15] unifies in a single model any type of conventional or central catadioptric systems respecting the constraint of unique effective viewpoint [5]. The central catadioptric geometry can be described as a composition of two mappings on the sphere. The first mapping is a perspective projection of a 3D world point onto a unit sphere, centered in the single viewpoint. The second mapping is a projection onto the image plane from a point on the z-axis of the sphere, at a distance from the projection center that depends on the mirror shape. This model is illustrated in Fig. 2, where the first mapping produces the intersection point Ps on the sphere. It follows a projection from a variable point on the z-axis, defined

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

3

Fig. 2. The unified sphere projection model for central catadioptric image formation.

by the parameter n, to the image point U. The distance between the image plane and the center of the sphere is defined by the parameter v. The unified sphere model has therefore two parameters n and v determined according to the geometric characteristics of the conical section of the mirror. Table 1 gives these parameters where e indicates the eccentricity of the conical section, 4p its latus rectum and f is the focal distance [15]. In [39], the authors proposed an extension to include fisheye models substituting the sphere by a quadratic surface. A later extended version of [15] to perspective cameras with non-linear lens distortions was proposed in [7]. A more recent approach to camera projection modeling and analysis is presented in [19]. The new enhanced unified camera model applies to catadioptric systems and wide-angle fish-eye cameras. This new model has analytic properties that allow it to be easily used in modeling of stereo systems.

3.2. Mathematical background We here introduce some mathematical definitions and concepts, we will use for computing generic spatial-color metrics. We refer the reader to [23] and [29] for more details on Riemannian geometry.

Definition 1. A manifold M is a metric space with the following property: Each point, p ∈ M, has a local neighborhood that is homeomorphic to Rn .

n

Parabolic Hyperbolic & elliptic

1 √

Planar

0

v e e2 +p2

Definition 2. A Riemannian metric on a smooth manifold M is a 2tensor field g that is symmetric (i.e. g(X, Y) = g(Y, X)) and positive definite (i.e. g(X, X) > 0 if X = 0). A Riemannian metric determines an inner product on each tangent space Tp M, which is typically written < X, Y >:= g(X, Y) for X, Y ∈ Tp M. A manifold M together with a given Riemannian metric g is called a Riemannian manifold and denoted (M, g). For instance, on Rn , the Riemannian metric is given by the standard inner product g(X, Y) = X.Y for all X, Y ∈ Tp Rn , for all p ∈ Rn . We call Rn with this Riemannian metric, the Euclidean space. Definition 3. A Riemannian embedding is a smooth map f : M → N between manifolds (M, glm ) and (N, hij ) such as glm = f ∗ hij , i.e. g(X, Y) = h(Df(X), Df(Y)) for all tangent vectors X, Y ∈ Tp M and all p ∈ M. Definition 4. Considering a manifold (M, glm ), the squared arc length ds2 is given by the metric glm on this manifold, ds2 = glm dxl dxm

1 − 2p e(1−2p) √

using the Einstein summation convention1 .

Let us consider the following embedding f that describes the mapping between the two Riemannian manifolds (R2 , glm ) and S2 ⊂

e2 +p2

f

(1)

3.3. Generic metric design

Table 1 The parameters n and v depending on the mirror shape. Mirror

Examples of two-dimensional manifolds (called surfaces) include the Euclidean plane R2 and the sphere S2 .

1

Identical indices that appear one up and one down are summed over.

4

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

(R3 , hij ) representing the image plane and the sphere, equipped respectively with Riemannian metrics glm and hij : R2 → S2 ⊂ R3   (x, y) → xs = w1−1 (x, y), ys = w2−1 (x, y), zs = w3−1 (x, y)

f :

(2)

(3)

−→ NU

(4)

with the coordinates N = (0, 0, n), Ps = (xs , ys , zs ), U = (x, y, v). This projection leads to the following equations: ⎧ ⎪ ⎪ ⎨xs = ax ys = ay ⎪ ⎪ ⎩z = a(v − n) + n

 



− 1 − n2 xy (n − v)2 + 1 − n2 x2

  2

2 2 2 − 1 − n xy (n − v) + 1 − n y

(10)

where,

where l, m = {1, 2} (referring to x, y-spatial variables) and i, j = {1, 2, 3}. ∂l f i = ∂xl f i is the partial derivative of fi with respect to xl variable. The projection w : S2 → R2 from the point N, belonging to the z-axis, maps a point Ps on the sphere to a point U in the image plane. The line through N and Ps intersects the hyperplane z = v as shown in Fig. 2. For a non-zero a, the image point U is characterized by, −→ NPs = a

glm = g

where f(x, y) = w −1 (x, y) = (xs , ys , zs ) maps a point from the catadioptric image plane to the sphere. Once the embedding f is constructed, the metric glm over R2 can be induced by hij as the pullback metric defined as a symmetric 2tensor: glm = hij ∂l f i ∂m f j

And consequently, the inverse metric tensor is,

(5)



g=

2 (n − v)2 + x2 + y2

 

2

(n − v)2 n(n − v) + (n − v)2 + 1 − n2 x2 + y2

(11)

As shown on Fig. 3, that illustrates the metric coefficients in the case of a hyperbolic mirror, this tensor field is anisotropic and depends on the position of the point in the image. The Riemannian metric allows distances and angles to be measured. We will use it in the following to compute the geodesic distance map from the boundaries of the free-space and to derive the shortest path between two points. 4. Riemannian distance transform As catadioptric images are reinterpreted as Riemannian manifolds, the distance transform should be posed in a Riemannian framework managed by the induced metric. We hereafter define the geodesics in the catadioptric image plane and a numerical approximation scheme for their computation. 4.1. Geodesic distance function Geodesics are locally length-minimizing curves that define the shortest path between points in the space. They are naturally induced by the Riemannian metric on the considered manifolds.

s

Being located on the surface of the unit sphere, the point Ps ∈ S2 satisfies: x2s + y2s + zs2 = 1

(6)

By substituting Eq. (5) into Eq. (6), we get:

a± =

n(n − v) ±



(n − v)2 + (x2 + y2 )(1 − n2 ) (n − v)2 + x2 + y2

(7)

Inserting a into Eq. (5) allows to obtain the coordinates (xs , ys , zs ) of a point belonging to the surface of the sphere according to the (x, y) spatial coordinates, and equivalently the formula of the inverse projection w −1 with respect to (x, y) and the parameters n and v. Considering the metric hij on the surface of the sphere as locally Euclidean, we can deduce the pullback Riemannian metric on the catadioptric image plane using the definition in Eq. (3). After calculations, the components of the Riemannian metric tensor take on the form of the following symmetric matrix, glm = f





1 − n2 xy (n − v)2 + 1 − n2 y2

  2

2 2 2 1 − n xy (n − v) + 1 − n x

(8)

where,

 

2

n(n − v) + (n − v)2 + 1 − n2 x2 + y2 f=





2

(n − v)2 + 1 − n2 x2 + y2 (n − v)2 + x2 + y2

Definition 5. Let (M,G) be a Riemannian manifold equipped with a Riemannian metric G = {glm }l,m=1. . . dim(M) . A curve c : I ⊆ R → R ˙ is parallel along parameterized by t ∈ I is a geodesic if and only if c(t) ˙ c(t), that is, at every point p ∈ c(t) the covariant derivative of c(t) is ˙ zero: ∇ c c(t) = 0. The most familiar examples of geodesics are the straight lines in the Euclidean geometry and geodesics on the sphere. In the case of a spherical surface, geodesics are great circles, i.e., the shortest path between any two points A and B is the arc of a circle passing through A and B and whose center and radius are the ones of the sphere (Fig. 4). Definition 6. Given a Riemannian manifold (M,G) and c(t) the geodesic joining two points p, q ∈ M. The smallest Riemannian distance dG (p, q) is defined as: ⎧T ⎫  ⎨ ⎬ c(0) = p ˙ ˙ )dt ; s.t. G (c(t), dG (p, q) = inf c(t) ⎩ ⎭ c(T) = q

(12)

0

One can fix a set of points S ∈ M and generalize the definition above to the distance map to this set. Definition 7. The Riemannian distance map to a set of points S ∈ M is defined as

(9)

DS (p) = inf dG (p, q) q∈S

∀p ∈ M

(13)

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

5

Fig. 3. Metric components {glm }l,m=1,2 = {g11 , g12 , g22 } on the image plane in case of a hyperbolic mirror with parameters n = 0.93 and v = 1.

The problem of optimal distance computation can be transformed into a problem solving a partial differential equation. It was shown [4] that the distance map DS is the unique weak solution, that is usually called the viscosity solution [10], of the time-dependent Hamilton-Jacobi equation:

∂ u(p, t) = ∇ G u G ∂t

(14)

√ where ∇ G u G = ∇uT G−1 ∇u is the norm of the gradient in the Riemannian metric space (see [25] for a detailed mathematical view of the Hamilton-Jacobi equation and the distance function on Riemannian manifolds). The distance DS from the set S can also be considered as the arrival time of a front, propagating from S and moving at velocity controlled by the metric G. In this way, the evolutive problem (14) can be replaced by a stationary Hamilton-Jacobi equation, called the Eikonal equation: ∇G

x G = 1

(15)

where x represents the time at which the front reaches the point p ∈ M. The metric G plays an important role and may take different forms depending on the application. In the specific case of a weighted Euclidean metric by a scalar function F where G = F2 dij , one can recover the well-known Euclidean Eikonal equation: ∇x = F

(16)

We denoted here by . the standard Euclidean norm. 4.2. Anisotropic fast-marching Different fast algorithms have been proposed in order to approximate numerically the solution of the Euclidean version of the Eikonal Eq. (16). The local solvers can be classified into three categories. The first approach is to use upwind finite differences [20,34,35]. The second one is to use dynamic programming techniques, based on a semi-Lagrangian numerical scheme where the problem is recast as an optimal control problem [6]. The last strategy is to use finite element techniques [24] (see [1] for an extensive comparison). One of the early attempts to deal with the Hamilton-Jacobi equation within a framework of anisotropic propagation is due to [4,18], generating anisotropic meshes and automatic refinement in some particular regions. Later, [13] derived a Riemannian level set formulation for image segmentation using the mean curvature flow, where the induced metric is related to the local structure tensor of an image. The most famous solver of the Eikonal equation is the FastMarching (FM) algorithm introduced by Sethian [34]. It operates similarly to the Dijkstra algorithm for finding the shortest path on a network. Kimmel and Sethian [20] have developed a version of the Fast-Marching algorithm for a surface M ⊂ Rn with a specific metric. Their algorithm operates on triangulated meshes, and treats the triangles as being locally flat and equipped with an isotropic metric (Fig. 5). The same algorithm was used [30,31] with anisotropic metrics, to compute geodesic distances and paths for many applications as surface and shape processing, mesh generation, and shape

Fig. 4. Example of geodesics: The straight line between points A and B in the Euclidean plane (left), and the arc between points A and B of a great circle on the sphere surface (right).

6

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

Fig. 5. (a) Triangulated mesh. (b) Example of propagation using the Fast-Marching algorithm on triangulated domains. Source: Source [30].

comparison. We use this algorithm, with an optimal finite difference scheme on triangulated domains, in order to solve the Riemannian Eikonal equation directly in the catadioptric image plane considered as a Riemannian manifold.

robot equipped with a paracatadioptric omnidirectional sensor that captures color images with a resolution of 800 × 600, while navigating in an outdoor real-world environment. 5.2. Free-space form

5. Free-space extraction and path planning We consider here the problem of free-space detection and distance map construction using catadioptric imaging only. As we will see, the considered Riemannian framework yields to a global setting for tailor-made metrics depending on the application. For a global overview of the method, Algorithm 1) summarizes the three-steps process for Riemannian distance maps and geodesic path generation.

Algorithm 1. Path planning using catadioptric color images.

The first step is to find where the robot is free to move, namely, the free-space. The model is based on the assumption of the road to be approximately of constant intensity value in the image. The edge map allows in this way to delimit the free-space. One can thus initiate a circular curve around the dark spot in the center of the image (which is the projection of the support that holds the mirror above the lens) and propagate it. This propagation is described by the following Riemannian Eikonal equation:

∇xT G−1 c ∇x = 1

(17)

where Gc denotes the Riemannian metric of the space. We should carefully select Gc in order to obtain a distance map that yields to extract the free-space. As the edges permit delimiting this free-space, we choose to build Gc from the classical structure tensor. The latter was initially introduced by Di Zenzo in [11] and used to analyze the local structure in the image and to extract interest features. For a real-valued function Ii , with i = {1, . . . , m}, representing an image of m channels, the structure tensor could be obtained as the pullback metric of the following mapping: f : R 2 → Rm

(18)

(x, y) → (I1 , . . . , Im ) 5.1. Experimental setup and test images For evaluation purposes, the distance maps and the generated paths obtained with standard Euclidean and Riemannian metrics will be compared using synthetic and real catadioptric images acquired in indoor and outdoor environments with various topologies like corridors and crossroads. Synthetic catadioptric images were obtained by defining virtual scenes and including a hyperbolic catadioptric sensor in the POV-Ray vision ray-tracer simulator2 with the objective of generating images presenting various geometric topologies (corridors, intersections) of the free space. We also used real test images that we extracted from the dataset [26]3 acquired by a mobile

2

http://www.povray.org/ 3 The dataset of catadioptric images is available at https://box.xlim.fr/d/ 188e4048714c4f10b597/

We then deduce the expression of the classical structure tensor given by the symmetric matrix, m i i  (Ixi )2 Ix Iy i=1 m i i i 2 I I i=1 x y i=1 (Iy )



m

T=

i=1 m

(19)

It can be decomposed as a weighted sum of its eigenvectors (v1 , v2 ), T = k1 v1 vT1 + k2 v2 vT2

(20)

Considering k1 ≥ k2 , the direction of k1 , namely v1 , is maximally aligned with the gradient and indicates the prominent local orientation. Fig. 6 represents the eigenvalues k1 and k2 of the structure tensor for an example of omnidirectional image.

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

7

Fig. 6. Eigenvalues k1 and k2 of the structure tensor (19) for an omnidirectional image.

Since we want to underline where the evolving curve crosses the boundary of the road, one can construct a metric Gc from the structure tensor T that would serve as edge detector, by Gc = k1 v1 vT1

(21)

The metric Gc incorporates the local gradient information and the problem is formulated using the Riemannian Eikonal Eq. (17). The viscosity solution of Eq. (17) is the geodesic distance map Dc . This map is approximately zero in homogeneous flat regions (mainly the road), and tends to have large values along the edges. In this way, the structure tensor is used as a metric that leads to measure the distance between colors. The boundary of the free-space is then extracted using the geodesic distance map Dc and a specified threshold. One main advantage of this approach is to define a unified formulation to compute distance transform that efficiently includes the multidimensional aspect of the image. Fig. 7 shows, for synthetic and real catadioptric images, the color-metric distance map and the free-space depicted as white regions. 5.3. Distance transform and geodesic paths Optimal path planning consists in determining the shortest path between two given configurations while respecting a given number of criteria and constraints. The A∗ algorithm is the most famous to solve this problem. It is known for its simplicity and privileged for its speed. A recent alternative to the algorithm A∗ is based on the resolution of the Eikonal equation. For this purpose, different numerical schemes have been proposed, such as, the Fast Marching Method (FMM), the Fast Sweeping Method (FSM), and the Fast Iterative Method (FIM).Their main advantage is to offer a good quality trajectory, smoother and more accurate. Vision-based path planning has been largely explored for perspective images. However, only a few works deal within the framework of catadioptric omnidirectional vision, and almost all methods are based on algorithms developed for conventional images. In the following, we introduce a new formulation of the distance transform from free-space boundaries, and how it will be used in path planning between any two points in the free-space domain of catadioptric images. 5.3.1. Distance transform from free-space boundaries Given the free-space form, each point of this domain is assigned a value equal to its distance to the nearest point belonging to the freespace boundary. This process of computing the distance transform

is usually based on the Euclidean metric and produces a distance map which values are zero on the boundaries and increase when moving far from it. The notion of distance from the boundaries in this case is obviously altered by the projective deformations of the mirror. We propose to solve again the Eikonal Eq. (22) based on an omnidirectional Riemannian metric defined by,

∇xT G−1 s ∇x = 1

(22)

where Gs is the catadioptric omnidirectional metric induced from the mirror geometry (8). This leads to the distance transform depicted Fig. 8c for some examples of catadioptric scenes. In order to show the efficiency of Riemannian processing on catadioptric images, we compare the Euclidean and Riemannian distance functions from the boundaries of the free-space, extracted by Eq. (17), by solving respectively the Euclidean and the Riemannian formulations of the Hamilton-Jacobi Equation. Comparative results are shown (Fig. 8a & b). In order to understand these results, we transform them into bird’s-eye views in Fig. 8c & d, where the remapped images are the orthographic projections of the ground plane. These views are convenient to the human vision system and allow better visualization. We note that the distance maps and their corresponding bird’seye views of the two methods, are different. At first glance, one can see that for the Euclidean metric, the output appears to match the expected distance from the boundary of the free-space shape. However, the projection into the ground plane produces a view that does not fit the real environment. While, for the Riemannian method, the statement is different. The bird’s-eye views prove that the Riemannian metric induced from the mirror geometry of the catadioptric camera, provides correct and more accurate results than the Euclidean ones. 5.3.2. Obstacle-free geodesic paths Classical FM methods provide the shortest path in terms of Euclidean distance, which could lead to risky situations because of its closeness to obstacles (see examples of distance maps in Fig. 8). The authors in [36] have proposed the FM Square (FM2 ) algorithm in order to solve this problem and generate reliable and safe paths. The idea behind this method is to compute two front expansions over a gridmap. The first FM expansion allows to obtain a first potential field which is the distance map from the boundaries of the obstacles. This potential field is then used to compute the second expansion,

8

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

Fig. 7. Free-space detection using anisotropic front propagation: (a) The original images followed by the photometric distance map (b), and finally, the corresponding extracted free-spaces (c).

from the target point to the initial point, generating a path following the maximum gradient direction. Motivated by the robustness and the simplicity of the FM2 method, we adapt it to catadioptric images and reset it in a Riemannian framework. The first potential field F1 , in our case, is the Riemannian distance function obtained as the viscosity solution of the Riemannian Eikonal Eq. (22). This distance map F1 is used to compute a second potential F2 which is the distance from the target point to the initial point. It is obtained by solving the equation:

∇xT (F12 dij )−1 ∇x = 1

(23)

The path c that the wavefront has followed from the origin (x0 ) to the target point (xg ) will be the shortest path in time. It is computed using a gradient descent method considering that the traveling speed along the path is the expansion velocity of the wavefront.  ∂ c(t) ∂t

∇F

) = − ∇F2 ((c(t))

c(0) = xg .

c(t)

2

∀t > 0

(24)

As the front expansion originates from the target point, the field F2 will have only one global minima at the target point. Hence, following the maximum gradient direction from the initial point, we

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

9

Fig. 8. Comparison between Euclidean (a) and Riemannian (b) distance maps and their corresponding bird’s eye views (c) and (d) respectively, for synthetic and real catadioptric images with free spaces of different topologies (corridors and crossroads).

will reach the originated target point, in a unique way, obtaining the trajectory. In order to analyze the generated paths, we consider, in Figs. 9 and 10, free spaces of simple shapes. In Fig. 9a, the free-space is a disk on a para-catadioptric image. Following the proposed Algorithm 1, we compute geodesics applying the Euclidean (Fig. 9b & c) and the proposed Riemannian (Fig. 9d & e) metrics respectively. It is well

known [16] that for a para-catadioptric image, a 3D straight line (which is the shortest distance in Euclidean space) projects as a circle on the catadioptric image plane. Hence, as expected, the geodesics obtained with the Riemannian metric for para-catadioptric images, are circular and reflect the shortest paths. The other example Fig. 10 depicts a corridor free-space projected on a para-catadioptric image. The geodesic path computed following Algorithm 1 is superimposed

Fig. 9. Comparison of trajectories obtained using Euclidean and Riemannian metrics on a disk-shaped free space centered in a synthetic para-catadioptric image. The Euclidean path is superimposed on the distance map from the boundaries of the free space (b) and on the distance map from the target point (c). The geodesic obtained with the proposed Riemannian metrics is superimposed on the distance map from the boundaries of the free space (d) and on the distance map from the target point (e).

10

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

Fig. 10. Path planning in a corridor-shape free space. The geodesic path is superimposed on the distance maps, computed respectively: (b) from the boundary of the free space (the second step in Algorithm 1) and (c) from the target point (the third step in Algorithm 1).

on the distance maps from the free-space boundaries and from the target point respectively. One can see that the path lies on the medial axis of the free-space form, the latter being characterized by the maxima (peaks) of the distance transform. As a result, this trajectory is safe and optimal in time. The proposed approach leads to compute paths that tend to navigate far from obstacles.

These results are confirmed by experiments using synthetic and real catadioptric images where a comparison of the computed paths between different pairs of points is demonstrated in Fig. 11. One can see that using the Euclidean metric, the generated trajectory is close to the boundary of obstacles, while it is expected to get away from them as far as possible. However, with the Riemannian metric, the

Fig. 11. Comparison of the generated paths between two points using Euclidean (a) and Riemannian (b) Fast Marching Square methods. The corresponding bird’s-eye views are shown respectively in (c) and (d) with the computed paths.

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

11

Fig. 12. Comparative results of the minimal distances along the paths using the Euclidean and the Riemannian metrics in the FM2 scheme. The Euclidean path (in blue) and the Riemannian one (in red) are superimposed on the original images. To be noted that the lengths of the paths being different, the path progression is normalized within [0, 1]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

produced path is optimal and short. It leads to a better safety distance from obstacles according to different roads topologies and situations. The Euclidean and Riemannian paths are embedded in distinctive metric spaces. In order to provide a more quantitative assessment of the paths’ safety, we projected both paths on the Riemannian distance maps, mirroring the real geometry of the sensor. Fig. 12 depicts the minimal distances to the closest obstacles along each path obtained with the FM2 scheme built on the Euclidean and Riemannian metrics respectively. We note that in all cases, the geodesic path presents higher distances to obstacles. This is also confirmed by computing the average of the distance map for several paths for the two methods (Fig. 13), showing that greatest safety distances are achieved by the Riemannian path planning. We deduce that the Riemannian adaptation is not an optional choice that could only outperform the conventional Euclidean metric-based algorithms; but it is a mandatory measure for ensuring safe and optimal tasks in mobile robot applications using catadioptric imaging. To be noted that the potential field computed by the proposed FM2 can also be used to compute a velocity potential field and to adapt the robot velocity according to its distance to free-space boundaries (obstacles), which provides a safe speed for the robot at any point in the environment.

6. Conclusion We have introduced in this paper a Riemannian framework to handle color catadioptric imaging for free-space detection and path planning applications. We have demonstrated that the Riemannian approach for color catadioptric images overcomes the limitations of the classical Euclidean one, and gives a powerful alternative to the classical calculus. In particular, the distance transform computation is recasted to involve an adequate Riemannian metric. The choice of the metric is specific to each task. It was, in our case, a color-metric for performing an improved free-space extraction, and a spatial-metric, induced from the mirror geometry, for computing Riemannian distance maps and geodesics paths. Hence, one of the main advantages of this framework lies in its generality and flexibility, which allows to extend the classical algorithms to manifolds and to endow the space with different metrics. Additionally, the processing is done on the catadioptric plane directly without need of unwrapping the images to any other representation. Experiments of free-space detection and optimal trajectory generation, using real and synthetic color catadioptric images, acquired in indoor and outdoor environments, showed promising results for mobile robot navigation. Our ongoing work deals with the complexity optimization of the algorithm to provide real-time trajectories and to extend the approach to 3D free spaces and robot navigation.

The distance map from obstacles average for several images

0.045

Euclidean method Riemannian method

CRediT authorship contribution statement

0.04

Fatima Aziz: Conceptualization, Methodology, Formal analysis, Software, Writing - original draft, Writing - review & editing. Ouiddad Labbani-Igbida: Conceptualization, Supervision, Writing original draft, Writing - review & editing, Funding acquisition, Project administration. Amina Radgui: Supervision, Writing - original draft. Ahmed Tamtaoui: Funding acquisition, Resources.

0.035

0.03

0.025

0.02

Declaration of competing interest 0.015

The authors declare no conflict of interest. 0.01

References 0.005

0 1

2

3

4

5

6

7

8

9

10

11

Fig. 13. Comparison of the average of the distance map to the closest obstacles for several paths using the Euclidean metric (in blue) and the Riemannian one (in yellow). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

[1] S. Ahmed, S. Bak, J. McLaughlin, D. Renzi, A third order accurate fast marching method for the Eikonal equation in two dimensions, SIAM J. J. Sci. Comput. 33 (5) (2011) 2402–2420. [2] F. Aziz, O. Labbani-Igbida, A. Radgui, A. Tamtaoui, Color-metric tensor for catadioptric systems, Image Processing (ICIP), IEEE Int. Conf. on, IEEE. 2016, pp. 1594–1598. [3] F. Aziz, O. Labbani-Igbida, A. Radgui, A. Tamtaoui, Generic spatial-color metric for scale-space processing of catadioptric images, Comput. Vis. Image Underst. 176–177 (2018) 54–69.

12

F. Aziz, O. Labbani-Igbida, A. Radgui, et al. / Image and Vision Computing 95 (2020) 103872

[4] P. Bagnerini, P. Hoch, M. Rascle, The Eikonal equation on a manifold. Applications to grid generation or refinement, Hyperbolic Problems: Theory, Numerics, Applications, Springer. 2001, pp. 109–118. [5] S. Baker, S.K. Nayar, A theory of single-viewpoint catadioptric image formation, Int. J. Comput. Vis. 35 (2) (1999) 175–196. [6] M. Bardi, I. Capuzzo-Dolcetta, Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations, Springer Science & Business Media. 2008. [7] J.P. Barreto, A unifying geometric representation for central projection systems, Comput. Vis. Image Underst. 103 (3) (2006) 208–217. [8] I. Bogdanova, X. Bresson, J.-P. Thiran, P. Vandergheynst, Scale space analysis and active contours for omnidirectional images, IEEE Trans. Image Process. 16 (7) (2007) 1888–1901. [9] C. Cai, X. Weng, Q. Zhu, Sea-skyline-based image stabilization of a buoy– mounted catadioptric omnidirectional vision system, EURASIP J. Image and Video Process. (2018). [10] M.G. Crandall, P.-L. Lions, Viscosity solutions of Hamilton-Jacobi equations, Trans. Am. Math. Soc. 277 (1) (1983) 1–42. [11] S. Di Zenzo, A note on the gradient of a multi-image, Comput. Vis. Graphics Image Process. 33 (1) (1986) 116–125. [12] J. Dong, H. Yu, L. Kong, Research on imaging model and unwrapping algorithm of catadioptrio-mnidirectional vision system, 2018 IEEE International Conference on Mechatronics and Automation (ICMA), IEEE. 2018, pp. 861–865. [13] R. San José Estépar, S. Haker, C.-F. Westin, Riemannian mean curvature flow, Int. Symp. on Visual Computing, Springer. 2005, pp. 613–620. [14] J. Gaspar, N. Winters, J. Santos-Victor, Vision-based navigation and environmental representations with an omnidirectional camera, IEEE Trans. Robot. Autom. 16 (6) (2000) 890–898. [15] C. Geyer, K. Daniilidis, A unifying theory for central panoramic systems and practical implications, European Conf. on Computer Vision, Springer. 2000, pp. 445–461. [16] C. Geyer, K. Daniilidis, Catadioptric projective geometry, Int. J. Comput. Vis. 45 (3) (2001) 223–243. [17] I. Guelzim, A. Amkoui, L. Koutti, Matching of omnidirectional images based on the geodesic distance, 2016 IEEE/ACS 13th International Conference of Computer Systems and Applications (AICCSA), IEEE. 2016, pp. 1–5. [18] P. Hoch, M. Rascle, Hamilton-Jacobi equations on a manifold and applications to grid generation or refinement, SIAM SIAM J. Sci. Comput. 23 (6) (2002) 2055–2073. [19] B. Khomutenko, G. Garcia, P. Martinet, An enhanced unified camera model, IEEE Robot. Autom. Lett. 1 (1) (2015) 137–144. [20] R. Kimmel, J.A. Sethian, Computing geodesic paths on manifolds, Proc. Natl. Acad. Sci. 95 (15) (1998) 8431–8435. [21] Y-J. Ko, S-Y. Yi, Catadioptric imaging system with a hybrid hyperbolic reflector for vehicle around-view monitoring, J. Math. Imaging Vision 60 (4) (2018) 503–511. [22] V. Krishnamoorthy, S. Agarwal, K. Venkatesh, Vehicle trajectory prediction using a catadioptric omnidirectional camera, 2016 International Conference on Advances in Computing, Communications and Informatics (ICACCI), IEEE. 2016, pp. 2761–2764.

[23] J.M. Lee, Riemannian Manifolds: An Introduction to Curvature, 176. Springer Science & Business Media. 2006. [24] F. Li, C.-W. Shu, Y.-T. Zhang, H. Zhao, A second order discontinuous Galerkin fast sweeping method for Eikonal equations, J. Comput. Phys. 227 (17) (2008) 8191–8208. [25] Mantegazza, Mennucci, Hamilton-Jacobi equations and distance functions on Riemannian manifolds, Appl. Math. Optim. 47 (1) (2002) 1–25. [26] R. Marie, O. Labbani-Igbida, P. Merveilleux, E.M. Mouaddib, Autonomous robot exploration and cognitive map building in unknown environments using omnidirectional visual information only, 2013 IEEE Workshop on Robot Vision, WORV 2013, 2013. pp. 191–196. [27] R. Marie, O. Labbani-Igbida, E.M. Mouaddib, Scale space and free space topology analysis for omnidirectional images, Robotics and Automation (ICRA), 2014 IEEE Int. Conf. on, IEEE. 2014, pp. 4451–4456. [28] P. Merveilleux, O. Labbani-Igbida, E.M. Mouaddib, Real-time free space detection and navigation using omnidirectional vision and parametric and geometric active contours., Robotics and Automation (ICRA), IEEE Int. Conf. on, IEEE. 2011, pp. 6312–6317. [29] P. Petersen, Riemannian Geometry, vol. 171. Springer. 2006. [30] G. Peyré, L.D. Cohen, Geodesic methods for shape and surface processing, Advances in Computational Vision and Medical Image Processing, Springer. 2009, pp. 29–56. [31] G. Peyré, M. Péchaud, R. Keriven, L.D. Cohen, et al. Geodesic methods in computer vision and graphics, Found. Trends Comput. Graph. Vis. 5 (3-4) (2010) 197–397. [32] K. Plataniotis, A.N. Venetsanopoulos, Color Image Processing and Applications, Springer Science & Business Media. 2013. [33] L. Puig, J.J. Guerrero, Scale space for central catadioptric systems: towards a generic camera feature extractor, Computer Vision (ICCV), IEEE Int. Conf. on, IEEE. 2011, pp. 1599–1606. [34] J.A. Sethian, A fast marching level set method for monotonically advancing fronts, Proc. Natl. Acad. Sci. 93 (4) (1996) 1591–1595. [35] J.A. Sethian, et al. Level set methods and fast marching methods, J. Comput. Inf. Technol. 11 (1) (2003) 1–2. [36] A. Valero-Gomez, J.V. Gomez, S. Garrido, L. Moreno, Fast marching method for safer, more efficient mobile robot trajectories, IEEE Robot. Autom. Mag. 1070 (9932/13), (2013). [37] J. Van De Weijer, T. Gevers, A.W. Smeulders, Robust photometric invariant features from the color tensor, IEEE Trans. Image Process. 15 (1) (2006) 118–127. [38] X. Wang, X. Bu, B. Yang, G. Li, J. Yu, Corrected unwrapping method based on the tangential and radial distortion centre for the panoramic annular image, IET Image Process. 9 (2) (2014) 127–133. [39] X. Ying, Z. Hu, Can we consider central catadioptric cameras and fisheye cameras within a unified imaging model, European Conf. on Computer Vision, Springer. 2004, pp. 442–455.