G Ronsmans
email: gronsman@ulb.ac.be

B Langerock

C Wespes

M De Mazière

D Hurtmans

C Clerbaux

P-F Coheur

FORLI-HNO 3 8-years time evolution of stratospheric HNO 3 columns: investigation of the drivers of variability and of the link to O 3

Introduction

Nitric acid (HNO 3 ) is one of the main species involved in the stratospheric ozone cycles. Until recently, the distribution of HNO 3 in the stratosphere was available only from limb-sounding instruments, such as MLS, MIPAS or ACE-FTS. The IASI nadir looking thermal infrared instrument is now providing HNO 3 concentration distributions with unprecedented spatial and temporal sampling since 2007. Here, we briefly review the IASI observational capabilities for HNO 3 (in terms of vertical sensitivity and errors) and show the results of a validation exercise made through a comparison with ground-based FTIR measurements. We also provide an analysis of the time evolution of HNO 3 concentrations from IASI measurements since the end of 2007, through global distributions as well as latitudinal time series of HNO 3 columns. The first results of a multivariate regression analysis are shown in order to highlight the various factors responsible for HNO 3 spatial and temporal variability. The capabilities of the IASI instrument also allows for a joint analysis between HNO 3 and O 3 evolutions throughout the years, also briefly described hereafter.

-Fast Optimal Retrieval on Layers for IASI = near-real time processing chain, at ULB -Retrieval range: 860-900 cm -1 ; Retrieval scheme: Optimal Estimation Method (Rodgers, 2000) -Inversion on 41 layers, with a single a priori profile and covariance matrix -Data kept if: DOFS > 0.9; RMS of the spectral residual < 3×10 -8 W.m -2. sr -1 .m -1 ; Cloud coverage < 25% -See Hurtmans et al. ( 2012

  ) for detailed information Validation o AVK: same shape for all layers, covering the entire altitude range from the troposphere to the stratosphere • Around 1 level of information in both equatorial and polar regions • Maximum sensitivity in the stratosphere around 10-20 km altitude, depending on the latitude o DOFS= from 0.9 (poles) to 1.2 (equator) o ERRORS: vary from 0-10% at polar latitudes to 35-45% at the equator, mostly due to the interference with water vapor absorption lines Total error on the total HNO 3 column (%) • Comparison of FTIR and IASI time series (bottom right figure): • Retrieved columns (5-35 km) of both datasets are within the error range of one another • Differences between FTIR and IASI data (I-F)/F vary between -30 and 30% • Mean of differences (bias) = 4.0 %  IASI slightly overestimates concentrations compared with FTIR • Standard deviation = 9.7 %  Bias > std  differences not significant, compared to variability Choice of 6 stations, all part of the NDACC network Example for Thule station (Groenland)

  Regridding and smoothing of the FTIR vertical profiles (bottom left figure) give good agreement with IASI profiles

	•
	Stations	Bias (%)	Standard	#	R
		Smoothed	deviation	pairs
		FTIR	(%)
	Thule	4.0	9.7	151 0.84
	Kiruna	8.6	11.9	206 0.81
	Jungfraujoch	13.9	9.6	583 0.91
	Izaña	9.2	9.8	256 0.74
	Lauder	16.3	11.9	80	0.81
	Arrival Heights -0.8	15.9	75	0.78
	All stations	10.5	11.5	1351 0.93

Spatial & temporal distributions

Contact


Acknowledgements

1-MEASUREMENTS

Twenty-two years of Arctic ozone depletion observations and simulations.

Is there a trend ?


Peter Hedman

Tobias Ritschel

George Drettakis

Gabriel Brostow

Scalable Inside-Out Image-Based Rendering

Keywords:

Figure 1: Images from our method rendered in 1080p at 55 Hz on an Nvidia Titan X GPU. Input is an RGB-D video and 298 high-quality photos of 'Dr Johnson's house', London. With no wheelchair access to this floor, curators were keen to have their rooms digitized.

Introduction

The ability to capture and reproduce virtual versions of real indoor places is essential for many applications, such as virtual navigation of real-estate, museums, games, and safety training. While it is possible to obtain 3D reconstructions [START_REF] Choi | Robust reconstruction of indoor scenes[END_REF]] and display them with texture [START_REF] Waechter | Let there be color! Large-scale texturing of 3D reconstructions[END_REF], this lacks view-dependent effects. Even with perfect geometry this approach looks artificial, as highlights are baked in or missing completely. Gradually, imagebased rendering (IBR) is becoming an effective way to achieve both realism and interactivity. High-quality results already exist for small objects [START_REF] Lensch | Image-based reconstruction of spatial appearance and geometric detail[END_REF]] and for outdoor environments [START_REF] Fehn | Depth-image-based rendering (DIBR), compression, and transmission for a new approach on 3D-TV[END_REF][START_REF] Goesele | Ambient point clouds for view interpolation[END_REF][START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF]. A key to recent success in IBR is the use of per-view input image information, such as custom meshes and super-pixel over-segmentation, which preserve depth boundaries even with imperfect 3D reconstructions.

Indoor scenes present specific challenges. While capturing outdoor scenes or views around an object involves "outside-in" viewing (see [START_REF] Wei | Multi-view depth map estimation with cross-view consistency[END_REF] for globally consistent 3D in that setting), interiors typically require an "inside-out" approach. For the latter, capturing similar numbers of photos yields far less overlap in scene coverage. This creates much larger parallax and many occlusions at close distances, resulting in two challenges: the need for higher-quality 3D reconstruction, and for capturing significantly more input photos, slowing down rendering when using per-view information.

We propose a new indoor IBR algorithm addressing both these challenges by combining indoor-friendly depth sensors and multiview stereo (MVS) for improved reconstruction, and a scalable rendering algorithm which uses mesh simplification and tiling to accelerate free-viewpoint IBR of indoor scenes.

Depth sensors allow easy 3D acquisition of indoor scenes, but have several limitations, e. g., depth maps that do not align well with image edges, and low resolution compared with RGB cameras. Depthsensor fusion algorithms [START_REF] Newcombe | KinectFusion: Real-time dense surface mapping and tracking[END_REF]Nießner et al. 2013;[START_REF] Zhou | Elastic fragmens for dense scene reconstruction[END_REF][START_REF] Choi | Robust reconstruction of indoor scenes[END_REF] provide high-quality 3D reconstruction, but can fail in flat regions such as walls, ceilings and the floor or when viewing incandescent lights. In contrast, MVS can achieve better results for textured zones and lights, and allows high-quality alignment of meshes with image edges. We carefully combine these two modalities to accurately align a global depthsensor mesh to separately-captured images, specifically addressingand significantly increasing -indoor IBR quality and speed.

One way to categorize previous free-viewpoint IBR methods is whether they use a global mesh [START_REF] Buehler | Unstructured lumigraph rendering[END_REF][START_REF] Eisemann | Floating textures[END_REF][START_REF] Goesele | Ambient point clouds for view interpolation[END_REF], or, for more recent approaches, forward-projection of per-view information [START_REF] Zitnick | Stereo for image-based rendering using image over-segmentation[END_REF][START_REF] Kopf | Image-based rendering in the gradient domain[END_REF][START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF][START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF]]. Per-view information compensates for insufficient accuracy of the global mesh, which can lack entire regions or contains incorrect geometry. For indoor scenes, depth sensors can estimate a globally "consistent" mesh [START_REF] Newcombe | KinectFusion: Real-time dense surface mapping and tracking[END_REF][START_REF] Choi | Robust reconstruction of indoor scenes[END_REF]], but still have errors at object boundaries, producing visual artifacts during IBR. Per-view information is still needed, but the number of images required for high quality is very high.

There are two key challenges for IBR of indoor scenes: combine the global 3D mesh with per-view geometry and render high-quality novel views using this geometry, without degrading performance. Our approach creates high-quality per-view meshes optimized to align with depth and image edges, which we subsequently simplify. Our rendering algorithm introduces a specialized tiling data structure, which compactly references the per-view data. We introduce an adaptive blending cost to combine several input views, and derive upper bounds of the cost for each tile. We use this cost to prioritize view-dependent data used to render the novel view, further accelerating rendering. In summary we present two main contributions:

• Careful merging of multi-modal sensor data for global scene reconstruction, which makes it easy to capture indoor scenes for IBR, followed by an algorithm to produce per-view meshes that respect view-specific depth and color discontinuities.

• A fast rendering algorithm allowing high-quality freeviewpoint IBR for indoor scenes, which uses a blending approach based on conservative cost bounds for quality, and uses mesh simplification together with tiled rendering to decouple speed from the number of input images.

We show results on a variety of complex indoor scenes, demonstrating that our approach provides high-quality results and can scale to a higher number of input images, compared to previous work.

Previous Work

Here we discuss the most closely related previous work in reconstruction and rendering. We refer to the excellent survey by [START_REF] Shum | Image-based rendering[END_REF] for the broader context of IBR.

Geometry Reconstruction for IBR Several recent IBR methods [START_REF] Goesele | Ambient point clouds for view interpolation[END_REF][START_REF] Eisemann | Floating textures[END_REF][START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF][START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF]] use multi-view stereo (MVS) [START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF][START_REF] Goesele | Multi-view stereo for community photo collections[END_REF]] to reconstruct the 3D geometry of the scene, starting only with a sequence of photographs as input. In most cases, these methods use outdoor scenes as examples. Individual objects can be reconstructed at real-time rates with sufficient accuracy for IBR using visual hulls [START_REF] Matusik | Image-based visual hulls[END_REF]].

Indoor scenes have been digitized for IBR, e. g., [START_REF] Sinha | Piecewise planar stereo for image-based rendering[END_REF] focused on piecewise-planar geometry, while [START_REF] Furukawa | Reconstructing building interiors from images[END_REF] modeled shape using the Manhattan-world assumption. With manual intervention and a carefully crafted user interface, [START_REF] Sankar | Capturing indoor scenes with smartphones[END_REF] reconstructed plausible indoor scenes, captured even on a mobile device.

For indoor scenes, consumer-level depth sensors [START_REF] Newcombe | KinectFusion: Real-time dense surface mapping and tracking[END_REF][START_REF] Henry | RGB-D mapping: Using Kinect-style depth cameras for dense 3D modeling of indoor environments[END_REF][START_REF] Choi | Robust reconstruction of indoor scenes[END_REF]] can be used to recover 3D information. Elaborate hardware has been used to methodically capture large interior spaces, room by room, producing data that is hard to use in interactive systems [START_REF] Bahmutov | Efficient large scale acquisition of building interiors[END_REF]]. In the Sun 3D data set [START_REF] Xiao | SUN3D: A database of big spaces reconstructed using SfM and object labels[END_REF]], 3D reconstruction is done using Structure from Motion and RGB-D images, but does not use digital camera images, which can result in imperfect reconstructions especially at depth silhouettes. Using stereo cues to improve depth sensor images is a well studied topic, surveyed by [START_REF] Nair | A survey on time-of-flight stereo fusion[END_REF]. These methods seldom fuse multiple depth images, and refine depth using precalibrated binocular stereo. In contrast, our per-view refinement method works with fused depth images and unstructured multi-view capture.

Existing 3D reconstruction methods suffer from various inaccuracies, including missing or spurious additional geometry, misalignment between geometry and image boundaries, and excessively smooth or noisy geometry; these result in distracting visual artifacts for IBR [START_REF] Stich | Perception-motivated interpolation of image sequences[END_REF]. Combining textures from multiple views on a reconstructed object is challenging [START_REF] Lensch | A silhouette-based algorithm for texture registration and stitching[END_REF].

Reconstruction-based solutions to these problems include warping images to match the geometry [START_REF] Zhou | Color map optimization for 3D reconstruction with consumer depth cameras[END_REF], hiding seams between images while texturing [START_REF] Waechter | Let there be color! Large-scale texturing of 3D reconstructions[END_REF], or super-resolution approaches [Goldlücke et al. 2014]. In contrast, we carefully craft the 3D reconstruction to align RGB images with a global mesh, typically created by a depth sensor, then build per-view meshes to capture image and depth discontinuities with the accuracy needed for IBR. Large laser-scanned scenes can be rendered efficiently by converting them into local geometries, either into meshes to be rasterized [Arikan et al. 2014] or depth maps displayed with ray tracing [START_REF] Arikan | Multi-depthmap raytracing for efficient large-scene reconstruction[END_REF]]. These methods apply texture by hiding seams between high-resolution photographs projected onto the local geometry, but only in Lambertian scenes. While they focus on performance, our approach refines global geometry into per-view geometries to improve quality and is designed to reproduce view-dependent effects such as highlights.

Image-based Rendering

From the outset, Image-Based Rendering methods (e. g., [START_REF] Mcmillan | Plenoptic modeling: An image-based rendering system[END_REF][START_REF] Levoy | Light field rendering[END_REF][START_REF] Gortler | The lumigraph[END_REF]) capture the entire visual information of the scene using photographs. The Light Field [START_REF] Levoy | Light field rendering[END_REF]] and the Lumigraph [START_REF] Gortler | The lumigraph[END_REF] represent the dense information of light rays in a scene using different parameterizations. The Lumigraph introduced the idea of using geometry to assist in the generation of novel views. View-dependent texture mapping assumes the existence of a global mesh and all appearance variation is stored in textures [START_REF] Debevec | Efficient view-dependent image-based rendering with projective texturemapping[END_REF]]. Surface light fields are another parametrization of a light field for scenes containing opaque surface with view-dependent appearance [START_REF] Wood | Surface light fields for 3D photography[END_REF]]. There have been several methods that compress and render such data efficiently [START_REF] Chen | Light field mapping: efficient representation and hardware rendering of surface light fields[END_REF][START_REF] Vanhoey | Robust fitting on poorly sampled data for surface light field rendering and image relighting[END_REF], while Gigaray lightfields [START_REF] Birklbauer | Rendering gigaray light fields[END_REF]] use caching and tiling to accelerate computation. These methods operate on different data than ours: typically a single object or small baseline image sequences. High-quality IBR can be achieved through dense capture and by restricting the virtual camera to move on a plane [Aliaga et al. 2002] or to follow a track [START_REF] Uyttendaele | Image-based interactive exploration of real-world environments[END_REF]]. Our method supports free-viewpoint navigation and renders novel views that are substantially different from the input images.

The Unstructured Lumigraph Rendering (ULR) algorithm uses a set of complex weights [START_REF] Buehler | Unstructured lumigraph rendering[END_REF] to overcome limitations related to lack of geometry, input camera position, orientation and resolution. If the ULR approach is used per pixel, the cost of visiting all input pixels for each output pixel is prohibitive. Floating Textures overcome ghosting in ULR-like algorithms by using optical flow to improve blending, and also improve visibility [START_REF] Eisemann | Floating textures[END_REF]. ULR allows some degree of "free-viewpoint navigation" based on the global mesh which is used to back-project pixels in the novel view into the input images. Recently, [START_REF] Pujades | Bayesian view synthesis and image-based rendering principles[END_REF] present an interesting analysis of ULR by modeling the uncertainty in geometric reconstruction, but it is not adapted to our real-time rendering context. Alternative solutions include Ambient Point Clouds [START_REF] Goesele | Ambient point clouds for view interpolation[END_REF], which degrade image quality with a "non-photorealistic" look when reconstruction is unreliable, guided by epipolar constraints. However, in all the above cases imprecise or incomplete geometry and misalignment with input images results in significant visual artifacts in IBR, which worsen for indoor scenes.

A different class of approaches is based on the idea of using viewdependent information that preserves image discontinuities. These reduce artifacts at occlusion boundaries, by forward-projecting the over-segmented input images into the novel view [START_REF] Zitnick | Stereo for image-based rendering using image over-segmentation[END_REF][START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF][START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF]]. [START_REF] Bhat | Using photographs to enhance videos of a static scene[END_REF] use forward projection to enhance an existing video, allowing them to enforce temporal coherence and hide seams between input views using offline graph cuts. [START_REF] Kopf | Image-based rendering in the gradient domain[END_REF] forward project gradients and solve a Poisson problem to generate the image in the novel view, enabling smooth interpolation of reflections. Forward projection-based methods [START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF][START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF] produce high-quality free-viewpoint navigation mainly in outdoor scenes, typically when looking at a "backdrop" such as a building facade or a wall. For indoor scenes, forward-warping per-pixel geometry of the many input views required quickly becomes the computational bottleneck, given that a one-megapixel input image becomes a one-million vertex mesh.

A different -but direct -approach is to produce a novel view as an optimization [START_REF] Fitzgibbon | Image-based rendering using image-based priors[END_REF], such that the novel-view image matches the statistics of a set of examples. This approach has been demonstrated to work well for small changes of viewpoint and for dense image sets; it is however unclear how well the method adapts to the more challenging scenario we consider, with much larger changes in viewpoint and interactive rates. Initial results indicate that Deep Learning could be applicable to IBR [START_REF] Flynn | DeepStereo: Learning to predict new views from the world's imagery[END_REF], but in its current state this method suffers from blurring and high computational cost. Commercially backed systems like Google Jump and Facebook Surround 360 are emerging, which achieve highquality results by limiting the user to follow pre-defined paths. This keeps quality very high and makes on-set capture relatively easy; extending to our free-viewpoint context is much more challenging.

Finally, the quality of IBR critically depends on the quality and quantity of the images used. Many redundant or useless images can sometimes be worse than a smaller number of well-chosen views. Our tiled rendering approach is related to the view planning problem, applied to IBR [Vázquez et al. 2003], which can also be solved by putting a user into the loop [START_REF] Davis | Unstructured light fields[END_REF]]. We develop a compact and efficient solution, combining a pre-processing phase with a fast run-time step.

Overview

Our goal is to achieve real-time free-viewpoint image-based rendering of indoor scenes, using a collection of high-resolution digital color photographs and RGB-D video from a consumer-level depth sensor as input (Fig. 2a). Unlike traditional methods for RGB-D reconstruction [START_REF] Newcombe | KinectFusion: Real-time dense surface mapping and tracking[END_REF][START_REF] Dai | BundleFusion: Real-time globally consistent 3D reconstruction using online surface re-integration[END_REF]] and texturing [START_REF] Zhou | Color map optimization for 3D reconstruction with consumer depth cameras[END_REF][START_REF] Waechter | Let there be color! Large-scale texturing of 3D reconstructions[END_REF], we aim to reproduce view-dependent appearance such as highlights, and to render accurate images even in regions where a global 3D reconstruction of the scene has missing or inaccurate data.

Accurately aligning the photographs with the RGB-D images (Fig. 2b) and creating a global 3D reconstruction of the scene (Fig. 2c) is necessary for consistent rendering. However, these steps alone are not sufficient to achieve high quality with back-projection IBR (e.g., [START_REF] Buehler | Unstructured lumigraph rendering[END_REF][START_REF] Eisemann | Floating textures[END_REF]), since object boundaries in the photographs seldom align with the global geometry. Rather than relying on improved global geometry, we make a deliberate trade-off: To respect object boundaries we sacrifice global agreement and create per-view geometry for each high-resolution photograph (Fig. 2d). Now, we can render novel high-quality images by forward projecting per-view geometries into the novel view, and blending them using view-dependent blend weights.

Even with high-quality geometry, hundreds of input photographs are required to faithfully reproduce the view-dependent appearance of indoor scenes during free view-point navigation, given the many occlusions and large parallax at close distances. However, iterating over all photographs when rendering novel views is very expensive, especially for forward projection where each photograph corresponds to a per-view mesh with more than a million vertices. To overcome this problem, we present a culling strategy with three elements: Simplified per-view geometry (Fig. 2e), a tiling data structure to only render potentially visible geometry (Fig. 2f), and worst-case cost bounds to avoid rendering per-view geometry likely to be discarded during view-dependent blending (Fig. 2g).

Combining the insights above, our approach works as follows: During offline pre-processing we first use Structure-from-Motion to align the photographs with the RGB-D images in the same coordinate system and create global geometry using off-the-shelf surface reconstruction (Fig. 2b & 2c,Sec. 4.1). Next, we create per-view geometry that aligns well with input image edges. We use the global geometry to guide an iterative region-growing multi-view stereo method that creates per-view depth maps for each input photograph (Fig. 2d, Sec. 4.2). We then convert the depth maps into a simplified mesh representation suitable for rendering (Fig. 2e, Sec. 4.3). Finally, we store the resulting meshes in a tiling data structure to implement our culling strategy (Fig. 2f, Sec. 5.2). During run-time, we query the data structure to find tiles of relevant per-view geometry needed to form the novel view (Fig. 2g, Sec. 5.2), and derive bounds on the blending cost per tile to allow prioritization. Then, we forward project these tiles into the novel view and blend them using adaptive blend weights that automatically balance blurring and banding (Fig. 2h, Sec. 5.1). We now proceed to explain the 3D reconstruction (Sec. 4) and the rendering (Sec. 5) phases in detail.

3D Reconstruction

In the absence of perfectly accurate global geometry, we instead seek global geometry with sufficient quality to serve as initialization for local per-view reconstruction. In general, image-based rendering methods strive to be robust against imperfect 3D geometry. However, some problems are difficult to fix with rendering alone, such as big holes, "flying" geometry and disjoint reconstructions.

Indoor scenes are harder than outdoor in many ways, but they do allow for the use of Kinect-style depth sensors. We use such a depth camera, along with a high-resolution digital color camera because taken together, they will allow us to reconstruct both flat-looking parts of scenes (e. g., ceilings and walls) and highly textured areas. In practice, good global reconstructions are still elusive. Scene completeness is one major challenge, and accuracy is the other. We found that global reconstructions under-or (less often) overestimate the volume of scene elements, so they fail to align with edges in the input views. Broadly speaking, global reconstruction is designed to reward the average agreement between the estimated 3D shape and all views. But even small inaccuracies in the camera poses conspire to erode 3D shapes. Our goal of free-viewpoint IBR is less forgiving about missed details, but quite compatible with global inconsistencies. We thus separate reconstruction into two goals: finding a global 3D structure that brings all data into the same coordinate system, and constructing local per-view geometry, which aligns well with image edges in just the nearby input images.

To achieve this, we first use the RGB-D images to reconstruct a consistent global geometry of the scene (Sec. 4.1), which we then refine into view-specific depth maps for each high-resolution camera using multi-view stereo (Sec. 4.2). Finally, we convert the depth maps into simplified meshes amenable to efficient rendering (Sec. 4.3).

Global Geometry

In this step, we align the RGB and RGB-D images and consolidate the depth information into one consensual surface geometry. We found this procedure to work best for our input data, in terms of scene completeness. However, other geometry reconstruction methods (e. g., BundleFusion [START_REF] Dai | BundleFusion: Real-time globally consistent 3D reconstruction using online surface re-integration[END_REF], developed in parallel) might be worthwhile replacements, once available.

RGB-to-RGB-D Registration

We estimate relative camera poses between the high-resolution photographs and the color component of the RGB-D images using standard Structure-from-Motion (SfM) software [START_REF] Moulon | Adaptive structure from motion with a contrario model estimation[END_REF]], including radial distortion correction of the high-resolution photographs. A fundamental limitation with SfM is that it is unable to estimate the scale of the scene. Fortunately, it does provide us with the 3D locations of feature points and their 2D projections into the images (specifically the RGB-D images). Using this information, we align the SfM reconstruction with the (metric) scale of the RGB-D images: We use linear regression in a RANSAC loop to find the metric scaling factor that aligns the depths of feature points, projected into the RGB-D images, with the values stored in their depth maps.

Surface Reconstruction

We project all pixels in the RGB-D images out into a unified point cloud. Note that any point between an RGB-D camera and the surfaces it can see (determined by its depth pixels) is likely to be an outlier; we call this a visibility conflict. Inspired by [START_REF] Furukawa | Accurate, dense, and robust multiview stereopsis[END_REF], we use this insight to design a visibility filter that removes invalid points from the point cloud: If a point causes more than nmax visibility conflicts, we remove it from the cloud. As the depth captured by an RGB-D camera is less accurate further away, nmax is made a decreasing function of depth z in meters, i. e., nmax = 10m/z . Finally, we estimate normals for the remaining points in the cloud using RANSAC plane-estimation, and fuse them into a global mesh using Screened Poisson Surface Reconstruction [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF].

Local Per-view Geometry

For high-quality IBR, pixel-accurate alignment of depth with luminance edges is much more important than the precise depth value [START_REF] Stich | Perception-motivated interpolation of image sequences[END_REF][START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF]]. We achieve this by creating a different depth map for every input view, which might deviate from the global mesh, but respects image edges.

We use the global mesh created in the previous step to guide the creation of these high-quality per-input-view depth maps, which we then use for rendering novel views. For this purpose we employ a region-growing multi-view stereo method to align the global geometry with edges in the high-resolution images. We designed this method to exploit the global mesh in three ways: 1) to provide a good initialization of the surfaces in the scene, 2) to derive a reliable reconstruction cost function for optimization and 3) to determine visibility when selecting views for evaluating photo-consistency. Starting from pixels associated with planes (depths and normals) (a), we update the plane at a pixel (in red) by first extending the planes of neighboring pixels (b). We replace the plane at a pixel if any of the extended planes has a lower reconstruction cost (c). This is performed for all pixels in the image and repeated for 10 iterations.

The objective for our region-growing method is to find a plane (i. e., a depth and a normal) with a low reconstruction cost (defined below) for each pixel in an input image. To this end, we initialize the planes with the global geometry and iteratively refine them using the spatial propagation scheme from normal diffusion multi-view stereo (ND-MVS) [START_REF] Galliani | Massively parallel multiview stereopsis by surface normal diffusion[END_REF]. At every iteration, we create a set of plane candidates for each pixel by extending the planes of its neighbors, Fig. 3. Then, we replace the plane at a pixel with the candidate that has the lowest reconstruction cost. In all experiments, we perform 10 iterations of propagation. Please see [START_REF] Galliani | Massively parallel multiview stereopsis by surface normal diffusion[END_REF] for exact details of the propagation scheme.

Our reconstruction cost c(i) incorporates prior information from the global mesh with multi-view photo-consistency. More precisely,

c(i) = cp(i) + cm(i) + cs(i), (1)
is a sum of a photo-consistency, mesh, and smoothness terms.

We evaluate the photo-consistency cp(i) at a pixel by reprojecting a small (3 × 3) patch around the pixel via the candidate plane i into other images, where we compare both the colors and image gradients using the formulation in [START_REF] Galliani | Massively parallel multiview stereopsis by surface normal diffusion[END_REF]]. We use the global mesh to reliably select these other images. First, we rank the images based on the overlap they have with the input image when both are projected onto the global mesh. Then, we reproject the patch into the top six images not occluded by the global mesh.

The mesh prior cm(i) encourages candidate planes to stay close to the global mesh M. This is important for texture-less regions, where photo-consistency does not provide any extra information. Formally, where pi is the world-space location of the candidate plane and m is a point on the global mesh,

cm(i) = min m∈M |pi -m| 2 (50 mm) 2 . ( 2
)
The smoothness term cs(i) encourages the resulting depth map to be smooth for neighboring pixels i and j of similar color. Formally,

cs(i) = j∈N (i) exp -|ci -cj| 2σ 2 c min 0.1, |pi -pj| 2 2 (100 mm) 2 ,
(3) where N (i) is the set of neighboring pixels to i and ci and cj are the RGB colors at pixel i and j. In all experiments, we set σc = 0.1.

As a post-process, we remove pixels at unusual depths with the visibility filter described in Sec. 4.1, using nmax = 8. Finally, we use single-view information to align the depth maps to image edges where multi-view reconstruction fails. We filter each depth map with a cross-bilateral weighted median filter [START_REF] Ma | Constant time weighted median filtering for stereo matching and beyond[END_REF]] (guided by color, σc = 0.1) with a 4 px standard deviation.

Representation

Instead of storing full per-pixel per-input view depth maps, we create a compressed, per-view polygonal mesh for faster drawing and reduced storage. To this end, we convert every pixel in the depth map to a 3D vertex and form a grid mesh by connecting neighboring vertices with triangles. To avoid connecting foreground and background layers at occlusion boundaries, we do not create triangles at edges where the depth difference between the neighboring pixels is larger than 10 cm. We call those edges splits in the mesh.

Meshes made from depth maps have saw-tooth artifacts at discontinuities, which worsen for surfaces at glancing angles. To avoid this, we replace the depth of pixels at splits by a cross-bilateral average of their neighborhood depths, guided by both RGB and depth.

Finally, off-the-shelf mesh simplification [START_REF] Garland | Surface simplification using quadric error metrics[END_REF] is used to reduce the number of triangles in the mesh. We run the simplification until the mesh has an edge longer than 75 pixels in its input view, which typically yields a 50× to 100× reduction in the number of triangles. We flag vertices next to a split as boundary conditions that should not be altered during simplification. As can be seen in Fig. 9, this preserves crisp occlusion boundaries.

Rendering

We can now use the simplified per-view meshes for indoor IBR by blending many input views. To achieve high quality, including view-dependent effects such as highlights, we introduce an adaptive blending cost (Sec. 5.1). Directly rendering the per-view meshes is however prohibitively expensive, so we introduce two acceleration methods (Sec. 5.2). First, we present a spatial tiling data structure, limiting the per-view meshes used to render a given novel view. The number of tiles drawn is further reduced using a priority scheme, by deriving a bound on the blending cost for each tile. Details of our GPU rendering implementation (Sec. 5.3) conclude this section. Consequently, the minimum cost is high and the filter bandwidth becomes wide, blending all input views with nearly equal weights.

Blending

We form the novel view by locally blending the input images. Choosing appropriate blend weights is key to capturing view-dependent effects while minimizing ghosting artifacts. We do this by first resolving visibility of the novel view. Then, we compute a per-pixel cost for every input view. Finally, we reconstruct the novel view by adaptively blending the visible input views based on this cost.

Resolving Visibility As we form the final image by blending multiple per-view geometries, we cannot resolve visibility with a traditional depth test that only recovers the front-most surface. Instead we use a fuzzy depth test (detailed in Sec. 5.3), allowing us to consider all per-view geometries that represent the same surface, even if they do not perfectly align.

Cost

We derive an IBR cost that prioritizes which input views to blend for every pixel in the novel view. This is a function cIBR(yi, yn, x) of the camera position yi for the input view i, the camera position yn for the novel view and a location x on the perview geometry. Similarly to ULR [START_REF] Buehler | Unstructured lumigraph rendering[END_REF]] our cost is a weighted sum of two terms cIBR(yi, yn, x) = (1 -γ)ca + γc d , defined as follows:

ca = α(yi → x → yn) (Angle term) (4) c d = max 0, 1 - ||yn -x|| ||yi -x|| (Distance term) (5)
The angle term is the angle between the direction vectors from x towards the camera locations yi and yn. This helps to reproduce view-dependent effects, as it prioritizes input views observing the scene from viewpoints similar to the novel view. The distance term depends on the ratio of the distance between the current position x and the input view yi resp. the novel view yn. This reduces blur The tiling procedure, for a 3D scene with a single room (a), acquired from three input views (yellow, red, blue) to rendered in the novel view (in black). With a top-down perspective (b), we see that the blue and red input cameras can be seen in the black novel view, while the yellow one does not. Our method partitions the scene into a grid and splits the geometry into tiles (c). Now, tiles (e. g., all yellow ones and part of red and blue) not visible in the novel view get removed (d). The table in (e) shows the priority of each input view tile (columns) for all grid cells visible in the novel view (rows). Table entries are blank wherever the input view does not have a tile in the corresponding grid cell.

(undersampling) in the image by penalizing input views far from the surfaces visible in the novel view. Input views closer than the novel view are not penalized, and we correct for oversampling using mip maps and anisotropic filtering. In all experiments, we set γ = 0.1.

Bandwidth Selection

We form the novel view by blending together the input views using blend weights determined by IBR cost. Specifically, we compute the blend weight wi = exp(-ci/σ) for an input view i by applying an exponential filter kernel to the IBR cost ci = cIBR(yi, yn, x). This implies a trade-off illustrated in Fig. 4. We can choose a wide filter bandwidth σ resulting in blurring (i. e., giving all views a similar weight) or a narrow bandwidth implying banding (i. e., giving only one view a very large weight). Our key observation is that no single weighting function can capture two conditions that happen in practice: At one extreme, if many similar but incorrect views are available we need to select a wide bandwidth, i. e., to blur many images. Going for a single, but imperfect view would produce a sharp, but wrong result. At the other extreme, one view is much better than the others. Here, only this input view should be selected and others should not receive any weight as blurring many views would spoil the good one.

Similarly to variable kernel density estimation, we achieve this tradeoff by letting the filter bandwidth depend on the local costs of the input views. We adaptively set the filter bandwidth by scaling the filter bandwidth with the IBR cost cmin of the best novel view. In other words σ = σ k cmin, where σ k (set to 0.33 in all experiments) controls how quickly we transition from blending the input views to selecting the best one.

Tiled rendering

Even with the simplified mesh representation, drawing the per-view geometry of every input view into the novel view is slow. The key idea here is to avoid iterating over all input views when forming the novel view. We achieve this by operating on tiles, i. e., entire groups of input triangles and novel-view pixels. Similar ideas have been used in real-time graphics with tile-based shading [START_REF] Olsson | Tiled shading[END_REF].

The idea is shown in Fig. 5. Consider an example of four views (black, red, blue and orange) observing a room as seen in Fig. 5a. The black camera is the novel view, the others are input views. Fig. 5b is a top-view of the same scene. For every pixel in the novel view (black), we need to compute IBR weights for all input views (yellow, red and blue) so we can determine which images to fetch RGB values from. With forward warping, this means that all triangles from each input view need to be drawn into the novel view.

To reduce the computational cost, we partition the scene into a regular 3D grid as shown in Fig. 5. As a pre-process (Fig. 5c), we associate all triangles in the input views with the grid cells they intersect. We refer to the unique pairing of a grid cell with an input view as an input tile, i. e., the collection of triangles from a single input view that intersect a grid cell. At run-time, we save effort by only rendering input tiles visible to the novel-view (tiling and culling). In Fig. 5d, we see that this allows us to ignore the yellow view altogether as well as parts of the red and blue views.

However, because of our view-dependent blending, only a few visible input tiles will actually contribute to the final image. We therefore strive to render a small subset of the input tiles that ensures a low IBR cost for all pixels in the novel view. To achieve this, we sort the input tiles according to their worst-case IBR cost and draw only the best ones in each visible grid cell. The exact number varies between grid cells, as we stop drawing input tiles when we determine that a cell has been sufficiently covered. Predicting the worst-case IBR cost is part of our contribution, which focuses the rendering effort to the input tiles that actually affect the final image.

We next discuss all steps in detail.

Tiling and Culling

We store the input tiles in a 3D grid with a resolution of 32 × 32 × 32 that covers the scene. During preprocessing, we create the input tiles by conservatively voxelizing the per-view geometry into the grid. As a consequence, each triangle in the input geometry may belong to multiple input tiles.

At run-time, we need to find the grid cells visible to the novel view. This is done in two passes. First, we draw the coarse global geometry into the framebuffer of the novel view. Second, we process all framebuffer pixels in parallel and mark all grid cells they intersect as visible using an atomic operation. To account for the mismatch between the coarse global geometry and the per-view input geometry, we perform this intersection conservatively by inflating each pixel in the framebuffer into a cube with a side length of 5 cm. Priorization To select which input tiles to render, we sort them within each visible grid cell according to their worst-case IBR cost, i.e., an upper bound for cIBR (yi, yn, x) between the input view yi and the novel view yn where x is allowed to be anywhere in the grid cell. Concretely, x is inside an axis-aligned bounding box (xmin, xmax), see Fig. To find an upper bound for ca = α(yi → x → yn), we exploit that the angles in a triangle sum to π radians. Consider the triangle between x, yi and yn. We find lower bounds αmin(yi) and αmin(yn)

x min x xc x max y i y n α(y i → x → y n ) r y i y n α(x c → yi → y ) n sin -1 (r/||y i → xc||) b) a)
Figure 6: We derive the upper bound on the angle α(yi → x → yn), where yi is an input view, yn the novel view and x an point in a grid cell (a, Eq. 6). We find the bound implicitly through lower bounds for the other angles in the triangle. Here, we use the bounding sphere around the grid cell (center xc, radius r) to bound α(x → yi → yn) (b, Eq. 7).

for the other two angles in the triangle and form the upper bound:

c max a = π -αmin(yi) -αmin(yn). (6)
Fig. 6b shows how we compute the lower bound αmin(yi) for one of the other angles α(x → yi → yn). First, we convert the bounding box (xmin, xmax) to its encompassing bounding sphere centered at xc with the radius r. Now,

αmin(yi) = α(xc → yi → yn) -sin -1 r ||yi -xc|| . (7)
We derive an upper bound for the distance term c d based on the maximum distance dmax(y) and minimum distance dmin(y) between a point y and a bounding box (xmin, xmax). Using dmin(yn) as lower bound for the numerator and dmax(yi) as an upper bound for the denominator, we see that

c max d = max 0, 1 - dmin(yn) dmax(yi) . (8)
Drawing Once all input tiles have been sorted, we determine how many are to be rendered from each grid cell. Ideally, we would conserve bandwidth by rendering as few input tiles as possible.

However, this may result in holes in the novel view; the triangles in an input tile are not guaranteed to cover a grid cell, e. g., due to splits (Sec. 4.3) or image boundaries. We refer to an input tile without any of these discontinuities as complete.

For each grid cell, we render the sorted input tiles until three complete tiles or a maximum of 12 input tiles have been rendered. We found 12 to be a good trade-off between completeness and compactness. In our experiments this typically reduces the number of drawn primitives to fewer than 10 % of the original per-view meshes.

Implementation

Four geometry passes are required to render a novel view. Deferred shading cannot be used to avoid multiple passes, as the color of every output pixel is determined by blending several input primitives.

As described in Sec. 5.2, the first pass renders the global geometry to find the visible grid cells. These cells are downloaded to the CPU, which prioritizes the input tiles to render in the remaining three passes. We store the geometry associated with all input tiles in a single large vertex buffer object, allowing us to perform the other passes using a single glDrawIndirect OpenGL instruction.

The final three passes implement the blending algorithm described in Sec. 5.1. First, we recover the front-most surface with a depth pre-pass, which enables us to perform a fuzzy depth test: In the remaining two passes any surface fragment far enough (10 cm) from the front-most surface is discarded. The second pass finds the minimum IBR cost cmin of all input views that project onto each pixel. This is required for the per-pixel filter bandwidth σ explained in Sec. 5.1. Finally, the last pass forms the image, blending together the input views using the bandwidths computed in the previous pass.

Experiments

We compare our approach to existing baseline algorithms, using published code wherever possible. Most existing rendering systems were meant for outside-in scenes, so results are presented on our own data. Here and in the supplemental video, we provide qualitative comparisons of our reconstruction and rendering phases, and quantitative results for speed tests.

Data Sets

We captured our own scenes as we needed both highquality photos and depth-maps with good coverage of all interesting surfaces in a variety of indoor environments. We captured 150-300 RAW photos using a digital camera (Sony NEX-C3 at 1228 × 816 or Canon EOS 550D at 1296 × 864). We also recorded an RGB-D video using the Asus Xtion PRO depth-camera, and sub-sampled in time to yield 150-450 640×480 RGB-D images. In a pre-processing step, images were color-harmonized using Adobe Lightroom.

In Fig. 1, 7, 8 and the supplemental video we show results in the following scenes: CREEPY ATTIC: A small (5 × 4 × 4 m 3 ) attic in an old building containing textiles and artwork. There is a prominent object (doll on a chair) in the middle of the room, which clearly shows the need for per-input view local geometry since under-or overestimated global geometry results in clear IBR artifacts. DORM ROOM: A small bedroom scene (4 × 4 × 5 m 3 ) in a student dorm with prominent textureless walls. MUSEUM: A preserved house from the 17th century. This is a large scene (16 × 6 × 5 m 3 ) exhibiting textureless walls, glossy tabletops and paintings behind mirror-reflective glass. PLAYROOM: A medium-sized (6 × 6 × 5 m 3 ) livingroom cluttered with toys for young children. Aside from large, textureless surfaces this scene contains many small geometric details that cannot easily be captured using geometry reconstruction alone. BOOK SHOP: A large (11 × 9 × 5 m 3 ) basement in a book shop with difficult occlusion characteristics (aisles separated by bookshelves). The books present a challenge for reconstruction and rendering as they are often non-diffuse and textured with high-frequency details. READING CORNER: A small (7 × 4 × 4 m 3 ) reading corner in an academic book store. The scene contains a leather chair with strong view-dependent effects.

Rendering Comparison

In the accompanying video and Fig. 8 We compare our algorithm to [START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF], and also to ULR [START_REF] Buehler | Unstructured lumigraph rendering[END_REF] with improved visibility akin to floating textures [START_REF] Eisemann | Floating textures[END_REF]]. Compared to [START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF], our mesh-based approach tends to preserve complex shape boundaries better since super-pixels are warped independently in their method. This is illustrated in the left side of Fig. 8, where the arm and leg of the doll show severe ghosting using [START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF]]. Compared to ULR, our approach preserves the shape of detailed objects (the doll, middle of Fig. 8) or hard-to-capture objects like the lamp (Fig. 8 on the right).

Geometry Comparison

In Fig. 9 we compare the geometry reconstructions from different components of our system. We see that the global reconstruction alone is complete but tends to overor underestimate the size of foreground objects. We include ND-MVS [START_REF] Galliani | Massively parallel multiview stereopsis by surface normal diffusion[END_REF]] (Local only) as an upper bound for the quality can be achieved using local geometry reconstruction alone, i. e., without the mesh prior and the improved visibility tests in Sec. 4.2. This method uses stochastic search to find the planes in the scene as it does not rely on global geometry for initialization. In general it preserves details better than the global reconstruction, but is often incomplete (shown as white) as textureless regions are difficult to reconstruct using multi-view stereo. The geometry produced by our local per-view refinement (Sec. 4.2) is complete and aligns well with image edges. Note how it corrects occlusion boundaries grossly misestimated in the global reconstructions (circled in black):

The pillar (Fig. 9a); the cloth (Fig. 9b); the chair legs (Fig. 9c and9d); consistent vertical edges of a bookshelf (Fig. 9e) and the back of the armchair (Fig. 9f). Correctly reconstructing these boundaries is important for IBR, as getting them wrong causes foreground color to be displayed on background geometry or vice-versa.

Performance Tbl. 1 shows show statistics of our scenes and rendering. We measured the average frame time at 1080p in each scene as the virtual camera moves along a predetermined path. Without tiling, rendering takes 31-60 ms on our high-end machine (Desktop PC; Nvidia GTX Titan X), and 209-388 ms on the low-end machine (Laptop; Nvidia GTX 660M). Tiling provides a speedup of 1.6× to 3× depending on the scene and hardware. Interestingly, tiling reduces the number of rendered triangles more drastically than frame time. This is explained by the overhead from the tiling process itself, particularly from prioritization; the only CPU component of our rendering algorithm. Fig. 11 shows how the performance of our tiled algorithm scales with the number of input images, demonstrating a sub-linear trend where the frame time plateaus after 150 images.

Conclusion and Future Work

We have presented a new IBR method that first builds a high-quality global mesh from a depth sensor, and then registers high-resolution photos to this mesh. We then build per-input-view meshes guided by the global reconstruction. Our approach combines a global/local representation of the scene to introduce a novel IBR algorithm which introduces a new blending approach, a compact representation of the per-input-view meshes and uses a tiled rendering algorithm which scales well with the number of input images.

Limitations Our approach is ultimately limited by the quality of the initial capture; evidently SfM needs to succeed for all the input images for our method to work. We suffer from the same limitations as all 3D reconstruction methods, and in particular for scenes with glass, which are problematic for modern depth cameras [START_REF] Choi | A large dataset of object scans[END_REF]. Even in easier cases, feedback during capture would be helpful, since unobserved regions are usually completed as large blobby surfaces by the Poisson reconstruction. Capture problems could possibly be addressed with specialized on-set feedback about the capture progress, in the spirit of [START_REF] Davis | Unstructured light fields[END_REF]]. We show typical artifacts due to incorrect reconstruction in Fig. 10.

Handling much larger scenes with thousands of input photographs is currently a challenge. A clever disk and main memory caching scheme is required. Finally, our method is designed for relatively wide-baseline data sets. For light-field like densities, the CPU overhead for treating the tiling structure could become a bottleneck.

Future Work

Our tiling approach could be adapted to superpixelbased methods [START_REF] Chaurasia | Depth synthesis and local warps for plausible image-based navigation[END_REF][START_REF] Ortiz-Cayon | A Bayesian approach for selective image-based rendering using superpixels[END_REF], improving their performance and quality. Our approach is wellsuited to indoor scenes, which present specific difficulties, but are well-adapted to depth-sensor capture. However, it is very likely that our rendering method will benefit outdoor scenes, especially in the context of wide-angle viewing, e. g., for head-mounted displays, which requires a large number of input images.

  Overview of our algorithm: (a) The user captures high-quality inside-pointing-outward RGB photos, shown here as three different camera poses. Separately, they also record a low-resolution RGB-D video by walking around the scene. (b) The high-resolution photos are used to reconstruct a sparse 3D point cloud, to which the depth-maps are registered. (c) All the depth maps are fused into a global geometry, which smooths away many important details, but captures a consensual surface approximation. (d) Re-using the high-res photos, we compute a high-quality local mesh for each one. (e) The many local meshes are simplified to reduce triangle count. (f) The global geometry is partitioned into tiles (in red) that organize the scene's visibility with respect to input views. (g) At run-time, the per-view meshes are culled, leaving only those predicted as relevant for rendering a novel view (pictured as a black camera). (h) The relevant per-view geometry is rendered.

Figure 3 :

 3 Figure 3: The spatial propagation scheme of normal diffusion MVS.Starting from pixels associated with planes (depths and normals) (a), we update the plane at a pixel (in red) by first extending the planes of neighboring pixels (b). We replace the plane at a pixel if any of the extended planes has a lower reconstruction cost (c). This is performed for all pixels in the image and repeated for 10 iterations.

Figure 4 :

 4 Figure4: Adaptive bandwidth selection for two cases: One with a single good input view (top) and the other with many mediocre ones (bottom). In the "select" case, the first view matches the novel view well and has a low IBR cost. With a low minimum cost the filter bandwidth becomes narrow, assigning most of the blend weight to the best view. In "blend" case, all views have a high cost. Consequently, the minimum cost is high and the filter bandwidth becomes wide, blending all input views with nearly equal weights.

  Figure5: The tiling procedure, for a 3D scene with a single room (a), acquired from three input views (yellow, red, blue) to rendered in the novel view (in black). With a top-down perspective (b), we see that the blue and red input cameras can be seen in the black novel view, while the yellow one does not. Our method partitions the scene into a grid and splits the geometry into tiles (c). Now, tiles (e. g., all yellow ones and part of red and blue) not visible in the novel view get removed (d). The table in (e) shows the priority of each input view tile (columns) for all grid cells visible in the novel view (rows). Table entries are blank wherever the input view does not have a tile in the corresponding grid cell.

  6a. Using separately computed upper bounds for the angle term ca ≤ c max a and the distance term c d ≤ c max d , we form the upper bound cIBR(yi, yn, x) ≤ (1 -γ)c max a + γc max d .

Figure 7 :Figure 8 :

 78 Figure 7: Example results from the supplemental video (top to bottom): BOOKSHOP, PLAYROOM, DORM ROOM and READING CORNER.

Figure 9 :

 9 Figure 9: Different geometry reconstructions. In each scene we compare in the highlighted inset: High-quality RGB; Global reconstruction (Sec. 4.1); Our local-global reconstruction (Sec. 4.2); ND-MVS [Galliani et al. 2015] (Local-only), an upper bound for local-only reconstruction quality. Our local-global approach combines the best of both; complete reconstructions that align well with image edges (see the black circles).

Acknowledgements

Thanks to Dr Johnson's Museum, Judd Books, Waterstones on Gower Street, and Kathryn Emmett-Brostow for being so accommodating with our filming. Thanks to Rodrigo Ortiz-Cayon for computing the comparisons. The authors are grateful for the support of EU project CR-PLAY (no 611089) www.cr-play.eu, EPSRC EP/K023578/1, and project SEMAPOLIS (ANR-13-CORD-0003).

Table 1: Evaluation of our tiled IBR in six scenes with different levels of complexity. In each scene, we record the number of primitives needed when using our tiled method, and contrast it to a brute-force method rendering all local meshes. We render at 1080p and measure frame time in milliseconds, both on a high-end Machine A (Desktop with a Nvidia GTX Titan X) and a low-end Machine B (Laptop with a Nvidia GTX 660M). We fix the camera position in each scene and measure performance of our tiled rendering as the number of input images increases.

Scene


Gabin Personeni

Marie-Dominique Devignes

Michel Dumontier

Malika Smaïl-Tabbone

Adrien Coulet

Discovering ADE associations from EHRs using pattern structures and ontologies

Patient Electronic Health Records (EHRs) constitute an essential resource for studying Adverse Drug Events (ADEs). We explore an original approach to identify frequently associated ADEs in subgroups of patients. Because ADEs have complex manifestations, we use formal concept analysis and its pattern structures, a mathematical framework that allows generalization, while taking into account domain knowledge formalized in medical ontologies. Results obtained with three different settings show that this approach is flexible and allows extraction of association rules at various levels of generalization.

INTRODUCTION

Adverse Drug Events (ADEs) occurs unevenly in different groups of patients. Their causes are multiple: genetic, metabolic, interactions with other substances, etc. Electronic Health Records (EHRs) have been successfully used to detect ADEs [START_REF] Lependu | Pharmacovigilance Using Clinical Notes[END_REF]. We further hypothesize that mining EHRs may reveal that subgroups of patients sensitive to some drugs are also sensitive to others. We propose a method to identify frequently associated ADEs in these subgroups. Because the way ADEs manifest and are reported is complex and highly variable, we use an extension of formal concept analysis called pattern structures [START_REF] Ganter | Pattern Structures and Their Projections[END_REF] in combination with ontologies to enable generalization. We experimented on a dataset of EHRs from patients diagnosed with Systemic Lupus Erythematosus (SLE), a severe autoimmune disease. Such patients frequently experience ADEs as they often take multiple and diverse drugs indicated for SLE or derived pathologies [START_REF] Vasudevan | Established and Novel Treatments for Lupus[END_REF]. The SLE EHRs were extracted from STRIDE, an EHR data warehouse of Stanford Hospital and Clinics [START_REF] Lowe | STRIDE-An Integrated Standards-Based Translational Research Informatics Platform[END_REF].

MATERIAL AND METHODS

Data Corpus

Our data corpus is a set of 6,869 anonymized EHRs of patients diagnosed with SLE. It documents about 451,000 hospital visits with their relative dates, diagnoses encoded as a ICD9-CM phenotype codes (International Classification of Diseases, Clinical Modification) and drug prescriptions as a list of their ingredients, represented by their RxNorm identifiers. We use the term "drug" to denote an active ingredient rather than a commercial drug. As we aim to mine frequently co-occurring ADEs, we first need to identify such events in EHRs, and then select patients with at least two ADEs.

We first establish a list of ADE candidates, considering each patient EHR separately. From each two consecutive visits in the EHR, we extract the set of drug 𝐷 𝑖 prescribed during the first visit and the diagnoses 𝑃 𝑖 reported during the second. The interval between the two consecutive visits must be less than 14 days, as it is reasonable to think that a side effect should be observed in such a time period after prescription. Moreover, Table 1 shows that increasing that interval does not significantly increase the number of patients in our corpus. An ADE candidate 𝐶 𝑖 is thus a couple of sets 𝐶 𝑖 = (𝐷 𝑖 , 𝑃 𝑖 ). We retained in 𝑃 𝑖 only phenotypes reported as a side effect for at least one drug of 𝐷 𝑖 in the SIDER 4.1 database of drug indications and side effects [START_REF] Kuhn | The SIDER database of drugs and side effects[END_REF]. We then remove candidates where 𝑃 𝑖 is empty. Furthermore we remove an ADE candidate (𝐷 1 , 𝑃 1 ) if there exists for the same patient another ADE candidate (𝐷 2 , 𝑃 2 ) such that 𝐷 1 ⊆ 𝐷 2 : indeed, reiterated prescriptions of drugs may indicate that they are safe for this patient.

2

After filtering, we obtain a total of 3,286 ADEs from 548 patients presenting at least 2 ADEs. In Table 2, we present examples of ADEs that could be extracted from the EHRs, that will serve as a running example.

Formal Concept Analysis and Pattern Structures

Formal Concept Analysis (FCA) [START_REF] Ganter | Formal Concept Analysis: Mathematical Foundations[END_REF]) is a mathematical framework for organizing a dataset as a concept lattice, i.e., a hierarchical structure in which a concept represents a set of objects sharing a set of properties. It offers facilities for data mining of association rules.

In classical FCA, a dataset is a set of objects, where each object is described by a set of binary attributes. Pattern structures generalize FCA so that it can be performed on a set of objects with descriptions of any nature, such as sets, graphs, intervals, annotations from ontologies [START_REF] Ganter | Pattern Structures and Their Projections[END_REF][START_REF] Coulet | Using pattern structures for analyzing ontology-based annotations of biomedical data[END_REF]. A pattern structure is a triple (𝐺, (𝒟,⊓), 𝛿), where:  𝐺 is a set of objects, 𝒟 is a set of descriptions, 𝛿 is a function that maps objects to their descriptions.

 ⊓ is the meet operator defining a partial order ≤ ⊓ on elements of 𝒟 , such that, 𝑋 ⊓ 𝑌 is the most specific description that is more general than both 𝑋 and 𝑌 . Consequently, 𝑋 ≤ ⊓ 𝑌 denotes that 𝑌 is more specific than 𝑋 , is by definition equivalent to 𝑋 ⊓ 𝑌 = 𝑋 . Generalization on object descriptions is performed through the use of the meet operator. This is illustrated for our particular operators in Section 3. In pattern structures, the derivation operator . □ defines a Galois connection between sets of objects and descriptions, as follows:

𝐴 □ = ⊓ 𝑔∈𝐴 𝛿(𝑔) for a set of objects 𝐴 𝑑 □ = {𝑔 ∈ 𝐺 | 𝑑 ≤ ⊓ 𝛿(𝑔)}

for a description 𝑑

A pattern concept a pair (𝐴, 𝑑) verifying 𝐴 □ = 𝑑 and 𝑑 □ = 𝐴.

In our case study, 𝐺 is the set of patients that are related through 𝛿 to the description of their adverse drug events from 𝒟 . Section 3 describes different experiments using pattern structures, each providing their own definition of the triple (𝐺, (𝒟,⊓), 𝛿).

Medical Ontologies

We use two medical ontologies, considering only their class hierarchy, to uncover more general phenotype and drug descriptions: ICD9-CM that describes classes of phenotypes and the Anatomical Therapeutic Chemical Classification System (ATC) that describes classes of drugs. We use only the three most specific levels of ATC: pharmacological subgroups, chemical subgroups and chemical substances.

EXPERIMENTS

In this Section, we describe three experiments using pattern structures to extract association rules between ADEs. Each experiment defines a different representation of patient ADEs, making increasing use of ontologies.

First Pattern Structure Experiment

We define here the pattern structure (𝐺, (𝒟 1 ,⊓ 1 ), 𝛿 1 ) : objects of the pattern structures are patients, and a patient description of 𝒟 1 is a vector of sub-descriptions, with firstlevel ICD classes as dimensions. Each sub-description is a set of drug prescriptions, i.e., a set of sets of drugs. For instance, considering only the two ICD classes in Table 3: Sub-descriptions are associated to a first-level ICD class to represent ADEs: the patient presents a phenotype of that class after taking a prescription in that sub-description. We define sub-descriptions as sets of prescriptions, where none of the prescriptions are comparable to each other by the partial order ⊆. We then define the meet operator ⊓ 1 , such that, for every pair of descriptions (𝑋, 𝑌) of 𝒟 1 :

𝛿 1, ICD 580-629 (P1) = {{prednisone}, {acetaminophen}} 𝛿 1, ICD 710-739 (P1) = ∅
𝑋 ⊓ 1 𝑌 = max(⊆, {𝑥 ∩ 𝑦 | (𝑥, 𝑦) ∈ 𝑋 × 𝑌})
where max(≤ 𝑖 , 𝑆) is the unique subset of maximal elements of a set 𝑆 given any partial order ≤ 𝑖 . Formally, max(≤ 𝑖 , 𝑆) = {𝑠 | ∄𝑥. (𝑠 ≤ 𝑖 𝑥)} . In the present case, it retains only the most specific drug prescriptions in the description. For instance, given four drugs 𝑑 1 through 𝑑 4 :

{{𝑑 1 , 𝑑 2 , 𝑑 3 }} ⊓ 1 {{𝑑 1 , 𝑑 2 }, {𝑑 2 , 𝑑 4 }} = max(⊆, {{𝑑 1 , 𝑑 2 , 𝑑 3 } ∩ {𝑑 1 , 𝑑 2 }, {𝑑 1 , 𝑑 2 , 𝑑 3 } ∩ {𝑑 2 , 𝑑 4 }}) = max(⊆, {{𝑑 1 , 𝑑 2 }, {𝑑 2 }}) = {{𝑑 1 , 𝑑 2 }}
We only retain {𝑑 1 , 𝑑 2 } since {𝑑 2 } ⊆ {𝑑 1 , 𝑑 2 } and {𝑑 1 , 𝑑 2 } is the only ⊆-maximal element. Indeed, the semantic of {𝑑 2 } a prescription that contains the drug 𝑑 2is more general than the semantic of {𝑑 1 , 𝑑 2 }a prescription that contains both the drugs 𝑑 1 and 𝑑 2 .

Given that each patient has a description for each firstlevel ICD class, the meet operator defined for a single 3 description can be applied to a vector of sub-descriptions as follows:

𝛿 1 (P1) ⊓ 1 𝛿 1 (P2) = 〈𝛿 1,1 (P1), … , 𝛿 1,𝑛 (P1)〉 ⊓ 1 〈𝛿 1,1 (P2), … , 𝛿 1,𝑛 (P2)〉 = 〈𝛿 1,1 (P1) ⊓ 1 𝛿 1,1 (P2), … , 𝛿 1,𝑛 (P1) ⊓ 1 𝛿 1,𝑛 (P2)〉
Figure 1 shows the semi-lattice associated with this pattern structure and the data in Table 2. This example illustrates that information about the ADEs is quickly lost as we generalize. 2 using the pattern structure (𝐺, (𝒟 1 ,⊓ 1 ), 𝛿 1 ), where arrows denote the partial order ≤ ⊓ 1 .

Extending the Pattern Structure with a Drug Ontology

Using a drug ontology permits to find associations between ADEs related to classes of drugs rather than individual drugs. Thus we extend the pattern structure described previously to take into account a drug ontology: ATC. Each drug is mapped to its ATC class(es), as shown in Table 4. We define this second pattern structure (𝐺, (𝒟 2 ,⊓ 2 ), 𝛿 2 ) where descriptions of 𝒟 2 are sets of prescriptions with drugs represented as their ATC classes. In order to compare sets of classes from an ontology 𝒪, we define an intermediate meet operator ⊓ 𝒪 , for 𝑥 and 𝑦 any two sets of classes of 𝒟 :

𝑥 ⊓ 𝒪 𝑦 = max(⊑, {𝐿𝐶𝐴(𝑐 𝑥 , 𝑐 𝑦 ) | (𝑐 𝑥 , 𝑐 𝑦 ) ∈ 𝑥 × 𝑦})
where 𝐿𝐶𝐴(𝑐 𝑥 , 𝑐 𝑦 ) is the least common ancestor of 𝑐 𝑥 and 𝑐 𝑦 in 𝒪, and ⊑ is the ordering defined by the class hierarchy of 𝒪 . Thus max(⊑, 𝑆) is the subset of most specific ontology classes in 𝑆, and 𝑥 ⊓ 𝒪 𝑦 is the set of most specific ancestors of classes in 𝑥 and 𝑦 . From ⊓ 𝒪 we define the partial order ≤ 𝒪 which compares the semantics of a set of ontology classes, such that 𝑥 ≤ 𝒪 𝑦 ⇔ 𝑥 ⊓ 𝒪 𝑦 = 𝑥 and 𝑥 ≤ 𝒪 𝑦 denotes that 𝑦 is a more specific set of ontology classes than 𝑥 .

We then define the meet operator ⊓ 2 such that for every pair of descriptions (𝑋, 𝑌) of 𝒟 2 :

𝑋 ⊓ 2 𝑌 = max(≤ 𝒪 , {𝑥 ⊓ 𝒪 𝑦 | (𝑥, 𝑦) ∈ 𝑋 × 𝑌})
This pattern structure allows generalization of ADEs involving different drugs that share a pharmacological subgroup. For instance:

𝛿(P1) ⊓ 2 𝛿(P4) = 〈{{H02AB07}, {N02BE01}} ⊓ 2 {{H02AA03}} , ∅〉 = 〈{H02A} , ∅〉
This vector represents the closest generalization of the descriptions of patients P1 and P4, and can be read as: drugs of the class H02A (corticosteroids for systemic use, plain) are associated with a phenotype in the ICD class 580-629, no drugs associated to the ICD class 710-739.

Extending the Pattern Structure with a Phenotype Ontology

We define a third pattern structure that permits the use of both ATC and ICD for better specialization of phenotypes compared to the previous experiment. Here, we use only the two most-specific levels of ICD, as the previous experiment already covers the most general level. Table 5 describes the data representation used with this pattern structure. Here, ADEs are represented as vectors 〈𝐷 𝑖 , 𝑃 𝑖 〉 with two dimensions: the set of drugs 𝐷 𝑖 associated with the set of phenotypes 𝑃 𝑖 . A patient description is then a set of such vectors. We define the pattern structure (𝐺, (𝒟 3 ,⊓ 3 ), 𝛿 3 ), where descriptions of 𝒟 3 are sets of ADEs. We first define an intermediate meet operator ⊓ 𝐴𝐷𝐸 on our ADEs representations:

𝑣 𝑥 ⊓ 𝐴𝐷𝐸 𝑣 𝑦 = 〈𝑥 𝐴𝑇𝐶 , 𝑥 𝐼𝐶𝐷 〉 ⊓ 𝐴𝐷𝐸 〈𝑦 𝐴𝑇𝐶 , 𝑦 𝐼𝐶𝐷 〉 = { 〈𝑥 𝐴𝑇𝐶 ⊓ 𝒪 𝑦 𝐴𝑇𝐶 , 𝑥 𝐼𝐶𝐷 ⊓ 𝒪 𝑦 𝐼𝐶𝐷 〉 if both dimensions
contains non-root classes, 〈∅, ∅〉 otherwise.

The operator ⊓ 𝐴𝐷𝐸 applies the ontology meet operator ⊓ 𝒪 on both dimensions of the vector representing the ADE, using either ATC or ICD as 𝒪 . Both dimensions of the resulting vector needs to contain non-root ontology classes for it to constitute a representation of an ADE. If it is not the case, we set it to 〈∅, ∅〉 to ignore that ADE in further generalizations.

We define the meet operator ⊓ 3 such that for every pair of descriptions (𝑋, 𝑌) of 𝒟 3 :

𝑋 ⊓ 3 𝑌 = max(≤ 𝐴𝐷𝐸 , {𝑣 𝑥 ⊓ 𝐴𝐷𝐸 𝑣 𝑦 | (𝑣 𝑥 , 𝑣 𝑦 ) ∈ 𝑋 × 𝑌})
Compared to ⊓ 2 , ⊓ 3 adds a level of computation with ⊓ 𝐴𝐷𝐸 which generalizes ADEs, and applies ⊓ 𝒪 to an additional ontology: ICD9-CM.

RESULTS AND DISCUSSION

The three experiments result in three concept lattices, from which we extract Association Rules (ARs). An AR is identified between two related concepts in the lattice, with descriptions 𝛿(𝑙) and 𝛿(𝑟) such that 𝛿(𝑟) is more specific than 𝛿(𝑙). Thus, such an AR comprises a left part 𝐿 = 𝛿(𝑙) and a right part 𝑅 = 𝛿(𝑟)\𝛿(𝑙) , and is denoted 𝐿 → 𝑅 . Empirically, we only retain ARs with support (number of patients verifying the rule) at least 5, and confidence (ratio of patients verifying 𝐿 that also verify 𝑅 ) at least 0.75 . Table 6 presents a few statistics about this process in our three experiments. As expected, this process generates a large amount of rules, among which ARs serving our goal of identifying associations between ADEs must be identified. We therefore define a filter on ARs as the conjunction of the following conditions. (i) The right part 𝑅 of the AR contains at least one ADE, noted as 〈𝐷 𝑅 , 𝑃 𝑅 〉 for which there does not exist an ADE in the left part 𝐿, 〈𝐷 𝐿 , 𝑃 𝐿 〉, such that either 𝐷 𝑅 and 𝐷 𝐿 are ≤ 𝒪 comparable, or 𝑃 𝑅 and 𝑃 𝐿 are ≤ 𝒪 comparable. This condition ensures that the right part of the rule introduces new drugs and phenotypes unrelated to those of the left part, i.e., the association between ADEs is not trivial. (ii) As patients in the corpus are treated for Systemic Lupus Erythematosus rules must not include related phenotypes (ICD class 710 and descendants).

Finally, we present an example of AR obtained in the third pattern structure experiment, with support 10 and confidence 0.77: {〈{C08DB01}, {ICD 428.0}〉} → {〈{A02B}, {ICD 427.31}〉}

This rule means that 77% of patients that present congestive heart failure (ICD 428.0) after prescription of diltiazem ( C08DB01 ), also present atrial fibrillation (ICD 427.31) after prescription of a drug for peptic ulcer and gastroesophageal reflux disease (A02B). This rule holds true for 10 patients in our corpus. The complete set of filtered rules for each of the three experiments is available online (http://www.loria.fr/~gpersone/ade-assoc/).

In summary, we explore in this paper an approach based on pattern structures to mine of EHRs for commonly associated ADEs. It permits a detailed representation of ADEs, enriched with medical ontologies. This representation could be further enriched with additional properties of drugs and phenotypes, such as drugs targets annotated with Gene Ontology classes.

A large amount of ARs can be extracted from our concept lattices. We automatically selected a subset of these ARs by excluding rules that do not fit the scope of the study. We now need to prioritize these ARs with respect to their importance in terms of cost and risk of the phenotypes present in their right part.

One limitation of this work is the absence of temporal relationships between ADEs. We did not consider that aspect because the order of occurrence of ADEs can vary between patients. However, in cases of interest, this order can be checked in patients EHRs as pattern structures concepts retain patient identifiers as well as their description.

Figure 1 .

 1 Figure 1. Semi-lattice representation of the data in Table2using the pattern structure (𝐺, (𝒟 1 ,⊓ 1 ), 𝛿 1 ), where arrows denote the partial order ≤ ⊓ 1 .

Table 1 .

 1 Number of patients with at least 2 selected ADEs and number of ADEs for these patients, for different maximum inter-visit interval in days.

	Interval (days)	1	2	6	10	14	18	22	26	30
	|Patients|	434 461 498 526 548 555 558 564 576
	|ADEs|	2,396 2,587 2,902 3,110 3,286 3,388 3,454 3,501 3,621

Table 2 .

 2 Example of a corpus containing 3 patients and 6 ADEs.

	Patient	ADEs
	P1	({prednisone}, {ICD 599.8}) ({acetaminophen}, {ICD 599.9})
	P2	({prednisone}, {ICD 599.8}) ({prednisone}, {ICD 719.4})
	P3	({prednisone, acetaminophen}, {ICD 599.9}) ({acetaminophen}, {ICD 719.4})

Table 3 .

 3 Example of representation of patient ADEs for (𝐺, (𝒟 1 ,⊓ 1 ), 𝛿 1 ), with two first-level ICD classes: diseases of the genitourinary system (580-629), and of the musculoskeletal system and connective tissue (710-739).

		ICD 580-629	ICD 710-739
	Patient P1	{{prednisone}, {acetaminophen}}	∅
	Patient P2	{{prednisone}}	{{prednisone}}
	Patient P3	{{prednisone, acetaminophen}}	{{acetaminophen}}

Table 4 .

 4 Example of representation of patient ADEs for (𝐺, (𝒟 2 ,⊓ 2 ), 𝛿 2 ).

	Class labels: H02AA03 is desoxycortone, H02AB07 is prednisone,
	N02BE01 is acetaminophen.
		ICD 580-629	ICD 710-739
	Patient P1	{{H02AB07}, {N02BE01}}	∅
	Patient P2	{{H02AB07}}	{{H02AB07}}
	Patient P3	{{H02AB07, N02BE01}}	{{N02BE01}}
	Patient P4	{{H02AA03}}	∅

Table 5 .

 5 Example of representation of patient ADEs for (𝐺, (𝒟 3 ,⊓ 3 ), 𝛿 3 )

	Description

Table 6 .

 6 Statistics about the lattice building and AR extraction processes implemented in Java.

	Experiment	1	2	3
	Lattice size (millions of concepts)	1.9	2.3	2.5
	ARs extracted (millions)	5	7	9
	ARs retained after filtering	772 1,907 913