Liuxin Zhang
email: zhangliuxin@bit.edu.cn

Yumo Yang
email: yangyumo@bit.edu.cn

Yunde Jia
email: jiayunde@bit.edu.cn

State-Based Visibility for 3D Reconstruction from Multiple Views

Estimating visibility is one of the most important steps for multi-view reconstruction using volumetric scene representation. In this paper, we propose a simple approach to estimating visibility based on the current state of a scene which is implicitly represented as the zero level set of a function. Our method can determine the regions of a 3D space visible to several given viewpoints efficiently. Using this visibility method, we also present a new variational formulation for multi-view reconstruction. We cast the multi-view reconstruction problem as an optimization of a novel energy functional amenable for minimization with an Euler-Lagrange driven evolution. The proposed algorithm has been applied to both synthetic and real datasets with promising results.

Introduction

Volumetric multi-view stereo reconstruction, originally introduced by Seitz group [START_REF] Seitz | Photorealistic scene reconstruction by voxel coloring[END_REF][START_REF] Kutulakos | A theory of shape by space carving[END_REF], has recently been shown to produce 3D models from photographs or video sequences with fairly high quality [START_REF] Vogiatzis | Mult-view stereo via volumetric graph-cuts[END_REF][START_REF] Esteban | Stereo and silhouette fusion for 3D object modeling from uncalibrated images under circular motion[END_REF]. The basic principle in volumetric reconstruction is to find a classification for all elements (voxels or grid points) within a discrete volume space whether they belong to the surface of a 3D object or not.

Basically, the surface of a 3D object is nearly Lambertian, i.e., only elements belonging to this surface should have a consistent appearance in the input images (be photo-consistent). Therefore, the most central aspect of all the existing volumetric multi-view stereo reconstruction techniques is the estimation of the so called photo-consistency of a given element.

We assume that all the views (cameras) are well calibrated within the global world coordinate system, i.e., for each point in the world space, it is possible to determine the coordinates of its projection onto each view. In this case, before estimating the photo-consistency of a given point, we need to compute the visibility of this point with respect to all the views. This visibility problem can be described as follows: given a state of a scene, determine whether a point is occluded by other points from the view. This problem is one of the most important steps for multi-view reconstruction based on a volumetric scene representation.

In this paper, we propose a simple approach to visibility estimation based on the current state of a scene which is implicitly represented as the zero level set of a function. Our method can obtain the visibility information of each point in two or three dimensional space to all the views through an efficient recursion strategy. Using this visibility estimation, we also present a novel model for 3D reconstruction from multiple views in a variational framework. In contrast with some other multi-view reconstruction methods, the main benefit of our method is its independence from initialization.

In the following, we first review related work in Section 2 and discuss our visibility algorithm in Section 3. In Section 4, we describe the reconstruction method and its conversion to an optimization problem with efficient implementation by an Euler-Lagrange driven evolution. The experimental results are shown in Section 5 to validate our proposed approach and conclusions are made in Section 6.

Related Work

The overview of state-of-art algorithms for multi-view stereo reconstruction reported by Seitz [START_REF] Seitz | A comparison and evaluation of multi-view stereo reconstruction algorithms[END_REF] shows a massive domination of volumetric method. Nearly all of them use a photo-consistency measure such as normalized cross correlation (NCC) or sum of squared differences (SSD) to evaluate how consistent a reconstruction is with a set of images. Visibility plays a very important role in volumetric reconstruction because the only requirement to compute this photoconsistency is that camera visibility is available.

Lempitsky et al. [START_REF] Lempitsky | Oriented visibility for multiview reconstruction[END_REF] state that visibility estimation is arguably the most difficult problem in 3D reconstruction from multiple views. To estimate the true visibility of some surface element, one needs to know the true scene geometry and vice versa. Currently there are numerous algorithms for approximately solving this chicken-and-egg problem, which can be mainly divided into two categories: state-based visibility and oriented visibility. State-based visibility can be considered as a type of global visibility algorithm which determines the visibility information by the state of the whole scene. It is often used in the approaches which reconstruct the scene geometry based on iterative process [START_REF] Kutulakos | A theory of shape by space carving[END_REF][START_REF] Faugeras | Variational principles, surface evolution, PDE's, level set methods, and the stereo problem[END_REF], and at each moment of time, a point is considered visible to a viewpoint if it is not occluded with current scene configuration (Fig. 1 left). By contrast, oriented visibility is a type of local visibility algorithm. The visibility of a patch to a viewpoint is estimated by using the patch position and normal direction. In oriented visibility, a local patch is considered visible to a viewpoint if and only if it locates in front of this viewpoint according to a predefined angle from the normal direction (Fig. 1 right). Although this local method is simple, sometimes when the global scene geometry is more complex, such as a surface with high curvature (Fig. 2 left) or a surface with several unattached parts (Fig. 2 right), it is not reliable to compute the visibility of a local patch only from its orientation because of the self-occlusion.

Our visibility algorithm falls into the first category, and we take the global scene geometry to obtain more reliable visibility information. Our work is similar to that of Tsai et al. [START_REF] Tsai | Dynamic visibility in an implicit framework[END_REF], who describe an implicit ray tracing technique to X S X N f Fig. 1. Two approaches to visibility approximation (Lempitsky et al. [START_REF] Lempitsky | Oriented visibility for multiview reconstruction[END_REF]). On the left is a state-based visibility, the visibility information of a point X with respect to all the views is determined by the current global scene geometry S. On the right is an oriented visibility, a local patch (X, N) is considered visible from the viewpoints within a predefined angle φ from the normal direction N.

determine the visibility of a single viewpoint in two or three dimensional space. In their algorithm, a Cartesian grid is used to discretize the space and the viewpoint is supposed to lie on a grid node. They organize the entire grid in a special sweeping sequence so that each grid node can obtain its own visibility information to the viewpoint without violation. Although their method can solve the visibility problem of a single viewpoint quickly and efficiently, it can not be used directly in the field of multiple views. It is impossible for us to establish a uniform Cartesian grid to ensure all the views being on the grid nodes at the same time. In our work, we organize all the grid nodes in a more reasonable way that does not need the viewpoint to lie on the grid node. So, our algorithm is able to be used in any case including multiple views. 3 State-Based Visibility in an Implicit Framework

The Level Set Method and Implicit Representations

The level set method proposed by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] has found many applications in recent years. The basic idea behind this method can be thought of using a higher-dimensional representation for the problem of front tracking. This method provides us with a robust representation of even a more complex surface and its evolution.

We assume that the whole scene (or its part we are interested in) is located within some bounding volume Ω ⊂ ℜ 3 . Each allowed scene configuration is characterized by some occupied subvolume M ⊂ Ω with the piecewise-smooth boundary S = ∂M (the scene surface). We define a Lipschitz smooth function φ : ℜ 3 → ℜ in Ω, which represents the scene surface as the zeros of φ, with the following properties:

   φ(X) < 0 for X ∈ M , φ(X) > 0 for X / ∈ M , φ(X) = 0 for X ∈ ∂M = S , (1)
where M denotes the closure of the set M. φ is commonly referred to as the level set function. With such a choice of φ, we no longer have to worry about the scene surface. We should pay more attention to the domain of φ, and the related calculations are independent of the location and topology of S.

Implicit Ray Tracing

We now set up the foundation of our algorithm and derive properties of ray tracing of a single viewpoint in an implicit framework.

Let V 0 denote the position of the viewpoint in Ω, and any other point in Ω is represented by X i (i = 1, 2, 3, . . .). Then, a point X 1 is considered invisible to the viewpoint V 0 if it lies inside the scene surface or if it is occluded by another point X 2 in the radial direction ---→ X 1 V 0 (Fig. 3). We observe further that the visibility status of points, which share the same radial direction, satisfy a causality condition: if a point is occluded, then all other points behind it in the same radial direction are also occluded. This means that a point X 1 is considered invisible to the viewpoint V 0 if and only if there is at least one point in the line segment ---→ X 1 V 0 that lies inside the scene surface. Under this condition, the idea of out algorithm is to start from each point in Ω, shoot a ray towards the viewpoint, and check whether there is a point being inside the scene surface in this ray.

Sometimes, with an explicit representation of the surface it could be difficult to determine whether a point is inside or outside the surface. In contrast, implicit surface representations include some very powerful geometric tools. We can determine which side of the surface a point is on simply by looking at the local sign of the level set function φ. That is, X is inside the surface when φ(X) < 0, 3. Ray Tracing. A point is considered invisible to the viewpoint only in two cases. The first case is that it lies inside the scene surface and therefore it is occluded by the scene surface. The second case is that although it lies outside the scene surface, it is still occluded by another point before it in the same radial direction. X ′ 1 is invisible to V0 in the first case and X1 is invisible to V0 in the second case.

0 V 1 X 2 X S 1 ¢ X Fig.
outside the surface when φ(X) > 0, and on the surface when φ(X) = 0. Thus, if we use another function ψ(X) to implicitly represent the visibility status of point X, the problem becomes to compute

ψ(X) = min ξ∈ --→ xv0 (φ(ξ)), ( 2
)
where ξ is a point on ---→ XV 0 . If ψ(X) is negative, then X is occluded. As a matter of fact, our visibility algorithm is an approximation of this formula.

Numerical Implementation

In order to implement Eq. ( 2) numerically, a uniform Cartesian grid defined as

{(x i , y j , z k )|1 ≤ i ≤ m, 1 ≤ j ≤ n, 1 ≤ k ≤ l} in 3D
space is used to discretize the closed volume domain Ω, and all of the functions and quantities are defined on this grid. We will use the conventional indexing scheme to enumerate the grid nodes in Ω; i.e., X i,j,k = X ll + ∆x(i, j, k), where X ll is the lower left corner of Ω, and ∆x is the mesh size. With this notation, we define a voxel to be the cube with vertices {X i+p,j+q,k+r : p, q, r ∈ {0, 1}} ⊂ Ω for some (i, j, k). We need to update Eq. ( 2) in order to take advantage of this Cartesian grid.

At each grid point X, we first find the voxel containing X. From this voxel, find the face F that intersects the ray segment ---→ XV 0 . If ---→ XV 0 intersects F at a point, we denote this intersection by X ′ (Fig. 4

left). If

---→ XV 0 lies in F, this problem is reduced to a 2D problem, and X ′ is the intersection of ---→ XV 0 and an edge of F (Fig. 4 right). In fact, X ′ can be regarded as the closest point to X in ---→ XV 0 within the accuracy of the Cartesian grid. Hence, we use a recursion to rewrite Eq. (2) as

ψ(X) = min(ψ(X ′ ), φ(X)). (3)
Note that if ψ(X ′ ) < 0 or φ(X) < 0, then ψ(X) will be less than zero in the end. It means that X is invisible if X ′ is invisible or if X lies inside the scene surface. Therefore, Eq. ( 3) is essentially the causality condition of visibility we mentioned above.

In most cases, X ′ does not lie on the grid node, and we need to interpolate the value of ψ from the grid points closest to X ′ . For example, we use ψ(A i ), i = 1, 2, for linear interpolation of ψ(X ′ ) in 2D case (Fig. 4 right) and use ψ(A i ), i = 1, 2, 3, 4, for bilinear interpolation of ψ(X ′ ) in 3D case (Fig. 4 left). So, the values of ψ on all the neighboring grid points around X ′ must have been computed before estimating the value of ψ(X). A key part of this visibility algorithm is to organize all the grid points in Ω in a special sequence to make this interpolation procedure feasible. We describe such a in detail below, and use a queue data structure denoted by Q to store this sequence for convenience. Suppose the grid in Ω is bounded by 0 and nx in the X direction, 0 and ny in the Y direction and 0 and nz in the Z direction. We first find the cube voxel C 0 in Ω containing the viewpoint V 0 , i.e., V 0 ∈ C 0 := [x i0 , x i0+1 ) × [y j0 , y j0+1 ) × [z k0 , z k0+1 ). We can therefor start from this voxel C 0 and store the grid points into Q following the ray directions outwards. In order to represent our idea simply, we use three parameters denoted by s1, s2 and s3 to describe multiple directions from C 0 for i, j and k index respectively. The relationship between s1 and i is as follows:

X 0 V ¢ X 1 A 2 A 3 A 4 A X 0 V 1 A 2 A ¢ X F F
   i = i 0 -1 to 0 when s1 = -1, i = i 0 to i 0 + 1 when s1 = 0, i = i 0 + 2 to nx when s1 = 1,
and the relationship between s2 and j, s3 and k is similar to this. In the first step, we put the grid points X i,j,k where (i, j, k) satisfy (s1, s2, s3)=(0, 0, 0) into Q. Actually, X i,j,k are the vertices of C 0 (Fig. 5a red points) in this case. In the second step, put the grid points X i,j,k where (i, j, k) satisfy (s1, s2, s3)=(0, 0, -1), (0, 0, 1), (0, -1, 0), (0, 1, 0), (-1, 0, 0) or (1, 0, 0) into Q (Fig. 5b blue points). In the third step, put the grid points X i,j,k where (i, j, k) satisfy (s1, s2, s3) =(0, -1, -1), (0, -1, 1), (0, 1, -1), (0, 1, 1), (-1, 0, -1), (-1, 0, 1), (1, 0, -1), (1, 0, 1), (-1, -1, 0), (-1, 1, 0), (1, -1, 0) or (1, 1, 0) into Q (Fig. 5c green points). Finally, put the remaining grid points X i,j,k where (i, j, k) satisfy (s1, s2, s3)=(-1, -1, -1), (-1, -1, 1), (-1, 1, -1), (-1, 1, 1), (1, -1, -1), (1, -1, 1), (1, 1, -1) or (1, 1, 1) into Q (Fig. 5c It is easy to see that each grid node in such a queue Q can be visited in a specified order that maintains the causality. Thus our algorithm reads:

for each X in the queue Q if X is the vertices of C 0 Set ψ(X) = φ(X); else Find the face F and X ′ ∈ F ; Update ψ(X) = min(ψ(X ′ ), φ(X)); end end
In this algorithm, it is unnecessary for us to restrict the viewpoint V 0 to the grid node of a Cartesian grid. So it can be used to solve the visibility problem easily in multiple views.

Variational Stereo with State-Based Visibility

Most of the multi-view reconstruction methods relying on state-based visibility need a good initialization of the scene surface [START_REF] Vogiatzis | Mult-view stereo via volumetric graph-cuts[END_REF][START_REF] Faugeras | Variational principles, surface evolution, PDE's, level set methods, and the stereo problem[END_REF][START_REF] Sinha | Multi-view reconstruction using photo-consistency and exact silhouette constraints: a maximum-flow formulation[END_REF]. If the current scene state is far from the true state, state-based visibility will approximate the true visibility with significant errors. Therefore, this type of global visibility algorithm is not thought to be a good choice in the situation when no good initialization is given [START_REF] Lempitsky | Oriented visibility for multiview reconstruction[END_REF]. In this section, we formulate a variational model for multi-view reconstruction based on our proposed visibility estimate. Compared with other methods, the main benefit of our model is its independence from the initialization. We show that our state-based visibility algorithm is still a good choice when no good initialization is given.

Energetic Formulation

The goal of our idea is to render a multi-view reconstruction as an energy minimization problem by assigning an energy cost based on visibility estimate. Since the surface of the unknown scene is supposed to be nearly Lambertian, patches belonging to the true surface should have similar colors in the viewpoints observing them. Therefore, the energy cost of a patch can be defined as

Φ(X, N) = 1 T m i,j=1 Φ ij (X, N), (4)
where T is the number of items in the summation, m is the number of visible views to the current patch and Φ ij between two visible views i and j is based on the normalized cross correlation (NCC) and taken as [START_REF] Faugeras | Variational principles, surface evolution, PDE's, level set methods, and the stereo problem[END_REF] Φ

ij (X, N) = 1 -N CC(I i , I j ). (5)
Its value range between 0 (best correlation) and +2 (worst). The term N CC(I i , I j ) is the normalized cross correlation of the projections of the patch in views i and j. We direct readers to [START_REF] Alejandro | Multi-view multi-exposure stereo[END_REF] more on how to compute this term. Summation in Eq. ( 4) is only computed for those patches which are visible to the two concerned views. Thus, estimating Φ requires the step of computing the visibility status of the patches for all views. With our visibility algorithm, we rewrite Eq. ( 4) to involve the visibility status as follows:

Φ(X, N) = 1 T n i,j=1 χ(i, X)χ(j, X)Φ ij (X, N), (6)
where n is the number of all views, χ(v, X) is the Heaviside function of ψ(X) to the view v: H(ψ(v, X)). We have

χ(v, X) = H(ψ(v, X)) = 1 if ψ(v, X) ≥ 0, i.e, X is visible to v, 0 if ψ(v, X) < 0, i.e, X is invisible to v. (7)
The overall surface energy cost for the scene configuration is then calculated by integrating patches' costs over the scene surface S:

E(S) = S Φ(X, N)dA. (8)
Once Φ been chosen, a surface S which minimizes Eq. ( 8) can be regarded as the final scene surface (reconstruction result). Faugeras et al. [START_REF] Faugeras | Variational principles, surface evolution, PDE's, level set methods, and the stereo problem[END_REF] use a gradient descent evolution to solve this problem. However, the integral in Eq. ( 8) is only calculated on the surface S, and therefore the gradient descent evolution is driven in terms of this surface. For this reason, the final result in [START_REF] Faugeras | Variational principles, surface evolution, PDE's, level set methods, and the stereo problem[END_REF] tends to be sensitive to the initial position of the surface S. Since the surface to be reconstructed is represented by a level set function φ, the energy functional in Eq. ( 8) can be rewritten as

E(φ) = Ω Φ(X, N)δ(φ)|∇φ|dXdY dZ, (9)
where δ(φ) is the univariate Dirac function. Compared with Eq. ( 8), Eq. ( 9) is directly formulated in the level set domain, and therefore the minimization of it leads naturally to the evolution of φ. By introducing the term δ(φ)|∇φ| [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] in the integrand, the integral in Eq. ( 9) can be done on the whole 3D space Ω. So, the final result of φ that minimizes this energy functional is not sensitive to where it has started its evolution. We finally write down the full energy guiding the reconstruction:

E(φ) =α Ω 1 2 (|∇φ| -1) 2 dXdY dZ + β Ω Φ(X, N)δ(φ)|∇φ|dXdY dZ + γ Ω δ(φ)|∇φ|dXdY dZ, (10)
where the first term is a penalizing term [START_REF] Li | Level set evolution without reinitialization: a new variational formulation[END_REF] which force the level set function to be close to a signed distance function as a numerical remedy for maintaining stable evolution. The third term is a smoothness term that can be used as a shape prior to improve the final reconstruction result. All the terms in the functional are weighted by constants α, β, γ.

Energy Minimization

The Euler-Lagrange equation for φ of the functional in Eq. ( 10) can be given by

∂E ∂φ = -α(∆φ -κ) -βδ(φ)(∇Φ • N + Φ • κ) -γδ(φ)κ. ( 11
)
where ∆ is the Laplacian operator. κ and N are the mean curvature and unit normal of the surface respectively, which can be calculated by

κ = div( ∇φ |∇φ| ), N = ∇φ |∇φ| . (12)
Employing an artificial time variable t > 0, the steepest descent process for minimization of the functional E(φ) is the following gradient flow:

∂φ ∂t = α(∆φ -κ) + βδ(φ)(∇Φ • N + Φ • κ) + γδ(φ)κ. (13)
In practice, the Dirac function δ(φ) in Eq. ( 13) is replaced by the function δ ε (φ) defined as

δ ε (φ) = 1 π ε ε 2 + φ 2 , (14)
with a small ε > 0. It gives a globally positive approximation to δ(φ). We use δ ε (φ) in Eq. ( 14) instead of the one in [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Li | Level set evolution without reinitialization: a new variational formulation[END_REF] because it can act on all level sets of φ, while the one in [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Li | Level set evolution without reinitialization: a new variational formulation[END_REF] can only act on a few level sets around φ = 0. In this way, Eq. ( 13) has a tendency to obtain a global minimum of E(φ) with δ ε (φ) independently to the initial value of φ. The final evolution equation of φ is hence given by

∂φ ∂t = α(∆φ -κ) + βδ ε (φ)(∇Φ • N + Φ • κ) + γδ ε (φ)κ. ( 15
)
The reconstruction result can be obtained from the zero level set of φ when this evolution is over.

Experimental Results

In this section, we present the experimental results of our method on both synthetic and real datasets. In Fig. 6, we validate our visibility algorithm in some synthetic scenes that are composed of several spheres and cuboids (Fig. 6a left and Fig. 6b left). In this situation, since the true scene geometry is known, it can be very easy to determine the real visibility information of the elements with respect to several given viewpoints through our algorithm. In order to show our method usable in any multiple views, we take two cameras with known position as the viewpoints in our experiments. The visibility results (Fig. 6a right and Fig. 6b right) demonstrate the performance of our approach. In Fig. 7, we apply our visibility algorithm to a number of reconstructions with the Middlebury datasets [14]. The data consists of two objects, a dinosaur and a temple (Fig. 7a), and three different sets of input images for each one, with viewpoints forming a sparse ring, a full ring, and a full hemisphere around the object. We perform our experiments on the dinoSparseRing (16 images) and the templeSparseRing (16 images) inputs. Eq. ( 15) is used to estimate the final 3D models through an evolution process and the initial zero level set is taken as a sphere. Fig. 7b illustrates some of those initial zero level sets (shown in blue) along with the reconstructed 3D models (shown in red). It is clear that our method can successfully reconstruct the surfaces for any initial level set function.

Conclusions and Future Work

In this paper, we have presented a state-based approach to visibility estimation. It is based on an implicit ray tracing and can be used to solve the visibility problem of multiple views efficiently. With such a global visibility estimate, we have also proposed a variational formulation for multi-view reconstruction. Our main advantage over other reconstruction methods relying on state-based visibility is the independence from the initialization due to the ability of our optimization to yield a global minimum of the energy functional. Future work will consist of improving the efficiency of the photo-consistency measure as well as developing better shape priors. Quantitative experiments are also needed to evaluate the accuracy and completeness of the reconstruction results. We also plan to incorporate strategies to apply our method to non-closed surfaces.

Fig. 2 .

 2 Fig. 2. The fallibility of oriented visibility in a scene with complex geometry.

Fig. 4 .

 4 Fig. 4. An interpolation in 2D and 3D space.

Fig. 5 .

 5 Fig. 5. A special sequence of the grid nodes

Fig. 6 .

 6 Fig. 6. Visibility results of two cameras applied to the synthetic scene with known geometry. The red surface represents the scene (occluders). The green and blue surface represent the boundary between visible and invisible regions to camera one and camera two respectively. (a) The visibility results of two cameras in a simple sphere scene. (b)The visibility results of two cameras in a scene with more complex geometry.

Fig. 7 .

 7 Fig. 7. The reconstruction results. (a) Some views of the dinoSparseRing and tem-pleSparseRing datasets. (b) The reconstruction results of dino and temple (shown in red) with several different initializations (shown in blue).

Acknowledgements

This work is partially supported by the Natural Science Foundation of China (60675021) and the Chinese High-Tech Program (2006AA01Z120).


Stephen Gould

Paul Baumstarck

Morgan Quigley

Andrew Y Ng

Daphne Koller

Integrating Visual and Range Data for Robotic Object Detection

The problem of object detection and recognition is a notoriously difficult one, and one that has been the focus of much work in the computer vision and robotics communities. Most work has concentrated on systems that operate purely on visual inputs (i.e., images) and largely ignores other sensor modalities. However, despite the great progress made down this track, the goal of high accuracy object detection for robotic platforms in cluttered real-world environments remains elusive. Instead of relying on information from the image alone, we present a method that exploits the multiple sensor modalities available on a robotic platform. In particular, our method augments a 2-d object detector with 3-d information from a depth sensor to produce a "multi-modal object detector." We demonstrate our method on a working robotic system and evaluate its performance on a number of common household/office objects.

Introduction

Endowing autonomous robots with the ability to perceive objects in their environment is a notoriously difficult problem in computer vision. One standard solution is to equip the robot with a video camera and run state-of-the-art, image-based object detection algorithms on each frame. However, there are a number of difficulties with this approach.

Simply finding each object-let alone classifying it-is a non-trivial task in cluttered real-world environments. Many state-of-the-art object detectors combine the process of object detection (localization) and object recognition (e.g., the slidingwindow approach), and only compute features local to the region being considered. While this approach works well for uncluttered scenes with well-framed objects, it fails when applied to real-world images as seen by an autonomous robot. 1An autonomous robot needs to operate in a world of small, occluded objects that are often ill-framed. Many visual percepts will not have been encountered before, and portions of the scene may contain confusing or conflicting visual cues when taken out of context. It has been shown that contextual information greatly improves the performance of vision-based object detectors especially when the objects of interest are small [START_REF] Fink | Mutual boosting for contextual inference[END_REF][START_REF] Murphy | Using the forest to see the tree: a graphical model relating features, objects and the scenes[END_REF][START_REF] Hoiem | Putting objects in perspective[END_REF]. However, even when this information is taken into account, the context extracted from image features is inherently 2-dimensional and can be unreliable in natural scenes with significant clutter.

Sensor specifications:

high resolution video camera (704×480 images at 15 frames per second);

laser range scanner (75 vertical scanlines per second at 0.5 • resolution);

laser pan/tilt unit (37.5 • /s horizontal pan);

odometry from navigation laser and wheel encoders. Exacerbating the above difficulties are the real-time processing constraints imposed by robotic applications. This, in general, requires that the vision system be simple, and, in particular, rules out the option of going to higher resolution unless a sophisticated attentional mechanism is employed.

The shortcomings of image-based object detection can be addressed on a robotic platform by exploiting the multiple sensor modalities available to the robot rather than relying on a single monocular or stereo-vision sensor mode. Since the robotic system is situated in its environment it can directly measure features that would help detection instead of inferring these quantities from 2-d image cues. For example, measuring temperature can help in the detection of people, and finding large supporting surfaces can help locate small objects. With the availability of cheaper and easier-to-use sensors, such as infrared cameras and lasers, it makes sense to leverage these different sensor modalities for vision.

In this paper, we propose to augment state-of-the-art robotic visual systems with input from a laser depth scanner and to combine the 2-d and 3-d sensor modalities to enhance object detection in cluttered real-world environments. We fuse the depth and visual data so that, for each pixel in the robot's field-of-view, we have color/intensity, depth (location in 3-d space), and surface normal information. This allows us to directly extract contextual cues (such as height above the ground) and 3-d features (such as object size). Figure 2 depicts these features for a simple office scene.

Our method consists of a number of stages which can easily be pipelined for real-time processing on a robotic platform. (See Figure 1 for a schematic of our dataflow). First, using super-resolution techniques, we combine frames from a video camera with a low resolution, or sparse, point cloud generated from a laser scanner to produce a high resolution, or dense, point cloud in the video camera's frame of reference. Next, we apply a sliding-window object detector at multiple image scales: for each location ((x, y)-position and scale) in the image, we compute local 2-d patch features and corresponding 3-d contextual features. The patch features are provided as input to a trained boosted classifier which, like standard image-only object detectors, estimates the probability of an object being at the location given just the image features. Finally, we apply a learned logistic classifier for each object to a feature vector comprising this image-only probability and the computed 3-d features. A distinct advantage of this approach is that the image-only classifier can be trained using standard state-of-the-art techniques (and on large datasets) while our "multi-modal" classifier still takes advantage of 3-d contextual information. We demonstrate our method on a number of small household/office objects in natural environments. Importantly, we show that by fusing information from the depth sensor we are able to significantly improve object detection rates over state-of-the-art image-only techniques in cluttered real-world environments.

Background and related work

Sliding-window object detection is a popular technique for identifying and localizing objects in an image. The method has been very successfully applied to face detection [START_REF] Viola | Robust real-time face detection[END_REF] and can be highly efficient when combined with a cascade of boosted ensembles (CoBE) classifier and simple Haar-wavelet features [START_REF] Brubaker | On the design of cascades of boosted ensembles for face detection[END_REF]. However, Haarwavelet-like features tend to perform poorly on other object classes, and researchers have developed other, more sophisticated discriminative features (e.g., histogram of oriented gradients [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], biologically inspired (visual cortex) features [START_REF] Serre | Object recognition with features inspired by visual cortex[END_REF], or patchbased features [START_REF] Torralba | Contextual models for object detection using boosted random fields[END_REF]) at the cost of processing speed. When applied to multiple object classes, features can be shared amongst the detectors, thereby amortizing the cost of feature extraction [START_REF] Torralba | Sharing visual features for multiclass and multiview object detection[END_REF].

A number of researchers have shown that context can significantly improve object detection accuracy, especially when the objects are small [START_REF] Fink | Mutual boosting for contextual inference[END_REF][START_REF] Murphy | Using the forest to see the tree: a graphical model relating features, objects and the scenes[END_REF][START_REF] Torralba | Contextual models for object detection using boosted random fields[END_REF][START_REF] Torralba | Contextual guidance of eye movements and attention in real-world scenes: the role of global features in object search[END_REF][START_REF] Hoiem | Putting objects in perspective[END_REF]. However, little work has focused on 3-d context. Notable exceptions are the innovative works of Hoiem et al. [START_REF] Hoiem | Putting objects in perspective[END_REF] and Leibe et al. [START_REF] Leibe | Dynamic 3d scene analysis from a moving vehicle[END_REF] who infer the camera location and scene geometry from a single 2-d image or stereo video stream, respectively. These works reconstruct rough geometry of street scenes (pedestrians and cars) and cannot, for example, be used for estimating 3-d features of small objects. Instead of trying to infer 3-d structure, we propose to measure it directly using a laser. Some novel works use 3-d point clouds from depth sensors for detecting objects. Nuchter et al. [START_REF] Nuchter | Accurate object localization in 3d laser range scans[END_REF] use separate reflectance and depth information from a single sensor to detect known objects. Other works focus on discovering geometric primitives [START_REF] Rabbani | Efficient hough transform for automatic detection of cylinders in point clouds[END_REF], or detecting large novel objects but without recognition [START_REF] Chen | Object detection for a mobile robot using mixed reality[END_REF]. Unlike these works, we detect and recognize small objects from sparse depth information in real-time.

Sensor fusion and scene representation

In this section, we describe how we process the raw data from the robot's sensors (Figure 1) into a representation that unifies the modalities of the camera and the sparse depth sensor. In particular, we describe how the sparse laser measurements are used to infer depth at the (much higher) spatial resolution of the camera.

Ideally we would like to measure the depth (and surface normal) at every pixel in the image so that we can make use of 3-d geometric cues in our object detector. Unfortunately, the laser range scanner does not support the same resolution as the video camera. 2 Thus we need to resort to super-resolution techniques in order to infer the depth at every pixel. The resulting high resolution depthmap can then be used to estimate surface normals. State-of-the-art methods for super-resolution include MAP inference on a pairwise Markov random field (MRF) [START_REF] Diebel | An application of markov random fields to range sensing[END_REF][START_REF] Torres-Mendez | Reconstruction of 3d models from intensity images and partial depth[END_REF] and bilateral filtering [START_REF] Yang | Spatial-depth super resolution for range images[END_REF], and are based on the intuition that depth discontinuities usually coincide with color discontinuities in the image.

Our method is similar to the MRF model of Diebel and Thrun [START_REF] Diebel | An application of markov random fields to range sensing[END_REF]. 3 However, instead of encoding a preference for fronto-parallel planes (implicit in their formulation) we allow for arbitrarily sloped planar surfaces. Thus our method can be thought of as reconstructing a first-order approximation to each surface rather than a zeroth-order one. We also use a robust Huber penalty instead of the more commonly used ℓ 2 penalty. 4 A quantitative comparison measuring mean-square reconstruction error on a hold-out set of points showed that, on average, our model performed better than that of Diebel and Thrun [START_REF] Diebel | An application of markov random fields to range sensing[END_REF] on our office scenes-details omitted due to space constraints.

Concretely, let the image pixel intensities be {x i,j | (i, j) ∈ I}, the laser depth measurements be {z i,j | (i, j) ∈ L} and the reconstructed/inferred depth for every pixel be {y i,j | (i, j) ∈ I} where I indexes the image pixels and L ⊆ I indexes the laser measurements (projected onto the image plane). Two MRF potential functions are defined-the first penalizes discrepancy between measured and reconstructed depths, while the second encodes a preference for smoothness:

Φ ij (y, z) = h(y i,j -z i,j ; λ) (1)
Ψ ij (x, y) = w v ij h(2y i,j -y i,j-1 -y i,j+1 ; λ) + w h ij h(2y i,j -y i-1,j -y i+1,j ; λ) (2)
where h(x; λ) is the Huber penalty function, and

w v ij = exp{-c x i,j-1 -x i,j+1 2 } and w h ij = exp{-c x i-1,j -x i+1,j 2
} are weighting factors indicating how unwilling we are to allow smoothing to occur across vertical and horizontal edges in the image as in [START_REF] Diebel | An application of markov random fields to range sensing[END_REF]. We can now define our super-resolution MRF as

p(y | x, z) = 1 η(x, z) exp    -k (i,j)∈L Φ ij - (i,j)∈I Ψ ij    (3)
where k specifies the trade-off between measurement reconstruction and smoothness, and η(x, z) is the normalization constant (partition function). We initialize our depth estimates y by first projecting measured depths into the camera plane and performing quadratic interpolation to obtain an initial depth estimate at every pixel. We then find the MAP solution by minimizing k (i,j)∈L Φ ij + (i,j)∈I Ψ ij using the L-BFGS algorithm [START_REF] Liu | On the limited memory method for large scale optimization[END_REF]. The process is illustrated in Figure 3. Optimization is fairly quick, taking less than 10 seconds on a 320×240 image. Most of the progress towards the optimal solution is made in the first few iterations and high quality depth estimates can still be obtained by stopping early. Furthermore, in a robotic application where data is continuously being streamed, the algorithm can be initialized from the previous frame (with motion compensation) resulting in much faster convergence.

Once we have the high resolution depthmap we reconstruct the location of each pixel in 3-d space (relative to the camera's reference frame) { Xi,j ∈ R 3 | (i, j) ∈ I} by projecting a ray through the camera plane at each pixel location (i, j) and scaling by the inferred depth. We also infer point normals, {n i,j ∈ R 3 | (i, j) ∈ I, ni,j = 1}, from the local neighborhood around each point. 5 Figure 4 shows sample results from our super-resolution procedure. The tuple (x i,j , Xi,j , ni,j ) ∈ R 7 forms the basis of our 3-d scene representation.

Multi-modal object detection

In this section, we describe how we use the unified scene representation of Section 3 and a sliding-window approach to perform multi-modal object detection. Specifically, we discuss 2-d and 3-d feature extraction, and describe our method for combining them into a single probabilistic model.

Image-only object detectors

Sliding-window object detection is a simple yet effective approach to simultaneous localization and recognition of objects in images. The approach involves scanning the image with a fixed-size rectangular window and applying a classifier to the subimage defined by the window. The classifier extracts image features from within the window (sub-image) and returns the probability that the window (tightly) bounds a particular object. The process is repeated on successively scaled copies of the image so that objects can be detected at any size. The (x, y)-location of the window and scale σ of the image implicitly defines the bounding box, or set of pixels, belonging to the candidate object B(x, y, σ) = {(i, j)} ⊆ I.

Our image features are similar to Torralba et al. [START_REF] Torralba | Sharing visual features for multiclass and multiview object detection[END_REF] (and other works by the same author). Here a dictionary of image patches is constructed at training time by randomly selecting small patches from a set of training images. Associated with each patch is a spatial mask over which the patch response is valid. The mask is derived by taking a small rectangular region around the patch's original (image) location (see Figure 5(d)). In addition to patches taken from the (intensity) image, we also extract patches from the edge-filtered (gradient) image. 6 Thus our dictionary is comprised of the triplet of patch, spatial mask, and image operator (intensity or gradient). Figure 5 shows example training images and a patch dictionary.

We scan the image at multiple locations and scales. The response for feature f at location (x, y) and scale σ is given by v f (x, y, σ) = max w f |T f (I σ ) ⊗ g f | where I σ is the image at scale σ, T f (•) is the image transform associated with feature f (intensity or gradient), g f is the image patch, w f is the spatial mask, and ⊗ is the normalized cross-correlation operator.

Given the dictionary of patch features and a set of training images we learn a gentle-boost classifier [START_REF] Friedman | Additive logistic regression: a statistical view of boosting[END_REF] over two-split decision stumps for each object class. At run time, we apply the boosted classifier to each candidate location (x, y) and scale σ to obtain the probability (based on image features alone) of that location containing the given object, P image (o | x, y, σ). We then use the log-odds ratio

f 2d (x, y, σ) = log P image (o | x, y, σ) 1 -P image (o | x, y, σ) (4)
as a feature in our multi-modal object detector. Note that these detectors do not make explicit use of any 3-d information and therefore can be trained and evaluated using standard state-of-the-art techniques. Briefly, we use 100-200 positive and 20,000 negative training examples for each object class. The positive examples were downloaded from the web or manually collected using a digital camera. The negative examples were collected by randomly snipping rectangles from a five-minute video sequence. All examples were scaled to 32 pixels for the smaller dimension. (See Figure 5.)

We construct our patch dictionary by randomly selecting 10 (intensity and gradient) patches from each positive training sample. The patches vary in size from 4 × 4 to 16 × 16 pixels, and we fix the spatial mask w f to 7 × 7. We then train the object detectors in two stages for improved efficiency: first we select 2,000 6 We convolve the intensity image with 3 × 3 Sobel kernels to obtain horizontal and vertical gradient images G h (x, y) and Gv(x, y). The edge-filtered image is then given by the gradient magnitude negative training examples at random and train a boosted classifier for 50 rounds; next, we trim our patch dictionary to remove all patches not used by this classifier; finally, we retrain using all 20,000 negative training examples to obtain our final image-only object detectors. The resulting patch dictionary typically contains 50-75 entries compared to the 500 entries used in [START_REF] Torralba | Sharing visual features for multiclass and multiview object detection[END_REF].

I edge (x, y) = p G h (x, y) 2 + Gv(x,

3-d features

For each candidate location (x, y) and image scale σ, we compute 3-d features by taking the bounding rectangle implicitly defined by the location and scale B(x, y, σ) = {(i, j)} ⊆ I and conceptually projecting that region into the scene as shown in Figure 6. We shrink the set of points enclosed by the bounding box by removing the (i, j) corresponding to points Xi,j that lie in the outer 5% (in either the x-, y-, or z-directions) of the points in B. This removes outliers and most background points. The resulting shrunken set of points B ′ is then used for computing the features. Let { Xi,j ∈ R 3 | (i, j) ∈ I} be the location of each pixel in 3-d space (relative to the camera's frame of reference), and let {n i,j ∈ R 3 | (i, j) ∈ I, ni,j = 1} be the estimated surface normal vector for each pixel. The centroid and covariance of the points in B ′ are

µ X = 1 |B ′ | (i,j)∈B ′ Xi,j and Σ X = 1 |B ′ | (i,j)∈B ′ ( Xi,j -µ X )( Xi,j -µ X ) T (5)
where |B ′ | is the number of points in B ′ . A similar computation gives µ n and Σ n , the mean and covariance over surface normal vectors. We now enumerate the object attributes and contextual cues captured by our 3-d features and provide some insight into why they are useful for improving object detection. Implementation details for how we capture each attribute/contextual cue as a vector-value feature are also described.

Height above ground. Many objects are located at consistent heights above the floor. For example, computer monitors are often found on desks which are a standard height above the floor; door handles are placed at a relatively fixed location on the door where most people find it comfortable to reach; and wall clocks are placed high enough so that they can be seen above other objects in the room. Since we calibrate the point cloud xz-plane to coincide with the floor, our height above ground feature is simply

f height-above-ground = µ y X (µ y X ) 2 T (6)
where µ y X indicates the y-component of µ X . Distance from robot. Objects farther from the robot are harder to see, not only because of their diminished size but also because of lighting effects and depth measurement inaccuracy. Although the probability of an object does not depend its distance from the robot, knowing this distance allows our model to compensate for the influence the above effects have on other features (e.g., the 2-d detectors):

f distance-from-robot = µ z X (µ z X ) 2 T (7)
where µ z X indicates the z-component of µ X . Surface variation and orientation. Most objects have significant surface variation and extend above their supporting surface. Obvious exceptions are flat objects, such as LCD monitors, which appear as a single plane in the 3-d point cloud and so are distinguished by their lack of surface variation. We compute the following vector-valued features to capture these attributes:

f surf-var =    Σ xx X + Σ zz X Σ yy X Σ xx n + Σ zz n Σ yy n    , f surf-orientation =   µ y n 1 -(µ y n ) 2 (µ y n ) 2   (8)
where Σ X and Σ n are the covariance matrices defined above. Object dimensions. One of the most defining attributes of an object is its size. Knowing this information allows a significant number of false candidates to be rejected. The width, height, and (projected) area of the object can be estimated by considering the projection of the bounding rectangle into 3-d space:

f width = max (i,j)∈B ′ Xx i,j -min (i,j)∈B ′ Xx i,j (9)
f height = max (i,j)∈B ′ Xy i,j -min (i,j)∈B ′ Xy i,j (10)
f area = f width × f height (11)
where Xx i,j and Xy i,j denote the x-and y-components of Xi,j , respectively. Our feature vector also includes the square of these terms. A measure of object depth is already captured by our surface variation feature and so is not repeated here.

The above features are assembled into a single 17-dimensional descriptor f 3d which captures the 3-d attributes of a possible object.

Multi-modal object detectors

Our multi-modal object detectors are simple binary logistic classifiers based on the 2-d and 3-d features, f 2d and f 3d , defined above. Note that, although this is a simple representation, the inclusion of squared terms for many of the features allows us to learn rich decision boundaries.

The probability for an object o existing at location (x, y) and scale σ in the image plane can be written as

P (o | x, y, σ) = q θ T 3d f 3d + θ T 2d f 2d + θ bias (12)
where q(s) = 1 1+e -s is the logistic function and θ = {θ 3d , θ 2d , θ bias } are the learned parameters of the model. The parameters θ 2d and θ 3d trade-off between 2d and 3d features while the bias term θ bias models the prior prevalence of the object.

We learn one multi-modal model per object class on data (images and corresponding point clouds) collected from static scenes. We do this rather than using video sequences to avoid the introduction of bias during training and evaluation due to the high correlation between consecutive video frames. The images are annotated with a bounding box and class label for each object of interest. We construct a training set by first running the image-only object detectors over each image, keeping all detections with a probability above 0.001. Any detection overlapping by more than 50% with a groundtruth annotation is used as a positive example while all other detections are used as negative examples. If for any groundtruth annotation there was no overlapping detection, we further run the image-only object detector on the bounding box for that annotation and add it to our positive examples. Finally, for each training example, we extract the 3-d features as detailed above and learn the parameters θ of our logistic classifier using Netwon's method so as to maximize the log-likelihood of the training set. We use ℓ 2 regularization to prevent over-fitting, with weight chosen by cross-validation on the training set.

Experimental results

We collected 420 static images and corresponding sparse point clouds of cluttered scenes. The scenes contained a number of small objects (coffee mugs, disposable cups, monitors, wall clocks, door handles, and ski boots) which we would like to detect, as well as distractors (see Figure 7). A number of the scenes were extremely challenging and the authors even had trouble identifying objects in some images during groundtruth labeling because of image resolution. We performed kfold cross-validation and report the aggregate performance over the hold-out sets. On each fold we learned the model parameters and regularization weight using the training set. The data in the hold-out set was only used for testing.

We evaluate our performance by comparing against image-only object detection. Figure 8 shows precision-recall curves for our learned object detectors. The multi-modal detectors (solid red curve) are consistently superior to the image-only detectors (dashed blue curve).

In general, 3-d features significantly help, especially when strong features such as size and location overcome large intra-class appearance variation (door handles and ski boots) or lack of discriminating visual features (computer monitors). The cup class performs badly for both detectors primarily due to its lack of distinguishing features at small scale. We found a large number of false-positive cups coming from other small objects such as coffee mugs. This suggests that a model which considers the location of other objects may improve accuracy for such classes.

In order to understand the contribution that each of our features makes to improving object detection, we evaluated the performance of the image-only detector augmented with each 3-d feature separately (see Table 9). Here we compare performance by measuring the maximum F 1 -score. As expected, object dimensions and height above the ground are the strongest individual features.

Finally, to gain an intuition for where our method works, we provide some representative results in Figure 7. The first row shows that our multiple sensors help when objects are partially occluded. Here the monitor is detected even though the left edge is not visible. With some visual features missing, the image-only detector cannot conclude the existence of the monitor. The second row shows how knowing object size can be used to reject the false-positive disposable cup despite the strong visual resemblance. In this scene our method also successfully detects one of the two door handles. The last two rows show how image-only detectors are easily confused by significant clutter or textured backgrounds. Interestingly, the multimodal detector incorrectly labels some trash as a coffee mug (bottom of the third scene). A patch on the whiteboard that resembles a door handle is also mistakenly labeled by our detector in the last scene.

Discussion

In this paper we proposed a multi-modal object detector for robots situated in their environments. We showed how robots can exploit 3-d data from a low resolution depth sensor-the most common mode of robotic perception after optical sensingby combining it with 2-d image data for robotic object detection in real-world environments. Our main contributions were two-fold. First, we showed how superresolution techniques allow us to obtain a high resolution scene representation consisting of pixel intensities, 3-d locations, and surface normals. Second, we showed how 3-d features and contextual cues derived from this scene representation can be combined with a state-of-the-art 2-d object detector to significantly improve detection accuracy. Importantly, our method works with any 2-d object detector.

Our experimental results showed that the multi-modal detector improves over a baseline 2-d detector for common household/office objects. The improvement is most striking for objects which lack distinguishing features or those with high intra-class variance. Object size and location (height above the ground) are the strongest 3-d features, both of which are easily derived from our multi-sensor scene representation. An exciting avenue for further work is the exploration of more sophisticated 3-d features while still restricting ourselves to real-time depth sensors which are inherently low resolution.

Our architecture also fits well with the needs of real-time robotics. The data processing pipeline allows for multi-threading and can be scaled up to an arbitrary number of object detectors. The slowest operation is the extraction of the patchresponse features used by our 2-d detectors (taking 10-20s per frame per object class). Performance can be improved by pruning parts of the scene based on 3-d cues (such as size or surface variation) which are quick to compute.

Vision-only object detection systems are plagued by a number of difficulties in real-world scenes, e.g., lighting, texture, occlusion, etc. Many of these difficulties can be overcome by augmenting visual perception with complementary sensor modalities (e.g., depth and infrared) using the methods described in this paper, providing a step towards robust object detection for autonomous robots.

Fig. 1 .

 1 Fig. 1. Our robotic platform, hardware specifications, and data flow: a sparse point cloud is reconstructed in software by knowing the relative pose between the pan/tilt unit and the video camera, and correcting for robot motion using odometry. The depth of every pixel is then inferred, and features from the 2-d image and the 3-d point cloud are combined by our probabilistic object detectors.

Fig. 2 .

 2 Fig. 2. Example 3-d feature responses for a window size of 40 × 40 pixels (approximately the size of the coffee mug in this scene). Light areas indicate stronger feature values. For example, the surface normal feature can be used for identifying supporting surfaces such as the table with horizontal support or the walls with vertical support.

Fig. 3 .

 3 Fig. 3. Results from our super-resolution MRF on a typical office scene (a). Shown are (b) initial interpolated depth estimates; and (c) final depth estimates at convergence (526 iterations).

Fig. 4 .

 4 Fig. 4. 3-d views of (a) the sparse point cloud from laser scanlines, (b) reconstructed (dense) point cloud, and (c) reconstructed surface normals (color represents normal direction).

y) 2 .Fig. 5 .

 25 Fig. 5. Example positive and negative training examples for the "mug" class are shown in (a). An example dictionary for the class "mug" is shown in (d). Adjacent columns represent patch g f and spatial mask w f for each dictionary entry. White in the spatial mask indicates size of patch, and gray boundary indicates valid response area.

Fig. 6 .

 6 Fig. 6. Illustration of a candidate object location in the image plane being projected into 3-d space.

  (a) Original scene (b) Super-res. depth (c) Image-only det. (d) Multi-sensor det.

Fig. 7 .

 7 Fig. 7. Representative scenes showing our results from multi-modal object detectors compared against image-only detectors. (Best viewed in color.)

Fig. 9 .

 9 Fig. 9. Comparison of maximum F1-score for image-only detectors augmented with individual 3-d features. The delta over image-only detectors is given in parentheses. The 3-d only results ( †) are provided for comparison and are calculated for the same set of candidate rectangles returned by the image-only detector (but the log-odds ratio feature is not used).

The best performing methods in recent PASCAL Visual Object Classes Challenge[START_REF] Everingham | The PAS-CAL Visual Object Classes Challenge[END_REF] (which contains images of difficult natural scenes) achieved less than 50% average precision for most object classes.

The video camera on our robot has a base resolution of 0.1 • and can optically zoom down to 0.01 • . On the other hand, modern (full scene scanning) depth sensors, such as the SwissRanger SR-3000, only have a resolution of approximately 0.5 • .

Note that in their setup the camera and laser were axis-aligned and the scene imaged by rotating through 360 • . Thus they did not have to deal with calibration issues between the two sensors nor occlusions.

The Huber penalty function, h(x; λ) = x 2 for -λ ≤ x ≤ λ and λ(2|x|λ) otherwise, is convex and thus we can still find the MAP solution exactly.

Here we compute the SVD of the 3 × 3 covariance matrix of points Xi ′ j ′ near to the point in question Xij and take the surface normal estimate nij to be the direction of the singular vector corresponding to the smallest singular value. We resolve ambiguity in the sense of the normal vector by taking the solution that points towards the camera.

Acknowledgements

Support from the Office of Naval Research under MURI N000140710747 is gratefully acknowledged.


Marcelo A Colominas
email: macolomi-nas@bioingenieria.edu.ar

Anne Humeau-Heurtier

Member, IEEE Gastón Schlotthauer
email: gschlotthauer@conicet.gov.ar

Orientation-Independent Empirical Mode Decomposition for Images Based on Unconstrained Optimization

Keywords: Empirical mode decomposition, unconstrained optimization, data-driven, non-stationary image. I

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

components plus a final trend:

X (x, y) = K k=1 α k (x, y) cos φ k (x, y) + A K (x, y), ( 1
)
where every local detail D k (x, y) = α k (x, y) cos φ k (x, y) is riding on a local mean A k = K i=k+1 D i + A K . Empirical Mode Decomposition (EMD) [START_REF] H U A N G | The mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] provides av separation of the form of Eq. ( 1) for 1D-signals. It locally separates, in a completely data-driven and unsupervised manner, signals into fast and slow oscillations. A post-processing via Hilbert Transform (HT) provides an estimation of the local frequencies and amplitudes. The global procedure of EMD plus HT is the so-called Hilbert Huang Transform (HHT). In the EMD jargon the local details are called Intrinsic Mode Functions (IMFs) or simply modes: functions symmetrically oscillating around zero, albeit with modulation both in amplitude and frequency. The final trend, usually a monotonic function, is called the residue.

The adaptivity and flexibility of EMD encouraged several researchers to adapt this method for 2D-data. Nunes et al. [START_REF] Nunes | Image analysis by bidimensional empirical mode decomposition[END_REF] developed a bidimensional EMD where they used morphological operators to detect regional maxima and radial basis function for surface interpolation. Damerval et al. [START_REF] Damerval | A fast algorithm for bidimensional EMD[END_REF] introduced a fast bidimensional EMD where a piecewise cubic interpolation of a Delaunay triangulation is used to find the envelopes. Linderhed [START_REF] Linderhed | Image empirical mode decomposition: A new tool for image processing[END_REF] introduced an image EMD where the envelopes are found through thin plate splines. Wu et al. [START_REF] Wu | The multi-dimensional ensemble empirical mode decomposition method[END_REF] proposed a multi-dimensional approach where they decomposed rows and columns independently and combined the results based on a minimal-scale combination strategy. However, all of these algorithms rely on an iterative procedure called sifting process which lacks of convergence guarantees. To avoid this iterative process, several optimization-based strategies have recently been developed for 1D-signals [START_REF] Pustelnik | Empirical mode decomposition revisited by multicomponent non-smooth convex optimization[END_REF], [START_REF] Oberlin | An alternative formulation for the empirical mode decomposition[END_REF]. In [START_REF] Colominas | An unconstrained optimization approach to empirical mode decomposition[END_REF], the authors introduced an unconstrained optimization approach to EMD which is simpler and faster than the other proposals. The main advantages of the latter method compared to others are the following: the explicit computation of envelopes to find the local mean is not needed; there is no use of explicit spline interpolations; the proposed method provides an analytical solution. Moreover, the computational cost is similar to the one of EMD and only one parameter has to be set.

Nevertheless, such an approach for bidimensional data has not been proposed yet. We therefore herein introduce an orientation-independent EMD for images that is based on an unconstrained optimization problem. We show that, compared to other 2D-EMD-based methods, our framework has the advantages of being computationally efficient, orientation-independent, and leads to better performances for the decomposition of AM-FM images.

The paper is organized as follows. In Section II, we recall the unconstrained optimization approach to EMD for 1D-signals. Section III introduces the framework that we propose for the decomposition of bidimensional data. Experimental results for synthetic and experimental data are presented and discussed in Section IV. We end with a Conclusion.

II. UNCONSTRAINED OPTIMIZATION APPROACH TO EMD

The separation of a signal in detail (mode) plus approximation (trend or local mean) can be done in several ways. Traditional EMD [START_REF] H U A N G | The mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] uses the so-called sifting process, which heavily relies on cubic spline interpolations and does not have a guaranteed convergence. On the other hand, the optimization approach to EMD introduced in [START_REF] Colominas | An unconstrained optimization approach to empirical mode decomposition[END_REF] solves the following convex optimization problem to find the local mean a ∈ R N of a given signal x ∈ R N :

(P) min a ||P x (x -a)|| 2 2 + λ||La|| 2 2 , (2)
where matrix P x ∈ R N×N is an operator, depending only on the locations of local extrema of x, which models the penalization imposed on the mode d = xa at each extrema, matrix L ∈ R N×N is a second-order difference matrix, a and x are considered column vectors, and || • || 2 stands for the ℓ 2 -norm. This problem has a unique solution:

a * = (P T x P x + λL T L) -1 P T x P x x.

(3)

The matrix P x is a key element in the problem statement. It has as many non-zero rows as local extrema of the signal x. Let t l , l = 1,...,L, with

t 1 < ••• < t l < t l+1 < ••• < t L ,
b e the locations of the local extrema of x, which act as estimates of the local extrema of the mode d. The matrix evaluates the following function in every local extrema:

p l , d =d(t l ) + d(t l+1 )(t l -t l-1 ) t l+1 -t l-1 + d(t l-1 )(t l+1 -t l ) t l+1 -t l-1 , (4)
where the vector p l ∈ R N is the t l -th row of matrix P x , and it has non-zeros elements only on the t l-1 , t l , t l+1 -th positions:

p l =[... t l+1 -t l t l+1 -t l-1 ...1 ... t l -t l-1 t l+1 -t l-1 ...]. (5)
The goal of Eq. ( 4) is to compare the extremum d(t l ) with its mirror point on the would-be other envelope which is locally defined thanks to d(t l-1 ) and d(t l+1 ) [START_REF] Pustelnik | A multicomponent proximal algorithm for empirical mode decomposition[END_REF], thus favoring the symmetry of the mode. It can be clearly seen that if t l is a local minimum (resp. maximum), then both t l-1 and t l+1 are local maximum (resp. minimum). In this approach the IMF conditions (envelope symmetry) are evaluated only on the signal local extrema and not throughout the whole time span. Equations (3) and ( 5) allow to perform the decomposition in a deflationary scheme, which is outlined in Algorithm 1 [START_REF] Colominas | An unconstrained optimization approach to empirical mode decomposition[END_REF].

Algorithm 1 UOA-EMD for 1D Data (Signals)

III. PROPOSED FRAMEWORK

A. Bidimensional Approach

Let us consider now a function of two real variables

X = X (x, y), with X ∈ R N 1 ×N 2 . The countable set of isolated local extrema is L = z l = (x l , y l )/ X (z l )>X (z) ∨ X (z l )<X (z), 0 < |z -z l | <δ l ; (6)
this definition of local extremum may vary depending on whether a strict inequality is required or not. Our desire is to find an approximation A ∈ R N 1 ×N 2 such that the penalization on every local extremum of the desired mode D = X -A is minimized. To do that, we must apply Eq. ( 4) in a 2D manner. If z l is a local maximum (resp. minimum), let us define z i , i = 1, 2, 3, as the three closest non-collinear local minima (resp. maxima) [START_REF] Schmitt | 2D Prony-Huang transform: A new tool for 2D spectral analysis[END_REF]. Then, the penalization on z l is

P l , D =D(z l ) + c 1 D(z 1 ) + c 2 D(z 2 ) + c 3 D(z 3 ), (7)
where •, • stands here for the inner product between matrices,

P l ∈ R N 1 ×N 2 is a matrix with P l (z i ) = c i , i = 1, 2, 3,
c 2 = (x l -x 1 )(y 3 -y 1 ) -(x 3 -x 1 )(y l -y 1 ) (x 2 -x 1 )(y 3 -y 1 ) -(x 3 -x 1 )(y 2 -y 1 ) (8a) c 3 = (x 2 -x 1 )(y l -y 1 ) -(x l -x 1 )(y 2 -y 1 ) (x 2 -x 1 )(y 3 -y 1 ) -(x 3 -x 1 )(y 2 -y 1 ) (8b) c 1 = 1 -c 2 -c 3 . (8c)
To favour smooth solutions, we must use a roughening matrix. We choose a discrete approximation to the Laplacian operator L i, j ∈ R N 1 ×N 2 which consists in a matrix of the same size as X with all zeros except for the kernel

K = ⎡ ⎣ 0.510 .5 1 -61 0.510 .5 ⎤ ⎦ (9)
centered at the point i, j .W h e ni, j is at the border of the matrix, we simply cropped the kernel, and modify the central value so that the sum of the kernel entries is zero. With Eq. ( 7) and the Laplacian operator, we can formulate the following optimization problem:

(P1-2D) min A l/z l ∈L P l ,(X -A) 2 + λ i, j L i, j , A 2 , ( 10
)
in which the first term models the penalization in every local extrema and the second term favours smooth solutions. However, problem P1-2D can be difficult to solve. We must modify it. The inner product between matrices can be computed as A, B =Tr(A T B) = (vec(A)) T (vec(B)),w h e r e vec(•) is the classical vectorization operation, in which the columns of a matrix are vertically concatenated to form a vector; and Tr(•) stands for the trace of a matrix. Let us define the vectors

â ∈ R N 1 N 2 = vec(A) and x ∈ R N 1 N 2 = vec(X), the matrix PX ∈ R N 1 N 2 ×N 1 N 2 which contains (vec(P l )) T as rows and L ∈ R N 1 N 2 ×N 1 N 2 which contains (vec(L i, j )) T as rows.
The following optimization problem

(P-2D) min â || PX ( x -â)|| 2 2 + λ|| L â|| 2 2 , (11)
is equivalent to (P1-2D) and has a unique minimizer at

â * = ( PT X PX + λ LT L) -1 PT X PX x, (12)
from which the desired local mean A can be constructed. Problem (P-2D) arises as a natural 2D extension of (P), and paves the road for nD extensions, simply by defining the penalization tensor for every local extrema, and a roughening tensor on nD. It also constitutes a true genuine-2D approach, as a contrast to pseudo-2D approaches, because it does not separately constraints rows, columns and diagonals. The result is an orientation-independent image processing tool.

B. Implementation Issues

The sizes of the matrices involved in problem (P-2D) might look overwhelming at first sight, and therefore the inversion in Eq. ( 12) might appear as very expensive in terms of computational cost. Sparsity must be taken into account to significantly reduce the cost. Matrices P l have at most 4 non-zero elements, and therefore PX has a proportion of 4|L|/N 2 1 N 2 2 non-zero elements, with |L| being the cardinality of the set L,a n d|L|≪N 1 N 2 . The same applies to L i, j , which has at most 9 non-zero elements, and therefore L has a proportion of 9/N 1 N 2 of non-zero elements. Efficient routines to invert sparse matrices make our proposal a very low-cost method when compared with other bidimensional EMD implementations. From Eqs. ( 7)-( 9) and ( 12), we construct a deflationary scheme to decompose an image X ,a ss h o w n in Algorithm 2. We will call our proposal UOA-EMD, for unconstrained optimization approach.

IV. EXPERIMENTS AND RESULTS

A. 2D Dirac Function

Genuine 2D approaches must be orientation-independent. In order to test this property on UOA-EMD, we present as Algorithm 2 UOA-EMD for 2D Data (Images) Notations: we used the notation M(n;•) to denote the n-th row of matrix M. a first experiment the decomposition of one hundred noisy copies of a 2D Dirac function and the averaging of the final results, in an Ensemble EMD fashion [START_REF] Wu | Ensemble empirical mode decomposition: A noise-assisted data analysis method[END_REF]. Given ∈ R 64×64 , with (32, 32) = 1, we decomposed one hundred realizations of (i) = + βW (i) , with W (i) ∼ N (0, 1) and β = 0.001, into two modes and one residue.

The results are shown in Fig. 1. The latter figure confirms the orientation-independence of UOA-EMD.

B. AM-FM Images

In this subsection we present four different experiments decomposing AM-FM images, measuring both performance in component retrieval and computational time as a function of the only parameter λ. Moreover, we compare our results with those of the state of the art methods: Image EMD (IEMD) [START_REF] Linderhed | Image empirical mode decomposition: A new tool for image processing[END_REF], Prox-EMD [START_REF] Schmitt | 2D Prony-Huang transform: A new tool for 2D spectral analysis[END_REF] in its pseudo-2D implementation, and TV-G texture-geometry decomposition [START_REF] Gilles | Bregman implementation of Meyer's G-norm for cartoon + textures decomposition[END_REF].

IEMD uses sifting process to extract every mode, iterating until |e(x, y)| <ε , ∀(x, y),w h e r ee(x, y) stands for the mean envelope. The toolbox available on http://aquador.vovve.net/IEMD/ requires 0 <ε≤ 1. As ε becomes smaller, more sifting iterations are performed.

Prox-EMD achieves a decomposition by solving a convex non-smooth criterion which has texture and geometry terms. The toolbox is available on http://perso.ens-lyon.fr/nelly. pustelnik/logicielE.html.

TV-G texture-geometry decomposition provides a separation where the texture component is related with finer scales and thus can be considered as detail and the geometry component contents coarser scales and is related with the trend. The toolbox is available on http://www-rohan.sdsu.edu/ ~jegilles.

Both Prox-EMD and TV-G texture-geometry decomposition contain two parameters to be set: ν for the texture and ρ for the geometry. Based on [START_REF] Schmitt | 2D Prony-Huang transform: A new tool for 2D spectral analysis[END_REF], we use ν ={ 1, 5, 10, 50} and 1 ≤ ρ ≤ 20.

For the four methods (UOA-EMD, IEMD, Prox-EMD, and TV-G texture-geometry decomposition), we decompose the given image For all cases, we also show the decomposition results for the parameters that minimize the error, with the images being scaled to use the full colormap.

X = A + D into two components: X -→ [ Â, D].
1) Experiment 1: This experiment consists in decomposing X = X 1 + X 2 , with

X 1 = cos 6π(x + y) 2 + cos (18π(x -y)) X 2 = cos 1.5π(x + y) 2 , (13)
with

-1 ≤ x ≤ 1a n d-1 ≤ y ≤ 1.
The image has a size of 128 × 128 pixels. The results from the four methods are shown in Fig. 2. UOA-EMD and IEMD achieve very similar results, showing good recovering capabilities. On the other hand, Prox-EMD and TV-G show rather poor performances. Figure 3 presents both the errors when estimating the components, and the computational times. UOA-EMD gives the lowest error for two thirds of the considered ranges, and a constant computational time. IEMD evidences three different regimes: when the threshold value ε goes from 0 to 1, three different number of sifting iterations are used. Moreover, this first experiment shows that Prox-EMD and TV-G are unable to successfully recover AM-FM components for the image processed. Furthermore, the tuning of their two main parameters can be tricky. 2) Experiment 2: For this experiment we consider the decomposition of X = X 1 + X 2 , with

X 1 = cos 6π(x + y) 2 +cos (18π(x -y)) , if 0 ≤ y ≤ 1 cos 2π(x + y) 2 , if -1 ≤ y < 0 X 2 = cos 2π(x + y) 2 , if 0 ≤ y ≤ 1 cos 0.5π(x + y) 2 , if -1 ≤ y < 0 (14)
The bottom half of X 1 is equal to the top half of X 2 . Since X 2 is trend to X 1 , the same spatial frequencies are trend in one part of the image and detail in the other. A simple linear filtering is unable to separate the components. Results obtained with the four methods are shown in Fig. 4. UOA-EMD presents the best results, among others for the bottom half of the residue. The errors for Prox-EMD are close to those of UOA-EMD (Fig. 5). However, when analyzing the residue, Prox-EMD seems to favor piecewise constant solutions, which are not suited to catch AM-FM characteristics. The best performances for the computational times are obtained with UOA-EMD.

3) Experiment 3: In this experiment we introduce amplitude modulations: As before, the bottom half of X 1 equals the top half of X 2 . The decomposition given by the four methods are presented in Fig. 6. The errors and computational times are presented in Fig. 7. For this third experiment, results from UOA-EMD and IEMD are similar. The errors are of the same order of magnitude, except for ε = 1f o rI E M D where the error is significantly higher. Although errors for Prox-EMD are low for small values of geometry parameter ρ, the results deserve special attention. The Prox-EMD residue presents isolated squares of almost constant value, which are not present in the original component. The first mode of Prox-EMD presents in its bottom half a pattern of vertical, horizontal and diagonal lines, which do not appear in X 1 .

X 1 = ⎧ ⎪ ⎨ ⎪ ⎩ cos (24π(x + 0.1y)) , if 0 ≤ y ≤ 1 e -3(x 2 +y 2 ) (cos (6π(x + y)) + cos (6π(x -y))), if -1 ≤ y < 0 X 2 = ⎧ ⎪ ⎨ ⎪ ⎩ e -3(x 2 +y 2 ) (cos (6π(x + y)) + cos (6π(x -y))), if 0 ≤ y ≤ 1 e -3(x 2 +y 2 ) , if -1 ≤ y < 0 (15)
This effect is a clear consequence of a pseudo-2D approach, in which rows, columns and diagonal are separately constrained. The top half of TV-G first mode is almost equal to that of the original image X. Moreover, computational time remains almost constant for UOA-EMD, being thirteen times smaller than that of IEMD.

4) Experiment 4: As fourth and final experiment for this subsection, we masked the images from Experiment 3 with in order to test the behaviour of the methods on an image with a piecewise constant zone.

M = ⎧ ⎨ ⎩ 1, if x 2 + y 2 < 16 25 0, otherwise, (16)
The results are presented in Fig. 8. Same features as in the previous example appears. The errors (see Fig. 9) for IEMD increase, while the performance of UOA-EMD remains almost unchanged.

C. Real Data

Real standard images are widely used by image processors in order to test their performances. We herein present the decomposition of a 512 × 512 version of the image "Barbara" via the four methods. In every case, the image was decomposed into two modes plus residue. We used frequencies but conserves important features such as the folds of the woman's pants, of the cloth around her head and of the tablecloth. A portion of the original image is pointed out. The same zone in modes first and second is zoomed. In the zoomed zone of the first mode the stripes of the fabric can be appreciated, with no traces of the dark areas. The stripes of the fabrics are absent from this zone in the second mode, and only the folds can be seen. All the contrast of the original image is present in the residue.

On the other hand, the other three analyzed methods seem unable to capture such features. IEMD first mode includes not only the edges and fabric stripes but also coarser oscillations such as the folds of the fabrics. The second mode is too blurry for fine structures to be identified. It also contains most part of the image contrast, leaving the residue with almost no relevant information. Prox-EMD retrieves a first mode where some edges are hard to detect, such as those of the legs of the table. The second mode and the residue are difficult to analyze. While the second mode resembles a watercolour version of the image with no fine details, the residue is too blurry and no structure can be easily identified (we must remember that Prox-EMD does not provide an exact reconstruction). Finally, TV-G first mode is the most similar, except for the bright, to the original image among first modes. The legs of the table are significantly darker than the floor, proving more than just edge detection; the same can be said about the woman's chest. Second mode and residue provide little information about the original image (as with Prox-EMD, TV-G reconstruction is not exact).

D. Biomedical Application

We herein present the decomposition of a 256 × 256 medical image. The latter is extracted from the website of Harvard Medical School for Radiology at http://www.med. harvard.edu/AANLIB/cases/case1/mr1-tl4/029.html [START_REF] Wu | The multi-dimensional ensemble empirical mode decomposition method[END_REF]. This image, a T2-weighted magnetic resonance image (MRI), is shown in Fig. 11. It has been recorded in a woman suffering from a glioma. The decomposition of the data into three components obtained with UOA-EMD is shown in Fig. 11, where λ k = 1, 1 ≤ k ≤ 3 was used. We observe that the three components represent each a separate scale range: the first component corresponds to the highest frequencies (the finest texture); as the number of decomposition grows, the largest scales are discovered. The trend is observed on the residue. The first mode clearly reveals the contours of the lesion. These results show that UOA-EMD can also be useful for the medical field.

V. C ONCLUSION Local and auto-adaptive decompositions such as EMD have gained an increasing attention both for signals (1D data) and images (2D data). For the process of images, most EMD-based decompositions rely on an iterative procedure (sifting process) that has the drawback of lacking convergence guarantees. The recent optimization-based approaches to EMD, on the other hand, makes no use of such process, paving the road for mathematical foundations of EMD.

In this work we proposed a 2D extension of an unconstrained optimization approach to EMD (UOA-EMD), and compared its performance to three other decompositions methods: IEMD, Prox-EMD, and TV-G. Our results on synthetic and real data show that the proposed framework overcomes the existing algorithms: it is orientation-independent, shows a low temporal complexity, and gives better performances in decomposing AM-FM images.

For the artificial AM-FM images we confirmed the locality of the method. Applying the here proposed technique to a well-known standard image ("Barbara") the capabilities of our method, when different spatial frequencies are related to different structures, were illustrated. Finally, we showed a potential application of our method decomposing a biomedical image (MRI data).

The simplicity of the proposed framework leaves room for an nD extension (n > 2).

  and the coefficients c i are those necessary to compare the function D(z l ) with the plane intersecting D(z i ):

Fig. 1 .

 1 Fig. 1. Dirac function decomposition. Average results of one hundred decompositions of noisy 2D Dirac functions.

Fig. 2 .

 2 Fig. 2. Results for Experiment 1. Best results for the four analyzed methods when decomposing an image mixing components described in Eqs. (13).

For

  UOA-EMD and IEMD, we have X = Â + D;f o r Prox-EMD and TV-G we have X ≈ Â + D. Therefore, to measure recovering capabilities of UOA-EMD and IEMD we can use either the error for A or the error for D,s i n c e||A -Â|| F =| |D -D|| F ,w h e r e|| • || F stands for the Frobenius norm. For Prox-EMD and TV-G we use (||D -D|| F +||A -Â|| F )/2.

Fig. 3 .

 3 Fig. 3. Errors and computational times for Experiment 1 (see text for the details).

Fig. 4 .

 4 Fig. 4. Results for Experiment 2. Best results for the four analyzed methods when decomposing an image mixing components described Eqs. (14).

Fig. 5 .

 5 Fig. 5. Errors and computational times for Experiment 2 (see text for the details).

Fig. 6 .

 6 Fig. 6. Results for Experiment 3. Best results for the four analyzed methods when decomposing an image mixing those of Eqs. (15).

Fig. 7 .

 7 Fig. 7. Errors and computational times for Experiment 3.

Fig. 8 .

 8 Fig. 8. Results for Experiment 4. Best results for the four analyzed methods when decomposing an image mixing those of Eqs. (15) masked with Eq. (16).

λ 1 =

 1 λ 2 = 1 for UOA-EMD, ε = 0.5a n dK = 2f o rI E M D , and ν 1 = ν 2 = 1, ρ 1 = ρ 2 = 10 and K = 2 for both Prox-EMD and TV-G. The decomposition needed around 5 minutes for UOA-EMD, five and a half hours for IEMD, 17 minutes for Prox-EMD and 50 seconds for TV-G. The results of the decompositions for the four methods are shown in Fig. 10. The first mode of UOA-EMD captures the highest frequencies. The stripes of the woman's clothes and the tablecloth are perfectly conserved, along with the edges of the rest of the image. The second mode of UOA-EMD catches lower

Fig. 9 .

 9 Fig. 9. Errors and computational times for Experiment 4.

Fig. 10 .

 10 Fig. 10. Decompositions of Barbara. Results of the decompositions via the four methods.

Fig. 11 .

 11 Fig. 11. Decomposition of a magnetic resonance image with UOA-EMD, from fine to coarse scales.