Srikumar Ramalingam

Sofien Bouaziz

Peter Sturm

Matthew Brand

Geolocalization using Skylines from Omni-Images

We propose a novel method to accurately estimate the global position of a moving car using an omnidirectional camera and untextured 3D city models. The camera is oriented upwards to capture images of the immediate skyline, which is generally unique and serves as a fingerprint for a specific location in a city. Our goal is to estimate global position by matching skylines extracted from omni-directional images to skyline segments from coarse 3D city models. Our contributions include a sky segmentation algorithm for omni-directional images using graph cuts and a novel approach for matching omni-image skylines to 3D models.

Introduction and Previous Work

The skyline has long been a source of fascination for photographers, thought to identify the city as uniquely as a fingerprint. In this work, we take a step further and investigate its uniqueness for any given geospatial location. We propose a simple, robust and fast method to compute the geospatial location by matching the skyline imaged on a fisheye image with the skyline synthesized from a coarse 3D model of a large city. In other words, we solve the pose estimation problem for an omnidirectional camera in the actual world coordinate system. Pose estimation, although extremely challenging due to several degeneracy problems and inaccurate feature matching, is a well researched topic in computer vision. However most existing solutions are only proven on the smaller scale of the laboratory setup.

In the last few years, there has been an increasing interest in inferring geolocation from images [START_REF] Robertson | An image-based system for urban navigation[END_REF][START_REF] Zhang | Image based localization in urban environments[END_REF][START_REF] Yeh | Searching the web with mobile images for location recognition[END_REF][START_REF] Jacobs | Geolocating static cameras[END_REF][START_REF] Hays | Im2gps: estimating geographic images from single images[END_REF]. In [START_REF] Robertson | An image-based system for urban navigation[END_REF], Robertson and Cipolla showed that it is possible to obtain geospatial localization by matching a query image with an image database using vanishing vertical direction. Zhang and Kosecka showed accurate results in the ICCV 2005 computer vision contest ("Where am I?") using SIFT features [START_REF] Zhang | Image based localization in urban environments[END_REF]. Jacobs et al. used a novel approach to ge- olocate a webcam by correlating its images with satellite weather imagery at the same time. Hays and Efros used millions of GPS-tagged images from the web for georeferencing a new image [START_REF] Hays | Im2gps: estimating geographic images from single images[END_REF]. In contrast to most of these approaches that leverage on the availability of these georeferenced images, we use coarse 3D models from the web for geospatial localization: like georeferenced images, a large repository of coarse 3D models already exists for major cities in the world. Several 3D reconstruction algorithms have been proposed for the reconstruction of large urban scenes [START_REF] Teller | Calibrated, registered images of an extended urban area[END_REF][START_REF] Cornelis | Fast compact city modeling for navigation pre-visualization[END_REF][START_REF] Gallup | Real-time plane-sweeping stereo with multiple sweeping directions[END_REF]. Koch and Teller proposed a localization method using a known 3D model and a wide angle camera for indoor scenes by matching lines from the 3D model with the lines in images [START_REF] Koch | Wide-area egomotion estimation from known 3d structure[END_REF]. In contrast to their work, our work relies only on the skylines for geolocalization. Although, the idea of using skylines has been explored earlier in [START_REF] Meguro | Development of positioning technique using omni-directional ir camera and aerial survey data[END_REF][START_REF] Stein | Map-based localization using the panoramic horizon[END_REF], our approach is significantly different from them. An expensive infrared camera is used in [START_REF] Meguro | Development of positioning technique using omni-directional ir camera and aerial survey data[END_REF] instead of a visible one. There are several important differences between our approach and [START_REF] Stein | Map-based localization using the panoramic horizon[END_REF]. First, a human user input is required to extract the skyline in the case of [START_REF] Meguro | Development of positioning technique using omni-directional ir camera and aerial survey data[END_REF], whereas we develop an automatic algorithm using graph cuts. Second, a hash table is precomputed to match the skylines with the 3D model, whereas our algorithm synthesizes fisheye images on the fly 1 for matching. We show our basic idea of matching skylines from 3D models and omni-directional images in Figure 1.

Our main goal is to obtain highly accurate geolocation for an image. It has been shown formally that omnidirectional cameras can give much better accuracy in motion estimation than perspective cameras [START_REF] Daniilidis | The coupling of rotation and translation in motion estimation of planar surfaces[END_REF][START_REF] Brodsky | Directions of motion fields are hardly ever ambiquous[END_REF]. In the case of small rigid motions, two different motions can yield nearly identical motion fields for classical perspective cameras, which is not the case for omnidirectional images. Omnidirectional cameras have also proven advantageous in applications such as video conferencing, augmented reality, and surveillance.

In this work we used a fisheye lens with a field of view of about 183 • . However, we would like to propose a method that could work with all kinds of omni-directional cameras. For that reason, we chose to use a generic calibration approach which treats every camera as a mapping between a pixel and its corresponding projection rays. This general model has been used in various works [START_REF] Grossberg | The raxel imaging model and ray-based calibration[END_REF][START_REF] Neumann | Polydioptric camera design and 3d motion estimation[END_REF][START_REF] Pajdla | Stereo with oblique cameras[END_REF][START_REF] Pless | Using many cameras as one[END_REF][START_REF] Ramalingam | Towards complete generic camera calibration[END_REF][START_REF] Swaminathan | A perspective on distortions[END_REF], and is best described in [START_REF] Grossberg | The raxel imaging model and ray-based calibration[END_REF], where properties other than geometric ones are also considered. It was recently shown that the generic calibration algorithm outperforms standard parametric approaches in the case of very high distortions [START_REF] Ramalingam | Towards complete generic camera calibration[END_REF][START_REF] Dunne | A comparison of new generic camera calibration with the standard parametric approach[END_REF]. We use the calibration information to synthesize virtual fisheye views from the 3D model and match these with real fisheye images. To the best of our knowledge, we are not aware of any work in vision which synthesizes omni-directional images from 3D models for pose estimation. The synthesis is done using a pixel shader program implemented on a Graphics processing unit (GPU) for generating a large number of fisheye images in real time.

The various buildings blocks of this system are already novel contributions in their respective domains. First we propose a graph cuts based algorithm embedded with parameter learning for robustly detecting sky in omnidirectional images. The most closely related work, which already gives excellent results, is the geometric labeling problem to classify a given image into sky, buildings and ground [START_REF] Hoiem | Recovering surface layout from an image[END_REF][START_REF] Ramalingam | Exact inference in multi-label crfs with higher order cliques[END_REF]. However, some of their prior assumptions do not hold true in our model. First, the camera is not perspective. Second, the image is not taken by a person with the optical axis approximately parallel to the ground. In fact, the optical axis is perpendicular to the ground. Third, their algorithm expects the sky to be mostly on the top of the image. On the contrary, our images always have the sky at the center of the image. As a result of these differences, we chose to develop a tailor-made sky detection algorithm for upward looking omni-directional cameras. We then propose a novel method to match skylines from omni-directional images by synthesizing fisheye images using existing 3D models and calibration information, and its application to geolocalization.

Overview of our Algorithm

We show the various stages of our algorithm in figure 2. The first step is to detect the region corresponding to sky in the omndirectional images. In order to do this we use a graph cuts based algorithm, which is generally used to minimize a discrete quadratic pseudo-boolean energy function. We give the details of the features used to compute unary and pairwise costs of the energy function. The parameters of the function are usually fixed manually. As we expect our algorithm to work at various lighting conditions and places, we use a parameter learning algorithm to automatically compute the best set of parameters. The exact details of the algorithm are shown in section 2.1. The second stage is the calibration algorithm, which is used to compute the 3D projection rays corresponding to every image pixel in the omni-directional image. We also propose enhancements of the calibration algorithm, using plumb-line constraints to significantly improve its performance, cf. section 2.2. Third, we show a novel method to synthesize virtual fisheye views at various locations in a coarse 3D model. This uses pixel shaders executed on a GPU, and is very fast.

The next stage is the skyline matching algorithm where we compute the chamfer distance between the synthesized and real fisheye images. Finally in section 3, we show experimental results on fisheye images captured in the city of Boston. We also show a perturbation analysis on the shape of the skylines to study the robustness of our method and its superiority over perspective images.

Sky Detection

Given an omni-directional image, which is circular in the case of a fisheye model, we want to classify the pixels into sky and rest. This can be seen as a segmentation problem with two labels. The features that can be used for this segmentation can vary from simple RGB colorspace components to a wide variety of features like gradients, straight lines, vanishing points, etc. Our approach has two modules: a parameter learning method and a discrete optimization algorithm. In our problem we use graph cuts, which is both fast and highly successful in various vision problems like stereo, segmentation, and image restoration [START_REF] Szeliski | A comparative study of energy minimization methods for Markov random fields[END_REF]. An energy function, involving binary variables in unary and pairwise terms, is represented using a weighted graph whose minimum cut (computed using the maxflow algorithm) yields an energy-minimizing partition of the variables. The minimum cut separates the set of nodes in the graph into two sets, one belonging to a source node and one belonging to a sink node. Formally, let G = (V, E) be a directed graph with non-negative edge weights and two special nodes, namely, the source S and the sink T . The st-mincut algorithm partitions the set of vertices in V into two disjoint sets V S and V T , such that S ∈ V S and T ∈ V T . The special nodes S and T correspond to sky and rest respectively. We briefly introduce the energy function, whose parameters we are interested in learning. Let x i ∈ B = {0, 1} where i = 1, 2, ..n represent boolean random variables. In our problem, each of these variables could represent the boolean labels (sky and rest) of the pixels in the omnidirectional image. We use quadratic pseudo-boolean functions for representing our energy functions. These are nothing but energy functions of boolean variables that map a boolean vector to a real value and thus the name pseudoboolean. Let θ denote the parameters in our energy function. The parameter vector θ consists of the unary terms θ i;a and the pairwise terms θ ij;ab , where i, j = 1, 2, .., n and a, b ∈ B. These parameters are also referred to as unary and pairwise potentials. In contrast to many vision algorithms where these parameters are manually fixed, we compute them automatically. The unary parameter θ i;a can be seen as a pseudo-boolean function f : B → R that gives the cost when x i = a. Similarly, a pairwise parameter θ ij;ab is a quadratic pseudo-boolean function f : B 2 → R that gives the cost when x i = a and x j = b. The function mapping partitions to energies is then

E(x|θ) = i∈V {θ i;0 (1 -x i ) + θ i;1 (x i )} + (i,j)∈E {θ ij;00 (1 -x i )(1 -x j ) +θ ij;01 (1 -x i )x j + θ ij;10 x i (1 -x j ) +θ ij;11 x i x j } (1)
Our goal is to learn the parameters (Θ) automatically for our problem. We optimize the discriminating power of the model by estimating parameters that maximize the difference between ground truth labellings and all other labellings of a small number of manually labeled examples. Our method is similar to the maximum-margin network learning method using graph cuts [START_REF] Szummer | Learning crfs using graph cuts[END_REF]; we generate "near-miss" labeling, then estimate a parameter vector that maximizes the margin separating true labellings from the near misses in a convex optimization. This is done repeatedly; the process converges to an optimal parameter vector in a small number of iterations. The theoretical properties of our algorithm will be published in a companion paper.

In our current implementation problem we obtain the unary likelihood for sky and rest using their color values. We first estimate a Gaussian model for the classes sky and rest and compute the mean and covariance using manually segmented ground truth images. For a new test image we obtain the unary likelihood by computing the Mahalanobis distance to the sky and non-sky classes as shown in figures 3(b) and (c). We will reformulate our energy function by decomposing the unary parameters θ i;a as follows:

θ i;a = θ p a + l i θ l a ( 2
)
Once we have the likelihood cost l i for every node the unary parameters θ p a and θ l a are dependent only the label a. Similarly we assume that the pairwise parameters θ ij;ab are also independent of the associated nodes i and j and replace them by θ ab . Due to the problem nature, we assume that the pairwise matrix is symmetric. i.e. θ 01 = θ 10 . We denote the new parameter vector which we want to learn as

Θ. Θ = θ p 0 θ l 0 θ p 1 θ l 1 θ 00 θ 01 θ 11 (3)
The parameter vector Θ is then estimated via the standard convex program for linear SVMs, using vectors of unary and pairwise statistics from true labellings and near miss labellings as positive and negative examples. In order to guarantee that the estimated model supports optimal inference, we augment the convex program with a constraint on the pairwise terms in the parameter vector that guarantees submodularity. In the binary case the constraint is simply

θ 00 + θ 11 ≤ 2θ 01 (4)
The submodularity condition is the discrete analogues of convex functions in continuous domains.

Calibration

We briefly explain the main idea behind the generic calibration approach [START_REF] Ramalingam | Towards complete generic camera calibration[END_REF] (see Figure 4). Three images of a calibration grid are captured by the general camera that needs to be calibrated. The images are taken from unknown viewpoints. Every image pixel observes three 3D points in different calibration grids. These 3D points are obtained in three different coordinate systems. However, these points are collinear if they are expressed in the same coordinate system. That constraint enables us to compute the motion between the views and eventually the projection rays for the image pixels. Let Q, Q ′ and Q ′′ refer to the known points on the first, second and the third calibration grids respectively. Let the coordinate system of the first calibration grid be the reference frame. The pose of the second and the third calibration grids with respect to the first one be given by (R ′ , t ′ ) and (R", t") respectively. Let O represent the unknown optical center. We stack the three points O, Q and Q ′ (after transforming it to the reference frame) in the following matrix.

    O 1 Q 1 R ′ 11 Q ′ 1 + R ′ 12 Q ′ 2 + R ′ 13 Q ′ 3 + t ′ 1 Q ′ 4 O 2 Q 2 R ′ 21 Q ′ 1 + R ′ 22 Q ′ 2 + R ′ 23 Q ′ 3 + t ′ 2 Q ′ 4 O 3 Q 3 R ′ 31 Q ′ 1 + R ′ 32 Q ′ 2 + R ′ 33 Q ′ 3 + t ′ 3 Q ′ 4 O 4 Q 4 t ′ 4 Q ′ 4     (5)
The collinearity constraint will force the determinant of any submatrix, whose size is size 3 × 3, of the above matrix to vanish. In other words, we obtain four constraints by removing one row at a time. A similar collinearity constraint can be obtained from the 3D points on first and third calibration grids. Using all these constraints we can extract the motion variables (R ′ , t ′ ) and (R", t"). See [START_REF] Ramalingam | Generic imaging models: Calibration and 3d reconstruction algorithms[END_REF] for details on extracting the individual motion parameters. In practical calibration approaches, it is not possible to extract the projection rays for every pixel. Homography or bilinear interpolation can be used to compute the projection rays for other pixels.

Fisheye Synthesis

In figure 5 we show the various stages in synthesizing a fisheye image. OpenGL and DirectX pipelines do not support non-linear projections. In order to synthesize a fisheye image we first generate five binary perspective images, corresponding to the views in the cubemap shown in figure 5(c). As our algorithm uses 3D coarse models without any texture, we generate a cubemap simply by rendering a 3D model colored entirely in white; the resulting binary image is black in sky regions. We then use our calibrated ray-table to map the cubemap to a fisheye image, as shown in figure 5(d). The black region is the predicted shape of the sky. Some of the fisheye images synthesized at different places in a 3D model of Boston's financial center are shown in figure 5(e).

As the previous section indicated, our ray-calibrated view synthesis has the advantage that it does not "bake in" error that would arise from using a parametric lens model. In addition, a pixel shader program is implemented in GPU to generate these fisheye images at a very fast rate. This allows us to generate accurate fisheye images on the fly; there is no need to store images in a large database.

Skyline Matching

Using the graph cut sky detection algorithm we segment the sky in fisheye images and obtain the skyline as shown in figure 6(a). The predicted skylines corresponding to different locations in the 3D model are obtained using the method suggested in the previous section. Since the real-image segmentation and the predicted image are both high quality, a simple chamfer distance suffices to score the skyline match. During the chamfer matching, we vary the pose parameters of our virtual fisheye camera and obtain new skylines at various locations in the 3D model. This is done using a multi-resolution strategy: we start with synthesizing and matching fisheye images which are distributed at intervals of 5 meters and refine the location until we have achieved a positive match that is precise up to several centimeters. We show some matching results in figure 8.

Experiments

Our real experiments were all carried out in the financial district of Boston, which is shown in figure 9(b). In all our real experiments we used a Nikon Coolpix E8 Fisheye lens with a field of view of 183 • to capture about 300 images. We implemented all our programs using C++ and OpenGL and parallelized the GPU and CPU operations for enhanced performance. Also, we use a simple multi-resolution strategy as explained in section 2.4.

Table 1 gives the values of the parameters learnt with the approach described in section 2.1. Figure 7 shows a few results for the first step of the system, detecting sky regions in fisheye images. In the following, we comment on matching and localization performance.

In all of our test images, the match converged to solutions with very low chamfer scores of about 2 pixels (average distance of points on synthesized skylines, to skylines extracted in real images).

In order to investigate how uniquely the skylines "fingerprint" their locales, are we selected 16 random locations on the financial district. For each of these locations we generated 27 points by translations along X, Y and Z directions, separated each by 25cm. For the total of 432 images, we obtained a confusion matrix by computing the chamfer distance for all possible pairs (cf. figure 9(d)). We observed that for every point, the first 27 correct matches corresponded directly to its nearby locations. We conducted the same experiment in Boston's North End where all the buildings are residential and short as shown in figure 9(c). Our experiments demonstrated that the skyline matching algorithm can be robust even if the buildings are short. Indeed, the confusion matrix of the North End locations (cf. figure 9(e)) indicates even better location fingerprinting than that of the financial district.

θ p 0 1 θ l 0 -1 θ p 1 -0.03708 θ l 1 -0.
In order to assess the benefit of omni-directional cameras, we also conducted the same experiment of computing the confusion matrix with a perspective model in the financial district. We used a model with a field of view of 90 • and synthesized images at exactly the same locations in the financial district (tall buildings) that were chosen for the fisheye experiment. The confusion matrix (cf. figure 9(f)) indicates very clearly the large number of ambiguities between different skylines.

Conclusion

Our main goal in this work is to assess the suitability of skylines and omnidirectional image cameras for accurate geospatial localization. Our results clearly demonstrate that it is possible to outperform GPS measurements in urban scenes, which are known to be extremely problematic for commercial GPS units.

Nevertheless, we foresee potential problems due to inaccurate 3D models and loss of edges due to too much or too little sunlight or certain weather conditions. However, we expect to resolve ambiguities in such cases using other vision algorithms: interest point matching/tracking between consecutive frames, structure from motion algorithms, Kalman filtering, and other priors based on street map information. It is important to note that the current GPS systems have been improved and robustified for several years by addressing various complicated issues that even includes theory of relativity.

In a few years, we believe that it is possible to build an inexpensive, robust and accurate vision based GPS.

Figure 1 .

 1 Figure 1. An illustration of our motivating idea. (a) A car automounted with an omni-directional camera facing towards the sky. (b) The image captured by the camera with skylines marked. We match these skylines to those synthesized from 3D city models to compute the geospatial location.

Figure 2 .

 2 Figure 2. The various stages of our algorithm.

Figure 3 .

 3 Figure 3. (a) Original image. (b) Likelihood for the sky. (c) Likelihood for the rest of the image. The brighter values correspond to higher likelihood.

Figure 4 .

 4 Figure 4. (a) Generic calibration using three calibration grids. (b) and (c) show an original fisheye image and samples of distortion corrected regions, obtained using the calibration result.

Figure 5 .

 5 Figure 5. Our algorithm to construct a fisheye image. (a) and (b) show the 3D model and aerial image of an urban scene where we synthesize a fisheye image. (c) The cubemap generated from a given point in the 3D model. Note that in practice, we generate binary images; the slight shading here is just for visualization. (d) Fisheye image created from the cubemap. (e) Examples of skylines extracted from the fisheye images.

Figure 6 .

 6 Figure 6. (a) The skyline detected in a real fisheye image. (b) The chamfer distance map of the skyline of (a), with the best match of synthesized skylines superimposed on it.

Figure 7 .

 7 Figure 7. Sky detection results: Original and segmented fisheye images are shown.

Figure 8 .Figure 9 .

 89 Figure 8. Chamfer distance between the skylines from real and synthetic images, for changes in translation and rotation parameters in (a-c) and (g-i) respectively. The change in shape of the skylines is also studied for small perturbations in translation (X in (d),Y in (e), Z in (f)) and rotation parameters (Rx in (j), Ry in (k) and Rz in (l)). We use color-coding to display the change in shape of the skylines. The skylines drawn in blue, green and red show the initial, intermediate and final shapes. One unit in the model space corresponds to 25cms for translation and 1 degree for rotation. [Best viewed in color]

Table 1 .

 1 The learned parameters of the energy function shown in equation[START_REF] Brodsky | Directions of motion fields are hardly ever ambiquous[END_REF]. Note that the pairwise parameters satisfy the submodularity conditions, as imposed in our learning framework. The accuracy figure in the last row is the percentage of correctly classified sky pixels in the training set.

		0935
	θ 00	-0.5638
	θ 01	0.2487
	θ 11	0.3151
	Accuracy (%)	94.2

Acknowledgments:

The authors would like to thank Dr. Jay Thornton for useful discussions in 2D to 3D registration and car navigation. Srikumar Ramalingam would like to thank Prof. Philip Torr for the various discussions in graph cuts and geometric labeling problems.