Nicolas Bonneel

Balazs Kovacs

Sylvain Paris

Kavita Bala

STAR -State of The Art

Keywords: Descriptors (according to ACM CCS): I, 3, 3 [Computer Graphics]: Picture/Image Generation-Line and curve generation

Figure 1: We evaluate state-of-the-art intrinsic image decomposition algorithms based on their ability to produce seamless, artifact-free results for image edits. To fairly compare different methods, we automatize the image editing process. Left to right, by Poisson-based inpainting of the reflectance layer, we remove a logo on a shirt using the metod of Barron et al. [BM15], we add a picnic blanket over a shadow with the method of Grosse et al. [GJAF09], and add a painting over colored shadows with the method of Bousseau et al. [BPD09].

Introduction

The rich visual world that surrounds us is the result of the complex interplay between light and matter. Light reaches an observer typically after several interactions with physical objects, each of them having a nontrivial effect on its spectrum, different light wavelengths being affected differently depending on the properties of the material involved. Further, cause and effect can be separated by a large distance when an object casts a shadow far away from its actual position. Eventually, all these effects are conflated into a single image. Considered in its full generality, the image formation process may seem impossible to invert because the phenomena involved are too diverse and too complex to be disentangled. Yet, human observers effortlessly identify shadows and recognize object colors under all but the most extreme lighting conditions. Intrinsic decomposition of digital images is a task inspired from this ability. Originally, Barrow and Tenenbaum [BT78] sought to characterize properties inherent in a scene such as the color and geometry of objects independent of viewing conditions. With time, the geometric aspects of this original goal became associated with stereo and multi-view reconstruction, and intrinsic decomposition took a more focused meaning becoming synonymous with reflectance estimation. In this manuscript, we follow this trend and define intrinsic decomposition as the task of separating the effect of the scene illumination from that of the material reflectance. This task relies on a simple image formation model that explains each pixel as the product of two RGB triplets, one for the light color and one for material reflectance. This model analyses only the last bounce of the light transport and makes a number of simplifying assumptions such as "all materials are Lambertian" and "no participating media", but it nonetheless provides a powerful means to reason about light and object colors that has proven to be useful to many image editing applications. For instance, recoloring an object in an image is a nontrivial task even for a uniform-color object because shading variations make some parts bright and others dark, possibly in a discontinuous way if the object exhibits sharp geometric features. The same task with an intrinsic decomposition is straightforward since the reflectance of the object is constant. The rest of this manuscript presents the concept of intrinsic image decomposition in more detail, describes the main existing algorithms to compute such a decomposition, and reviews the most common use cases.

While some approaches consider multiple images as input (from multiple viewpoints [CBLD11, LBP * 12, Laf12, LBD13, HWU * 14, Duc15, XLL * 16], varying illumination [Wei01, [START_REF] Matsushita | Estimating intrinsic images from image sequences with biased illumination[END_REF][START_REF]Rank-constrained pca for intrinsic images decomposition[END_REF], or different focal distances [START_REF] Saini | Intrinsic image decomposition using focal stacks[END_REF]), images with depth information [BM13, CK13, JCTL14, HGW15], multi-spectral images [START_REF] Shao | Extracting intrinsic images from multispectral[END_REF], videos [YGL * 14, BST * 14, KGB14, SYC * 14, MZRT16], lightfields [GEZ * 16, AG16] or even videos with depth [LZT * 12], this document focuses on the use of a single RGB image as input, and emphasizes image editing applications. This is motivated by the wide availability of this kind of data and the need for illumination-aware image editing tools.

Aside from pedagogical content, this document makes the following contributions. First, we evaluate priors commonly used in the literature in the context of realistic and complex scenes. Second, we introduce a small but realistic synthetic ground-truth dataset based on PBRT [PH10], Mitsuba [START_REF]Mitsuba renderer[END_REF] and LuxRender [START_REF] Vergauwen | Luxrender -gpl physically based renderer[END_REF] scenes. Third, we evaluate 10 recent methods on real photographs in the context of image editing applications.

Problem formulation

The intrinsic decomposition problem seeks to decompose an image into the product of illumination and reflectance layers. This section exposes the motivation behind this problem, as well as various assumptions and priors that numerically help solve for this decomposition.

Intrinsic Decomposition

In a Lambertian scene, material reflectance does not depend on viewing direction or illumination incidence. This simplifies light transport and allows for the writing of a simplified (yet exact) physicallybased image formation model. In this context, we have

I(x, λ) = ρ(x, λ) L(x, λ) (1)
where x is the pixel position, λ the light wavelength, I the rendered image, ρ the diffuse albedo and L is a term which depends on light and geometry. In the following, we will call ρ the reflectance, and L the illumination.

This model holds in the continuous image plane domain, but the spatial filtering and sampling that occur inside real cameras makes the equality in the equation above break on traditional pixel grids when geometry, textures or illumination vary within a single pixel, or in the presence of lens or motion blur. Cameras also often store non-linearly processed pixel values, for instance, to gamma correct, enhance or white balance images.

The inverse problem of forming an image using the model in Eq. 1 is the problem of recovering the reflectance and illumination given an already formed image. This process is called intrinsic decomposition. Eq. 1 makes the intrinsic decomposition problem precisely defined in terms of photometric quantities.

Other techniques also try to understand the role of lighting or textures in images, or relate to intrinsic decompositions: A major challenge of intrinsic image decomposition is that the image formation model I = ρL is ill-posed because if ρ 0 and L 0 satisfy the model, i.e., I = ρ 0 L 0 , then aρ 0 and 1 a L 0 also satisfy the model for any nonzero a, including the case where a is spatially varying (Fig. 2). Concretely, this means that the absolute ground-truth decomposition is unattainable unless additional absolute measurements are available, e.g., using a light meter, which is not the case This also enables the quantitative evalution of the accuracy. For instance, one can use computer-generated images for which ground truth is available as a benchmark. Instead of directly reporting the error between a given decomposition and the ground truth as is often the case in other contexts, one first searches for the multiplicative factor that minimizes the error, which is the value reported. We will discuss further how to evaluate results later in the manuscript. However, even considering relative accuracy, the image formation model leaves room for too much freedom: the model still holds when a is spatially varying. In fact, twice as many unknown values (ρ and L) as known quantities (I) have to be estimated. To attain satisfactory decompositions, priors and constraints are needed to reduce the number of unknown variables and disambiguate the problem. This represents the heart of intrinsic decomposition research, and is discussed next.

Common priors

Priors statistically model one's beliefs about intrinsic decompositions and help disambiguating decompositions. A number of priors have been introduced in the literature, as well as assumptions and user constraints, detailed below and summarized in table 8. In the following, we assume the camera sensor response curve has been properly taken care of (for instance, by directly working on raw images or using photometric calibration techniques [START_REF] Krawczyk | Photometric Calibration of High Dynamic Range Cameras[END_REF]).

Monochromatic illumination (MI).

Often, the illumination layer is assumed to be grayscale. Up to a white balancing step of the input photograph, this corresponds to the use of a single light color and reduces the illumination layer to a single scalar value per pixel L(x) instead of L(x, λ). This is the most common assumption (see Table 8), and only few approaches allow for colored lights [CSBC09, BPD09, CCFI14, BM15, TNY15] ( [BM13] for RGB-D images). It is interesting to note that the monochromatic assumption often only applies to the light color as seen by the camera. That is, the light itself need not be monochromatic as long as the integral of its energy distribution over wavelengths is the same for each color sensor. The main disadvantage of this assumption is that inter-reflections are most often colored [CRA11], and so, this assumption often fails to capture illumination effects that are due to indirect lighting. However, the monochromatic illumination assumption can be considered as a prior by designing a soft penalty for strongly colored illuminations while still allowing for colored lighting. Edge is either a reflectance or illumination change (EoI). This states that it is unlikely for a single image edge to result from both a change of illumination and reflectance at the same time [LM71]. This has been used to classify image gradients [Hor74, GJAF09, LSX09, TFA05]. However, this assumption is unlikely to hold at occlusion boundaries: in this case, a silhouette edge separates an occluding surface whose normal is parallel to the image plane (assuming orthographic projection) [START_REF] Barron | Shape, illumination, and reflectance from shading[END_REF] from an occluded surface whose normal can be arbitrary. This change in normal direction often coincides with a change in illumination, and in many cases, these two surfaces will have different reflectances as well. This prior can be seen as a strong version of the Retinex prior.

Chrominance variations are more likely to be reflectance variations (CR). When chrominance varies abruptly, this is more likely due to a change of reflectance than illumination [GJAF09, TFA05, ZTD * 12].This can be implemented as a varying gradient threshold for abrupt changes in chrominance [START_REF] Tappen | Recovering intrinsic images from a single image[END_REF][START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF], by weighting chrominance and luminance differently in the pairwise cost functions used in CRF-based frameworks [BBS14, ZKE15], or by introducing a quadratic penalty term which enforces smooth reflectance only in places where the chrominance is smooth [START_REF] Shen | Intrinsic images decomposition using a local and global sparse representation of reflectance[END_REF]. A similar observation concludes that if hue variations correlate positively with intensity variations, this is more likely due to changes in reflectance [JSW10].

Low rank reflectances (LRR). This assumes that locally, within a small neighborhood, reflectance values form a 2-d affine subspace of the RGB space [START_REF] Bousseau | User assisted intrinsic images[END_REF]. This prior can be seen as a reflectance smoothness prior.

Sparse reflectance values (SRV)

. This assumes that most color variations are due to illumination, and that, in fact, few different reflectance values make up a typical image. This assumption is expected to work best for photographs of man-made scenes. This is often implemented via a color clustering step [GMLMG12, BBS14, LYZ15], or a sparsity constraint on reflectance [START_REF] Shen | Intrinsic images decomposition using a local and global sparse representation of reflectance[END_REF], or even using superpixels [BHY15]. Alternatively, this prior can be cast in the realm of information theory. Some reflectance values are more likely (RML). Barron and Malik [START_REF] Barron | Shape, illumination, and reflectance from shading[END_REF] assumes some reflectance values are more likely than others. This is implemented by building a 1-d histogram of logreflectance values of a ground truth dataset, and using it as a prior. This prior helps disambiguate the overall light color from reflectance colors: a blue pixel will be more likely the result of a white light illuminating a blue object, rather than a magenta light illuminating a cyan object, if cyan reflectances are a priori less likely than blue reflectances.

Mean correlates with variance (MV). Under illumination variations, the local mean of pixel values within a neighborhood should vary in the same direction as the local variance [JSW10]

. This is due to the fact that illumination acts multiplicatively on reflectance, which can be detected via correlation analysis.

Planckian lighting (PL). Under skylight and a narrow-band camera sensor, it can be shown that pixels belonging to objects of the same reflectance form a single line in log-RGB space when varying the lighting condition [CPCN13, FDL04, LYZ15]. This is often implemented by clustering lines of log-RGB pixel values, similar to the color lines model [START_REF] Werman | Color lines: image specific color representation[END_REF] (though not performed in the log-domain), or via entropy minimization [FDL04].

Non-local constraints (NLC). The goal of this prior is to find pixels that are most likely of the same reflectance value within an image. The idea is to compare texture descriptors, such as pixel neighborhoods, and if two descriptors agree, they most likely represent the same structure repeated at a different location [STL08, ZTD * 12]. This introduces long-range reflectance constraints. The difficulty lies within the comparison function, which assumes neighbors can be compared in a way agnostic to illumination variations... a chickenand-egg problem! In practice, simply dividing pixel color values by their intensity (i.e., taking the pixel chromaticity) often serves as a good proxy [STL08, ZTD * 12].

User-defined constraints (UC). The difficulty of automatic intrinsic decomposition has led researchers to rely on the user to add constraints. This typically involves asking the user to mark pairs of pixels of similar reflectance or illumination, or brush areas of similar

Input image

Reflectance masks log-RGB scatterplot This hardly results in distinct lines. This prior may produce artifacts when used on images taken in uncontrolled settings.

reflectance or brush absolute illumination values [START_REF] Bousseau | User assisted intrinsic images[END_REF]. These constraints make the decomposition interactive instead of fully automatic. While this requires effort from users, this often yields better decompositions. This is simply accomplished by adding these constraints in a linear system [BPD09, BST * 14] or similarly adding a quadratic penalty term [START_REF] Shen | Intrinsic images using optimization[END_REF]. While this state-of-the-art report focuses on single image intrinsic decompositions for image editing applications, we will briefly mention other specific priors and assumptions that have been introduced in the context of more general intrinsic decompositions.

Numerical techniques

Various numerical techniques have been investigated to account for (part of) these priors, some of which were described in Sec. 2.2. The mathematical formulation of these priors and of the image formation model matters in practice and is key to designing practical algorithms. We next review a few standard approaches.

A log-transform often conveniently rewrites the product I = ρ × L into a sum log I = log ρ + log L [LM71, Hor74]. Changing the name of these variables (here, lowercase letters denoting log values), this can be written as i = r + l, simplifying a non-linear to a linear relationship. Numerical tools from linear algebra can then be used, particularly when priors can also be expressed as linear relationships (for example, as the result of the minimization of a quadratic energy). This often leads to sparse linear systems [START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF]. Interestingly, log-transforms also often render these methods robust to gamma correction as log I γ = γ log I. Priors evoking smoothness or sparsity can often be expressed using gradients -for instance, methods based on the Retinex theory may classify image gradients as belonging to either the illumination layer or the reflectance layer [START_REF] Tappen | Recovering intrinsic images from a single image[END_REF]. In conjunction with the log-transform, the intrinsic decomposition can be advantageously rewritten as

∇ log I = ∇ log ρ + ∇ log L.
When priors cannot be easily cast as linear constraints (or quadratic penalty terms), full non-linear solvers have been used, such as l-BFGS [LZL14, BM15].

Alternatively, a probabilistic approach can be taken via CRFs [BBS14]. Here, priors are expressed via probabilistic models as exponentially decreasing functions of some energies, whose joint negative log-likelihood is minimized for. In traditional CRFs, message passing is used to minimize the energy function, though Bell et al.

[BBS14] approximate it via high-dimensional filtering.

Finally, with the recent advances in machine learning and the availability of ground-truth datasets, learning-based approaches have emerged. These methods guide the decomposition by directly classifying image gradients that are propagated using a Markov Random Field (MRF) [START_REF] Tappen | Recovering intrinsic images from a single image[END_REF], more generally regress image filters that are propagated via a pseudo-inverse [START_REF] Tappen | Recovering intrinsic images from a single image[END_REF], or use full-fledged convolutional neural networks [START_REF] Takuya | Direct intrinsics: Learning albedo-shading decomposition by convolutional regression[END_REF] or deep belief networks [START_REF] Tang | Deep Lambertian Networks[END_REF]. Jointly learning depth and intrinsic decomposition via deep convolutional network has seen some success [SBD15], and joint estimation of shape, illumination (as a spherical low-frequency environment map) and intrinsic decomposition performs well [START_REF] Barron | Shape, illumination, and reflectance from shading[END_REF]. Recently, multiple works emerged which learn a reflectance prior from the pairwise judgments of Bell et al.

Typical applications

Decomposing an image into illumination and reflectance components has several applications. First, this allows for understanding scenes better. For instance, it could be used to understand the role of the illumination with respect to intrinsic reflectance color in the popular blue-black dress meme [LSHC15] (see Fig. 9). Although the actual dress is blue and black, 30% of people perceive this dress as white and gold, due to a different perception of the illuminant [LSHC15] and other biological factors [START_REF] Rabin | Blue-black or white-gold? early stage processing and the color of 'the dress[END_REF].

Input

Reflectance Illumination Reflectance Illumination Editing results are occasionally illustrated on textureless surfaces, e.g., to alter the color of a uniform object. While it is un-challenging, it is also primarily more easily performed via simple luminancechrominance adjustments! In fact, a luminance-chrominance decomposition is a correct and valid intrinsic decomposition for textureless uniform surfaces. We hence recommend comparisons against naive baselines, as they often perform reasonably well for many image editing applications.

Evaluation and Comparisons

Evaluating intrinsic decomposition methods is not a trivial task. The seminal paper of Grosse et al. [START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF] introduced the ground-truth dataset now known as the MIT dataset. This paper advocated for the LMSE metric, a mean-square error computed and averaged locally, within pixel neighborhoods. Bell et al. [BBS14] relied on perceptual experiments to determine how two pixels differ in their reflectance, and compare these judgments with the results of automatic algorithms, leading to the WHDR, weighted human disagreement rate, metric. Instead, in this paper we evaluate intrinsic decomposition methods based on their applicative ability. That is, we do not need a decomposition to be accurate, as long as it is sufficient to perform a given image editing task. We thus only evaluate the result of the image editing process. This section describes this approach, as well as a more classical evaluation using LMSE on a new ground truth realistic synthetic dataset.

Datasets

Ground truth datasets are important for evaluating intrinsic decomposition methods, but also to provide data for training approaches based on machine learning. Datasets often precisely meet the assumptions of methods being assessed. This report instead focuses on realistic image editing applications, and we believe slightly violating assumptions is reasonable in this context.

In particular, the MIT dataset [START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF] relies on isolated objects with black background to minimize interreflections, and a single directional light source. The reflectance component is obtained by coating the object with a white diffuse paint, while specularities are removed via a polarizing filter. This technique does not allow for colored indirect illumination as the white coating may not reflect the same colored light as the initial object. Sierra used this same technique to extend this intrinsic decomposition database [Ser15]. The MIT dataset has been widely used for benchmarking intrinsic decomposition methods, but has been deemed "not representative of the variety of natural objects in the world" [START_REF] Barron | Shape, illumination, and reflectance from shading[END_REF], and would be hardly useful for assessing methods in ecological contexts. Beigpour et al. have more recently extended this dataset [BKK15, BHK * 16] under multiple lighting and viewing conditions using the procedure of Grosse et al. [START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF]. The two datasets each contain 5 scenes under 17 illumination conditions, and contain ground-truth depth and specularity information. They consist of two objects resting on a planar surface and remain of moderate complexity. The second dataset features 6 view conditions. A similar approach taken by Tappen et al. [START_REF] Tappen | Estimating intrinsic component images using non-linear regression[END_REF] uses colors to capture reflectance and illumination independently. In practice, they color a piece of paper using a green pen: this color is invisible in the red channel of the image, and the red channel is taken as the ground-truth illumination layer. They use this technique to build a ground truth database of isolated sheets of crumpled paper.

Bell et al. introduced an extensive crowdsourced dataset of pairwise reflectance comparisons for photos 'in the wild' (i.e., for Flickr images taken from real world settings) [BBS14]. MTurk workers were asked to determine whether random pairs of points share the same reflectance, or if one point has a darker surface color. The IIW dataset comprising 5,230 photos, includes 875,833 reflectance comparisons, and currently is the largest database for benchmarking intrinsic image decomposition algorithms. However, despite its size, since the dataset provides ground-truth data for sparse pairs of points, it is not amenable to evaluating high-frequency errors that often occur in the illumination component, which can lead to artifacts in re-texturing applications. * 13] introduced a synthetic dataset by rendering 3d scenes under various lighting conditions -including Planckian lighting. They provide renderings for 8 isolated objects and 9 complex scenes. While these complex renderings are indeed more complex than the MIT dataset, they consist uniquely of textureless objects. Other scattered datasets exist in specific contexts. The MPI-Sintel dataset offers 23 rendered video sequences used for assessing optical flow methods [BWSB12]. In addition to optical flow ground truth information, this dataset also offers ground truth reflectance and depth (among other data), and has thus been used to assess intrinsic decomposition methods [START_REF] Takuya | Direct intrinsics: Learning albedo-shading decomposition by convolutional regression[END_REF]. Fig. 10 illustrates examples from all these datasets. This report instead provides a dataset of 18 reflectance images from photorealistic scenes from the PBRT [PH10], Mitsuba [START_REF]Mitsuba renderer[END_REF] and LuxRender [START_REF] Vergauwen | Luxrender -gpl physically based renderer[END_REF] renderers (see Fig. 5, 11 and supplemental material). The illumination component is then recovered by taking the ratio between the original image and the reflectance layer. However, on realistic images, several of our assumptions break, due to specular components, transmissive surfaces, subsurface scattering and defocus blur. As such, the illumination component may contain artifacts, such as residual reflectance. In addition, artists often simplify complex but small-scale geometries using flat textures on coarser objects. This occasionally yields overly smooth illumination layers, as light transport is not properly simulated on this geometry. They may also add high contrast materials that do not reflect the real world measured materials.

Beigpour et al. [BSV

In the context of image editing applications, we further assess recent intrinsic decomposition techniques on a set of 21 real photographs (see Sec. 4.3).

Ground-truth comparisons

With 9 images from our new realistic ground truth dataset, we first evaluate the results of various algorithms using a classical LMSE metric. For automatic methods, we experimented with several parameter sets (up to 24, for the color Retinex of Grosse et al. [START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF]) and for each image, we kept the result minimizing the LMSE. For interactive methods, we manually adjusted parameters interactively and added strokes to visually obtain the best result possible. Fig. 12 plots the LMSE of tested intrinsic decomposition techniques using box-and-whisker plots, and sorts these methods by decreasing average LMSE. In terms of LMSE, the method of Shen et al. [START_REF] Shen | Intrinsic images using optimization[END_REF] performs best on this benchmark. However, as we shall see in Sec. 4.4, this does not portray an accurate picture of the state-of-the-art in intrinsic image decomposition when one focuses on image editing applications. In fact, piece-wise accurate results -such as faithful reproduction of large areas of constant reflectance -typically yield low LMSE values, while computer graphics applications are less tolerant to localized mistakes. While specularities occasionally yield residual textures in the ground truth shading, we found that this had little impact in the computed LMSE. For instance, on the "Breakfast" scene (Fig. 5), manually correcting the shading layer changes the LMSE by approximately 0.1%.

Evaluation on image-editing applications

We assess several state-of-the-art methods on image editing applications. We use a database of 21 Creative Commons photographs downloaded from Flickr, spanning different but realistic contexts, and exhibiting interesting illumination and reflectance variations. These include portraits, interior and exterior scenes. For all of these images, we have determined a specific image editing task that an artist could perform via intrinsic decomposition. These include removing a logo on a t-shirt, a tattoo or makeup, smoothing out wrinkles or freckles, or altering a shadow.

To assess multiple intrinsic decomposition methods efficiently, and to fairly evaluate methods with the same image editing operation, we automate these edits. Particular care has been taken to account for the varied dynamic ranges and the intrinsic decomposition scale

Input image

Reflectance Illumination invariance. For instance, replacing a texture should not be performed by directly editing absolute reflectance values, since they may differ from one decomposition to another. Instead, we favored gradientdomain approaches, or filtering operations on the different layers.

Logo removal. We manually determine a rough mask for a logo to be removed, and solve for the Poisson problem ∆u = 0 within the masked domain and u = ρ outside, with Dirichlet boundary conditions. We use the solution of this problem as the modified reflectance ρ , and reconstruct the final image as I = ρ × L. This effectively removes the logo on successful intrinsic decomposition results (see Fig. 1).

Shadow removal. We apply the same process as for logo removal, but apply it to the illumination layer (see Fig. 16, fourth row of results, right column). Texture replacement. We inpaint a new texture in the reflectance layer by using Poisson Image Editing [PGB03]. Since complex high-frequency textures may hide artifacts in the processed result -a phenomenon called spatial frequency masking [Dal93] -we choose textures of relatively low frequency content. Further, textures may contain low-frequency illumination variations that do not correlate with the scene geometry. We hence high-pass textures containing residual low-frequency illumination. We avoid complex light-geometry interactions by integrating mostly planar objects on planar surfaces, such as carpets or paintings.

Input image Reflectance Illumination

Wrinkles attenuation. Wrinkles are mostly due to shadowing effects on the skin. We manually determine a rough mask for the skin area to be corrected. We blur both the mask and the illumination layer via Gaussian filtering. We linearly interpolate the input illumination with the altered illumination based on the blurred mask, and obtain the final illumination layer L . We reconstruct the final image as I = ρ × L .

We provide the code and data in supplemental material for benchmarking purposes.

Evaluation

With our set of automatically generated image-processed results for various intrinsic decomposition methods, this section evaluates their success. We deem a method successful if both of the following criteria are met:

• The effect has been achieved. That is, if the goal is to remove a logo, the final image should not contain the logo anymore. Indeed, it is easy for a method to be free of any visible artifact, but to miss its primary purpose (e.g., a luminance-chrominance decomposition). • The result is realistic. That is, images do not contain artifacts, and given a processed and unprocessed images, one cannot determine which one is processed. The method should thus not deteriorate the quality of the final result. Note that this criterion also depends on the realism of the image editing process. We have hence put significant research effort in minimizing artifacts that are inherently due to the automatic image editing process.

Materials. We decompose our dataset of 21 images with the methods of [GJAF09, BPD09, GRK * 11, SYJL11, GMLMG12, ZTD * 12, BST * 14, BBS14, ZKE15, BM15, TNY15], using the author implementations. For the method of Grosse et al. [START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF], we evaluate both the grayscale and color Retinex approaches, and keep the best performing result between an L 1 and L 2 gradient reconstruction. For the method of Bonneel et al. [BST * 14], we evaluate both the automatic and user-assisted approaches. For each algorithm, we downsample the input images to obtain reasonable computation time, memory usage and robustness, to the largest size the algorithm could handle. We experimented with multiple parameter sets, and kept the best result for each image. We additionally compute baseline decompositions as: 1) "Baseline reflectance" where the decomposed reflectance image is the chromaticity image, i.e. for an input pixel with RGB values (r, g, b), the output reflectance is r r+g+b , g r+g+b , b r+g+b . 2) "Baseline illumination" where the decomposed illumination image is constant 1. 3) "Baseline sqrt" where the decomposed illumination image is the square root of the grayscale image, i.e. for an input pixel with RGB values (r, g, b), the output grayscale illumination is r+g+b 3 . Our supplemental materials contain the downsampled input for each algorithm, as well as their intrinsic decomposition using the best parameter set.

Results. We note that regarding reflectance editing, most methods fail at completely removing textures from the illumination layer. This results in visible artifacts when removing a logo from a tshirt or inpainting an object in a photograph (see Fig. 14 and15). The method of Barron and Malik [START_REF] Barron | Shape, illumination, and reflectance from shading[END_REF] succeeds on few, but difficult, examples (see Fig. 1), and the color Retinex of Grosse et al. [START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF] works better on average but very rarely removes textures completely. The user-assisted approach of Bonneel et al. [BST * 14] sometimes succeeds but, conversely, tends to leave too much illumination in the reflectance layer. No method succeeds in removing a tattoo better than the best baseline decomposition, but we note that most methods have difficulties dealing with dark gray or black pigments. In fact, when removing textures, a simple gradient-domain inpainting of the input image often produces better results than intrinsic decomposition methods, as no residual texture pollutes the edited image and some illumination information is propagated from the mask boundaries.

The user-assisted method of Bousseau et al. [START_REF] Bousseau | User assisted intrinsic images[END_REF] is the only one to succeed in completely removing a strong cast shadow, even on simple geometries. This method handles colored illumination layers, which partly explains this success in addition to user cues. Regarding wrinkles removal, most approaches work reasonably well, but the baseline also succeeds in this case. The method of Narihira et al. [START_REF] Takuya | Direct intrinsics: Learning albedo-shading decomposition by convolutional regression[END_REF] does not produce a decomposition such that the product of the reflectance and illumination layers yields the input image, which causes significant artifacts and results in image edits that are worse than the baseline in all cases.

Our supplemental materials provides a combined view of all results. Fig. 16 shows all image edits obtained by the best performing user-assisted and automatic methods.

Other considerations

As priors are introduced, as well as heavier optimization routines, the speed of intrinsic decomposition techniques rarely meet realtime constraints. In fact, most methods require minutes and even sometimes hours to compute, even on low resolution (<1 mega-pixel) images. In most cases, we downsampled our test images to about 2 mega-pixels for this reason. However, most DSLR cameras now output photographs of tens of megapixels (e.g., Canon EOS series range from 18 to 50 mega-pixels), and even compact and phone cameras often come close to DSLRs in term of pixel resolution (e.g., Samsung Galaxy S7 Edge is 17 mega-pixels or Sony xperia z5 is 23 mega-pixels). With this amount of data, intrinsic decomposition techniques should be able to treat more than half-mega pixel images to be useful for image processing (though they could remain useful for vision applications or image understanding). Notable exception include the GPU framework of Meka et al. [MZRT16] that runs in a fraction of a second. Fig. 13 illustrates the running time with respect to image resolution for the tested algorithms. We did not time interactive methods for which most time is spent in user interactions. In practice, aside from user interactions, the method of Bousseau et al. [START_REF] Bousseau | User assisted intrinsic images[END_REF] takes between 5 and 30 seconds to solve for 0.5 to 1.8 mega-pixel images, and the method of Bonneel et al. [BTS * 15], initially designed for videos, approximately takes between 0.2 and 1 second for 1 to 2 mega-pixel images. We do not claim fair nor accurate times, as all the methods we have tested have been run on different machines, various implementations may differ and may not have been optimized for speed, and some methods make use of multi-threading or GPU, but our comparison gives a rough sense of computation times.

We believe speed is an important factor, as slow methods preclude interactive editing applications, fine parameter tuning, their adoption by artists, processing on large data such as frames of a video, and realtime vision applications. This could sometimes just be a matter of engineering and fine implementation tuning.

While certain intrinsic decomposition methods are more robust when working on high-dynamic range (HDR) images, and images with linear camera response, most images available on the web are not HDR images. For instance, even the widely used Flickr image search engine for photograph enthusiasts do not support images of more than 8 bit per pixel and color channel. We have thus evaluated existing methods on non-HDR images. However, HDR images are now more accessible, and intrinsic decomposition methods will likely benefit from this trend.

Finally, we have evaluated algorithms based on available implementations. We advocate for reproducible research and encourage authors to disseminate their code in addition to their research paper.

Challenges and future work

In the context of image editing, the quality of most intrinsic decomposition methods is not satisfactory as much reflectance is left in the illumination layer or vice versa, and in many cases their speed only make them useful for less than mega-pixel images. Intrinsic decompositions have met with limited success even when comparing to synthetic benchmarks. Such ground truth data is incredibly useful for validation and machine learning purposes. In the past, the lack of widely available realistic 3d scenes have constrained researchers to use very simple and unrealistic 3d scenes or isolated-object benchmarks. This may have biased machine learning approaches, and led us to think intrinsic decomposition is a near-solved problem. Recent datasets like IIW [BBS14] have expanded the evaluation to more real world setting using photographs, showing that more progress is required. In fact, in our experiments with 8-bit images, for both synthetic ground-truth comparison and image editing tasks, few methods perform better than a well-chosen baseline. The first challenge is to make intrinsic decompositions suitable for image editing applications, for images of reasonable size. We believe some widely used priors may be harmful for image editing, such as monochromatic illumination that is actually rarely observed in real-world scenes because of effects such as interreflection. Another difficulty with intrinsic decomposition is that it is an intermediate step as far as practical applications are concerned. This makes it difficult to rely on user-driven approaches in practice because they ask users to work on a task that is not the one in which they are interested. Ideally user intervention should not be needed or be as limited as possible, which raises its own challenges. On the other hand, fully automatic techniques are not yet accurate and/or fast enough to be used reliably on real-world high-resolution images.

Current quality metrics do not consider the application for these intrinsic decompositions. Direct comparison with a reference is interesting for many applications, such as intrinsically-guided segmentation or optical flow computation, but it may not be appropriate for image editing purpose. In fact, in this context small local errors degrade the perceived quality of the image edits much more than large low-frequency errors. As such, perceptual metrics are clearly necessary. While metrics like WHDR take a first step towards this goal by evaluating intrinsic decomposition results, more efforts are needed to focus on image editing applications.

Finally, the rise of augmented sensors -whether with additional depth information, lenslets to capture lightfields, or multi-spectral sensors -could alleviate the quality problem. While several approaches use depth information, very few approaches deal with lighfields (to our knowledge, only Garces et al. and Alperovich et al. [GEZ * 16, AG16] address static lightfield images) and multispectral sensors (to our knowledge, only Shao and Wang [START_REF] Shao | Extracting intrinsic images from multispectral[END_REF]). We further argue that intrinsic videos can now reasonably be handled via per-frame image intrinsic decompositions followed by temporal regularization [BTS * 15]. We believe that benefiting from multiple input modalities could be a promising direction to address our third challenge: the handling of more complex materials, such as transparent and glossy materials. To allow for better evaluation, comparisons and machine learning on complex, photorealistic renderings with complex materials, we are further extending our ground-truth dataset with other 3d renderings. In this extended dataset, we include multiple lighting conditions, normals, depth, position, irradiance and segmentations. This dataset is currently available at http://liris.cnrs.fr/~nbonneel/ intrinsicstar/ground_truth/.

Conclusions. Intrinsic image decomposition is a long standing problem of great importance to computer graphics and vision algorithms. Despite the inherent difficulty of solving this ill-posed problem, significant progress has been made through sophisticated algorithms and detailed datasets. More recently, learning from larger scale datasets is being used to solve the problem. However, when evaluated with the view of specific applications, like image-editing, additional research effort is needed to make intrinsic decomposition sufficiently accurate. We find that user input is often needed to achieve the best decompositions for image-editing applications. We also find that surprisingly simple baselines sometimes can be effective. We introduce new ground truth synthetic datasets, and advocate for the development of perceptual metrics, and more public datasets and algorithms to solve this important and challenging problem.

? ? × × × × × × × × ∼ × × × R ∼ × ∼ ∼ ∼ ∼ ∼ ∼ ∼ ∼ × EoI × ∼ CR ∼ ∼ ∼ ∼ ∼ ∼ LRR ∼ SRV ∼ ∼ ∼ ∼ ∼ ∼ ∼ RML ∼ ∼ ∼ MV × PL × NLC ∼ ∼ UC ∼ ∼ DD × × ∼ ∼ ∼ ∼ × ×
Figure 8: This table illustrates the use of different priors and constraints (whose acronyms are described in Sec. 2.2) as they were introduced chronologically. × represents strict constraints which cannot be violated, while ∼ is a prior or soft constraint. The methods of Tappen et al. [START_REF] Tappen | Recovering intrinsic images from a single image[END_REF][START_REF] Tappen | Estimating intrinsic component images using non-linear regression[END_REF] use a monochromatic illumination constraint, but comparisons shown in two papers exhibit colored illumination [START_REF] Bousseau | User assisted intrinsic images[END_REF][START_REF] Garces | Intrinsic images by clustering[END_REF]. We expect the use of many priors to improve decompositions but on more limited datasets. We compute the LMSE accuracy of various algorithms on our realistic synthetic dataset. These methods are sorted by decreasing average LMSE. The "Villa" And "San Miguel" scenes have consistently lower accuracy (i.e., higher LMSE), while the "Head" and "Babylonian City" scenes have higher accuracy. However, LMSE does not reflect usefulness for computer graphics applications. Figure 14: Decomposition and application results for various methods. Top to bottom: reflectance, illumination, image-edited result. We add a carpet to the first image, and remove the t-shirt's logo on the second and third images. Additional results can be seen in supplemental materials. The baseline consists of the best of three naive approaches.

Figure 2 :

 2 Figure 2: The intrinsic decomposition ambiguity. Top row. An input image I is decomposed into the product of reflectance and illumination layers ρ and L using [BPD09]. Bottom row. The product of these altered reflectance and illumination layers ρ and L exactly reconstructs the input image I.

Figure 3 :

 3 Figure 3: If the monochromatic illumination assumption is used, colored shadows cannot be removed from the reflectance layer. First row. The input image I is decomposed into reflectance and (colored) illumination layers ρ and L. Second row. We desaturate the illumination layer L , and recover the corresponding reflectance image ρ = IL , hence simulating the best intrinsic decomposition achievable assuming grayscale illumination. This new reflectance partly contains the colored shadow. The monochromatic illumination assumption is still widely used to date (table8).

Figure 4 :

 4 Figure4: This illustrates the gradient norm of the reflectance and illumination components. The last column shows the log-histogram of gradient occurrences. A smoother illumination layer would exhibit a steeper decreasing curve for illumination gradients than for reflectance gradients. This is however hardly the case in practice: in cluttered scenes, occlusion boundaries dominate and produce strong gradients in both the reflectance and illumination layers.

  Using compression-based complexity measures, Nicola et al. show reflectance has lower complexity than illumination [SF15]. Similarly, Barron and Malik minimize the entropy of the log-reflectance [BM15] to obtain parsimonious reflectance values.

Figure 6 :

 6 Figure6: A Planckian lighting prior assumes that under skylight, pixels of the same reflectance belong to a single straight line in log-RGB space. We use outdoor linear images and plot log-RGB color distributions of two constant reflectance surfaces (second column). This hardly results in distinct lines. This prior may produce artifacts when used on images taken in uncontrolled settings.

Figure 7 :

 7 Figure 7: Scribbles help disambiguate intrinsic decomposition. Using a single input photograph, we do not expect fully automatic methods to differentiate a picture of an object displayed on a wall from an actual 3D object. From left to right, input photograph, illumination layer obtained using the method of Bonneel et al. [BST * 14] without, and then with user interactions to remove illumination variations within each picture on the wall.

  Temporal consistency. Working on videos, reflectance can be assumed to remain temporally consistent across video frames [LZT * 12, BST * 14]. RGB-D videos allow for easier tracking via the 3d reconstructed scene [LZT * 12], when RGB videos would require an illumination agnostic optical flow in the case of moving shadows [BST * 14] (again, a chicken-and-egg problem). The ambient occlusion model. Using multiple images of the same object under various lighting conditions, Hauagge et al. [HWBS16] replace the illumination term in the intrinsic decomposition by a scaled ambient occlusion term. Ambient occlusion estimates the fraction of the hemisphere visible from any point in the scene, regardless of the lighting conditions.

  [BBS14] with convolutional neural networks: Narihira et al. [NMY15] learn a lightness classifier, Zhou et al. [ZKE15] integrate the learned priors into Bell et al. [BBS14]'s CRF framework, and Zoran et al. [ZIKF15] solve a quadratic program on super-pixels with these data-driven priors. Other preliminary work on deep architectures for intrinsic decompositions are under investigation [Vit15, SL16, LVVG16].

Figure 9 :

 9 Figure 9: Left. The decomposition of Bousseau et al. [BPD09] may help understand the color variations of this popular dress meme [LSHC15]. Right. Although being one of the best performing methods, the priors introduced by Barron et al. [BM15] swap the illumination and reflectance layers. In this example, the priors for the reflectance layer perhaps better describe the illumination layer and vice-versa.

Figure 10 :

 10 Figure 10: Sample images from existing intrinsic decomposition ground-truth datasets. First row. The database of Bell et al. [BBS14] provides reflectance comparison judgments between pairs of points in photographs of real-world settings (mostly indoor photographs, as illustrated by the overlayed graph. Second row. The MIT database [GJAF09] contains simple isolated objects. Their illumination is obtained by painting the object in white, and the reflectance is obtained by taking the image / illumination ratio. Third row. The MIT database has been extended to more complex scenes by Beigpour et al. [BKK15], exhibiting two objects under various illumination conditions. Fourth row. The Sintel dataset [BWSB12], originally for optical flow evaluation, provides reflectance ground truth videos from synthetic renderings. However, taking the ratio image hardly recovers a valid illumination layer. Fifth row. The database of Beigpour et al. [BSV * 13] contains 9 complex but textureless synthetic scenes.

Figure 11 :

 11 Figure 11: Three of our realistic ground-truth intrinsic decomposition results obtained with LuxRender [Ver07].

[

  FDL04] [TFA05] [TAF06] [STL08] [BPD09] [GJAF09] [JSW10] [SY11] [SYJL11] [ZTD * 12] [GMLMG12] [BBS14] [CCFI14] [ZKE15] [BM15] [BHY15] [TNY15] [ZIKF15] MI

Figure 12 :

 12 Figure 12: We compute the LMSE accuracy of various algorithms on our realistic synthetic dataset. These methods are sorted by decreasing average LMSE. The "Villa" And "San Miguel" scenes have consistently lower accuracy (i.e., higher LMSE), while the "Head" and "Babylonian City" scenes have higher accuracy. However, LMSE does not reflect usefulness for computer graphics applications.

Figure 13 :

 13 Figure 13: Computation time for tested automatic intrinsic decomposition methods with respect to image resolution. No trend line is shown for Gehler et al. [GRK* 11] since images were resized to roughly the same resolution. For clarity, we merged results of gray and color Retinex[START_REF] Grosse R | Ground-truth dataset and baseline evaluations for intrinsic image algorithms[END_REF] ; however they show a bimodal timing distribution: this is rather due to L1 and L2 reconstructions.

•

  Reflectance map extraction [HS79, RRF * 15]: This generalizes the intrinsic decomposition problem to the recovery of arbitrary reflectance functions. Rematas et al. offers a deep learning approach that recovers a full hemispherical reflectance function per pixel [RRF * 15] -the reflectance map. • Shadow extraction. This closely related problem consists of detecting and extracting shadows. Under low-frequency lighting conditions, shadows become softer and a precise definition of shadows becomes an issue. Relating this problem to that of intrinsic decomposition, Isaza et al. [ISR12] evaluates intrinsic decomposition methods to detect shadows. • Light estimation. This tries to uncover the lighting conditions of a scene, for example, by extracting directional light sources [LMGH * 13] or environment maps [LE10], from an image. This problem becomes difficult in the presence of localized light sources or inter-reflections, in the absence of a 3d geometric model of the scene.

• Specularity removal. This is a different intrinsic decomposition approach that separates the diffuse from the specular reflection components. While this also extracts intrinsic images in the sense of

Barrow and Tenenbaum [BT78]

, the term intrinsic decomposition now most often refers to the separation of the diffuse from the illumination components (though exceptions exist [BvdW11]). The interested reader may find further information in the survey of Artusi et al.

[START_REF] Artusi | A Survey of Specularity Removal Methods[END_REF]

.

• Color constancy. When a colored object is illuminated by light sources of different colors, for example in a sunset, or in indoor lighting, the object appears to humans as having kept its original color. A common photographic operation is to try and compensate for the light source chromaticity -a process often called "white balancing" or color constancy correction. This is typically achieved globally using a grey card or a colored chart

[START_REF] Cheng | Beyond white: Ground truth colors for color constancy correction[END_REF]

.

In the presence of multiple colored light sources at the same time, this operation can be performed locally [HMP * 08, BBPD12]. A perfectly local color constancy would recover the illuminant color at each pixel, which would correspond to the chromaticity of the illumination layer L(x, λ). • Texture-Structure decomposition. This separates the high frequency textural elements from lower frequency structures

[START_REF] Gilboa | Structuretexture image decomposition-modeling, algorithms, and parameter selection[END_REF]

. The definition of a texture, however, depends on the scale of the observed element. For instance, a forest canopy can belong to the structure if seen from sufficiently close, and it becomes part of textures if seen from sufficiently far. This decomposition has seen applications for intrinsic decompositions [JCTL14, BHY15], and image-based material editing [BBPA15].

The monochromatic illumina- tion assumption is still widely used to date (table 8).

	Retinex (R). Following a series of experiments involving colored
	light sources illuminating patterns resembling Piet Mondrian's paint-
	ings, Land et al. [LM71] determined that the human eye made sense
	of relative intensities to form a representation of material reflectance.
	In this representation, sharp transitions (or edges) are perceived
	as changes in reflectance properties, while smoother variations are
	seen as changes in illumination. They coined the term Retinex, from
	Retina and Cortex, for the eye representation of lightness -a percep-
	tual quantity correlated with reflectance. They further built a Retinex
	Machine reproducing the reflectance of a grayscale color wheel
	based on a 1D strip of sensors. This approach has been extended
	in 2-d by Horn [Hor74] who provides numerical tools and physical
	interpretations of the 2-d Retinex problem. The Retinex model fails
	for hard shadows or occlusion boundaries since they produce strong
	edges associated with illumination variations, but it works reason-
	ably well for smooth surfaces. This prior is often implemented by

c 2017 The Author(s) Computer Graphics Forum c 2017 The Eurographics Association and John Wiley & Sons Ltd.

Acknowledgments

We thank the authors of the intrinsic decomposition methods for having shared their implementation with us and Sean Bell for sharing his evaluation framework. We also thank Julie Digne for initiating the idea of this report, and Adobe for software donations. We would like to thank our funding agencies NSF IIS 1617861, and a Google Faculty Research Award. We acknowledge the use and adaptation of LuxRender scenes from Andrew Price (Kitchen scene), Peter Sandbacka (Hotel Lobby) and Simon Wendsche (School Corridor), PBRT scenes from Jay Hardy (White Room), Guillermo M. Leal Llaguno (San Miguel), Florent Boyer (Villa), Marko Dabrovic and Mihovil Odak (Sponza), and BlendSwap user Wig42 (Modern living room). Some of these PBRT resources were compiled by Benedikt Bitterli, and are availabe in supplemental material. Mitsuba scenes from Johnathan Good (Arabic, Babylonian and Italian Cities).

and Mitsuba [START_REF]Mitsuba renderer[END_REF]. A 3d RGB scatterplot of reflectance values hardly exhibits clusters, while quantizing reflectance values into up to 15 color clusters still shows some artifacts on complex scenes (no dithering was applied -see insets). The reflectance remains mostly flat for man-made scenes, but fails on the head model. The illumination was computed as the ratio between the input and reflectance images, and may reflect inaccuracies on glossy or refractive objects, or due to subsurface scattering.

Figure 15: Decomposition and application results for various methods. Top to bottom: reflectance, illumination, image-edited result. We add a carpet to the first image, and remove the t-shirt's logo on the second and third images. Additional results can be seen in supplemental materials.


Chiara Foglietta

Andrea Gasparri

Stefano Panzieri

A Networked Evidence Theory Framework for Critical Infrastructure Modeling

Keywords: Modeling, evidence theory, situational awareness, data fusion

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

A nation's critical infrastructure comprises complex, interdependent systems whose proper operation and interaction are essential to the welfare of society. Modeling the interdependencies between the different critical infrastructure sectors is a complex task, but this is vital to correctly analyze and predict cascading phenomena.

Interdependencies, from the viewpoint of data fusion, can be analyzed to discern critical events and failures. Events associated with critical infrastructures typically have low probabilities but high impact. The impact is exacerbated by effects that are not localized and tend to propagate to interconnected infrastructures. This implies that a "signature" of the event can be identified in the interconnected system and, sometimes, in independent agencies such as meteorological services or police departments. Therefore, there is a need to fuse together the available data and derive a common belief of the current situation. Developing a common belief facilitates efficient and effective decisions and plans of action. Since critical infrastructures are inherently distributed [START_REF] De Porcellinis | Modeling critical infrastructure via a mixed holistic reductionistic approach[END_REF], the fusion mechanism should also be distributed.

This paper provides an extension of the transferable belief model to facilitate the application of evidence theory. The existing transferable belief model provides a centralized technique for fusing data. However, a decentralized approach is required for interdependent critical infrastructures. This paper demonstrates the application of an extended decentralized transferable belief model to the domain of critical infrastructure protection and shows how data fusion can help develop more accurate beliefs about interdependent infrastructures.

Background

Evidence theory is a methodology that is commonly applied to data fusion problems. Data fusion seeks to combine data from heterogeneous sources or sensors to provide estimates of ongoing events [START_REF] Das | High-Level Data Fusion[END_REF].

Evidence theory stems primarily from the pioneering work of Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] and Shafer [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] and is often considered in the same light as Bayesian networks. The primary difference, however, is the ability to deal with uncertainty [START_REF] Ng | Uncertainty management in expert systems[END_REF]. The Dempster-Shafer theory [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] explicitly considers uncertainty and examines if the various sources of data are inconsistent or if there is an error in the modeling process. On the other hand, Bayesian networks use the recognition of input values as a likelihood from pre-determined patterns. The application of the two approaches depends on the type of knowledge that has to be represented and fused [START_REF] Jensen | Bayesian Networks and Decision Graphs[END_REF].

Several methods have been proposed for combining and correlating the available data in the context of the Dempster-Shafer framework. This work uses the methodology proposed by Smets [START_REF] Smets | The transferable belief network[END_REF], which extends the Dempster-Shafer framework by assuming that the correct answer might not be among the considered ones (i.e., it engages the open world assumption). Smets' approach also allows the computation of the amount of contradictory data in the value of the empty set.

The major limitation to applying evidence theory in a real context is the number of hypotheses required to model the application of interest. This can be explained by the fact that, from a computational perspective, the power set of the set of hypotheses has to be computed, causing the complexity to grow exponentially with the number of hypotheses. However, evidence theory has been successfully applied to a variety of practical problems using approximations (see, e.g., [START_REF] Hall | Generation, combination and extension of random set approximations to coherent lower and upper probabilities[END_REF][START_REF] Voorbraak | A computationally efficient approximation of Dempster-Shafer Theory[END_REF]).

Data fusion can also aid in impact assessment. In fact, limiting an impact assessment strictly to measured events can lead to heavy underestimation or in-correct estimation. For example, an isolated failure can propagate in a number of ways depending on the cause and the detection methods. Examples include a fire blast detected by a sensor and a computer virus detected by an intrusion detection system. In the first example, the fire blast propagation effects are associated with an interdependency model according to a spatial proximity pattern. In the second example, a computer virus may propagate to similar, and highly dispersed, telecommunication nodes.

Evidence Theory

Evidence theory provides a means to form a consolidated belief by correlating evidence from different sources [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A Mathematical Theory of Evidence[END_REF][START_REF] Smets | The transferable belief network[END_REF]. This section provides an overview of the evidence theory formalisms used in this work.

Let ω i represent a cause of system failure and Ω = {ω 1 , • • • , ω n } be the set of hypotheses containing known possible failures. This set is called the "frame of discernment." For example, the possible causes of failures of a critical infrastructure asset include sabotage, device failure, fault due to weather and a (cyber) denial-of-service attack. Note that the hypotheses are assumed to be mutually exclusive in evidence theory.

The power set of the frame of discernment is expressed as The transferable belief model [START_REF] Smets | The transferable belief network[END_REF] is derived from the basic belief mass function m:

Γ(Ω) = {γ 1 , • • • , γ 2 |Ω| },
m : Γ(Ω) → [0, 1].
This function, also called the "basic belief assignment" (BBA), maps each element of the power set to a value between 0 and 1. Each BBA is an atomic element in the transferable belief model. In fact, each sensor, agent or node must be able to assign the BBA values by some subjective assumptions or through algorithms that automatically determine the assignment. The BBA function is constrained by:

γa⊆Γ(Ω) m(γ a ) = 1 with m(∅) = 0.
The transferable belief model examines propositions accordingly as: "the true value of ω i is in γ a " where γ a ∈ Γ(Ω). For γ a ∈ Γ(Ω), m(γ a ) is the confidence that supports exactly γ a . This implies that the true value is in the set γ a ; however, due to the lack of additional data, it is not possible to support any strict subset of γ a .

In the case of different independent data sources, a rule is necessary to aggregate the data. Several rules of combination exist in the literature; the most widely used rules are Dempster's rule [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] and Smets' rule [START_REF] Smets | The transferable belief network[END_REF].

Dempster's rule of combination [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF], which was the first to be proposed, is a purely conjunctive operation. The rule strongly emphasizes the agreement between multiple sources and ignores conflicting evidence using a normalization factor as shown in the following equation:

Dempster{m i , m j }(∅) = 0 Dempster{m i , m j }(γ a ) = γ b ∩γc=γa m i (γ b )m j (γ c ) 1 - γ b ∩γc=∅ m i (γ b )m j (γ c ) ∀γ a ∈ Γ(Ω).
On the other hand, Smets' rule of combination [START_REF] Smets | The transferable belief network[END_REF] provides the ability to explicitly express the contradiction in the transferable belief model by letting m(∅) ̸ = 0. This combination rule, unlike Dempster's rule, avoids normalization while preserving commutativity and associativity:

Smets{m i , m j }(γ a ) = m i (γ a ) ⊗ m j (γ a ) ∀γ a ∈ Γ(Ω)
where

m i (γ a ) ⊗ m j (γ a ) = γ b ∩γc=γa m i (γ b )m j (γ c ) ∀γ a ∈ Γ(Ω).
The relation m(∅) > 0 can be explained in two ways: (i) open world assumption; and (ii) quantified conflict. The open world assumption, proposed by Dempster, reflects the idea that the frame of discernment must contain the true value. If the open world assumption is true, then the set of hypotheses must contain all possibilities. Under this interpretation, if ∅ is the complement of Ω, the mass m(∅) > 0 represents the case where the truth is not contained in Ω. Alternatively, the notion of quantified conflict means that there is some underlying conflict between the sources that are combined to produce the BBA. Hence, the mass assigned to m(∅) represents the degree of conflict. Specifically, it is computed as:

m i (∅) ⊗ m j (∅) = 1 - γa∈Γ,γa̸ =∅ (m i (γ a ) ⊗ m j (γ a )) .

Data Fusion in the Network Context

Consider a network of multiple agents described by an indirect graph G = {V, E} where

V = {v i | i = 1, • • • , n V } is the set of nodes and E = {e ij | (v i , v j )}
is the set of edges that represent a communication channel between the nodes. Edges are indirect and, thus, the existence of an arc e ij implies the existence of an edge e ji .

In this work, we assume that no central unit is available to perform data aggregation. Additionally, communications between nodes are limited to the neighbors of the node under consideration (i.e., nodes that are physically or directly connected to the node under consideration). These assumptions are reasonable for data fusion problems in the area of sensor networks. A direct consequence of the assumptions, however, is that Smets' rule of composition cannot be directly applied. Indeed, applying Smets' rule multiple times over the same BBAs leads to different outputs. Note that this could easily happen in a distributed context where communications are local and limited to the one-hop neighborhood.

Consider, for example, a scenario where it is necessary to identify the degradation of services in a telecommunications network. In the case of extensive delays in packet transmission, it may be necessary to determine if the situation is a temporary congestion or a persistent malicious attack. Network sensors (e.g., intrusion detection systems) can detect a denial-of-service (DoS) attack. The DoS attack may result in cascading effects within the network that are identified by other sensors. If the attack is conducted on a sufficiently large scale, entire regions can be compromised, with different sensors providing different inputs based on localized knowledge.

A general method to define a BBA allocation is to consider the reliability of the data source. Suppose that the data obtained from the source supports a set of hypotheses in Ω. Then, the subset γ a of the power set Γ(Ω), containing the set of hypotheses, receives a mass m(γ a ) equal to the reliability of the source. The remaining mass 1m(γ a ) is assigned to the universal set because no other data is available.

Table 1 shows the BBA assignments for a cause classification problem in a telecommunications network using the centralized approach. The network has three sources (nodes) that relay data to determine the probable causes of a network congestion incident. The table shows the BBA values assigned to the various network nodes. The frame of discernment is Ω = {a, b}, where hypothesis a is a denial-of-service attack (i.e., congestion of one or more network nodes that intentionally degrades the telecommunications network), and hypothesis b indicates congestion in the telecommunications network due to routing problems.

If we assume that Node 1 is coordinating with Node 2, then upon applying Smets' operator, the result is:

m 1 ⊗ m 2 = {0.44, 0.11, 0.44, 0.01}. Algorithm 1 : Gossip Algorithm Data: t = 0, s i (t = 0) ∀i = 1, • • • , n Results: s i (t end ) ∀i = 1, • • • , n
while end condition do Select an edge e ij ∈ E(t) according to e; Update the states of the selected nodes by applying the operator R:

s i (t + 1) = s i (t) ⊗ s j (t) s j (t + 1) = s j (t) ⊗ s i (t) Set t = t + 1 end
If Node 1 then communicates with Node 3 by exchanging data about the possible cause of the congestion, then the result obtained by centralized aggregation is equal to: m 12 ⊗ m 3 = {0.77, 0.095, 0.134, 0.001}.

Next, we consider the communications between Node 1 and Node 3 by applying Smets' operator. The result, which is different from the one obtained in the centralized system, is given by: m 123 ⊗ m 3 = {0.8828, 0.0767, 0.0404, 0.0001}.

If a decision on the cause of the fault has to be made, then the latter case results in a change of opinion; according to the latest aggregations, the cause is a denial-of-service attack (hypothesis a) instead of network routing congestion (hypothesis b).

As shown, Smets' rule of combination cannot be directly used in a distributed data aggregation context. The next section presents a distributed algorithm that can update the knowledge of all the nodes in a network.

Data Fusion Algorithm

The algorithm proposed by Gasparri, et al. [START_REF] Gasparri | A networked transferable belief model approach for distributed data aggregation[END_REF] provides the ability to divide the knowledge of each node into two parts: (i) data shared between two nodes; and (ii) localized data retained by each node.

In the decentralized algorithm, the network is described by an indirect graph G = {V, E}. A spanning tree T = {V, E} is derived, where E ⊆ E is available to all nodes. Note that the nodes must have the capacity to save data. Interested readers are referred to [START_REF]A distributed algorithm for constructing minimal spanning trees[END_REF] for details about spanning tree construction.

Communications between nodes are asynchronous and follow the Gossip Protocol [START_REF] Boyd | Randomized gossip algorithms[END_REF]. This protocol is formalized as Algorithm 1. The triplet {S, R, e} specifies the network such that: S is the set of local states of each node in the network.

R is local interaction rule, i.e., for each pair of nodes (i, j) such that e ij ∈ E, the following equation holds: R : R q × R q → R q where q is the number of elements called "focal sets." e is the process of edge selection for which e ij ∈ E(t) is the selected edge at time t.

To define the operator R, it is first necessary to introduce the operator ⊙. Consider two sets of BBAs:

m k = {m k (γ a ) | ∀γ a ∈ Γ(Ω)} and m i = {m i (γ a ) | ∀γ a ∈ Γ(Ω)}, such that m k = m i ⊗ m j .
The operator ⊙ can then be defined as:

m j = m k ⊙ m i mi k .
Starting with the power set element with the highest cardinality and examining elements with decreasing cardinality, the value of a BBA can be computed recursively as follows:

m j (γ a ) = m k (γ a ) - γ b ∩γc=γa,γa⊂γ b m j (γ b )m i (γ c ) γa⊆γ b m i (γ b )
.

It is now possible to introduce the operator R, along with ⊕, to aggregate BBAs of the nodes:

m i (t + 1) = m j (t + 1) = m i (t) ⊕ m j (t) = { mj i (t, γ a ) ⊗ mj i (t, γ a ) ⊗ mi,j (t, γ a ), ∀γ a ∈ Γ(Ω)}.
Note that the term mj i (t, γ a ) indicates the innovation of node i with respect to the node j, which can be calculated recursively using the operator S:

mj i (t, γ a ) = m i (t, γ a ) ⊙ mi,j (t, γ a ).
The element mi,j (t, γ a ) express the common knowledge (i.e., knowledge exchanged between two agents (i, j) after the last aggregation) such that:

mi,j (t, γ a ) = n = {0, 0, • • • , 0, 1}.
Gasparri, et al. [START_REF] Gasparri | A networked transferable belief model approach for distributed data aggregation[END_REF] have shown that the algorithm converges to the same result as the centralized aggregation algorithm. The convergence time of the distributed algorithm is related to the diameter d of the spanning tree. Additionally, the computational complexity of the R-operator is the same as that of Smets' operator.

Application Scenario

Applying the transferable belief model involves modeling a problem, specifying the frame of discernment and selecting the BBAs. Note that the assignment of BBAs is problem dependent and depends significantly on the source of information and knowledge about the system. In this work, we assume that experts provide advice on assigning the BBAs. Other possible approaches are described in [START_REF] Das | High-Level Data Fusion[END_REF][START_REF] Veremme | Belief assignment on compound hypotheses within the framework of the transferable belief model[END_REF], where the reliability of the data sources is considered along with the effect of mass assignment on the compound hypotheses.

Our scenario involves a supervisory control and data acquisition (SCADA) system that aggregates data from a large number sensors positioned at remote locations. Operators located at a control center manage operations by acquiring system data and updating system parameters. Note that the remote facilities all have data regarding critical events and that the data is generated from different sources.

Communications between the control center and the devices in the field occur via dedicated telecommunications circuits or shared public media such as the Internet. Note that the same communications channels can support information sharing among critical infrastructures to help discern the possible causes of failures. The concept of information sharing across sectors can help prevent cascading effects or prevent the impact from propagating to interconnected infrastructures that have not yet been affected.

Our scenario considers n = 5 interdependent critical infrastructures. Each infrastructure owner determines the BBAs in his/her infrastructure. To discern the cause of a fault, the BBAs have to be aggregated and shared among the five infrastructures. The frame of discernment is Ω = {a, b, c}, where hypothesis a represents a cyber attack, hypothesis b represents the failure of an isolated single unit, and hypothesis c represents a natural disaster (e.g., earthquake). Table 2 shows the BBA assignments for the five infrastructure nodes.

First, we evaluate the scenario using the centralized algorithm. The results as shown in Table 3. Each column presents the results obtained using Smets' operator, with column "NN 12" representing the aggregation of Nodes 1 and 2, column "NN 123" representing the aggregation of Nodes 1, 2 and 3, and so on. The last column "C-TBM" shows the results for the centralized transferable belief model.

Next, we evaluate the scenario using the distributed algorithm. Table 4 shows the edge selection data. Table 5 shows the output of the R-operator. Each column corresponds to the aggregation of two nodes as related to the sequential timing. After each exchange of knowledge between two nodes, the R-operator reveals that the two nodes have the same knowledge value as the operator output. For example, at t = 5, the nodes have the same value: {0.6334, 0.0451, 0.3070, 0.0125, 0.0014, 0.0, 0.0004, 0.0001} and are consistent with the centralized transferable belief model outputs. Note that the algorithm terminates when the nature of the edge selection process is known [START_REF] Gasparri | A networked transferable belief model approach for distributed data aggregation[END_REF]. The results demonstrate that the decentralized approach (Table 5) yields the same values as the centralized approach (Table 3).

Conclusions

The decentralized extension of the transferable belief model enables the application of evidence theory to data fusion in critical infrastructure applications. This is important because centralized computation is ill-suited to critical infrastructure applications due to the associated interdependencies. The case study shows that the decentralized approach produces the same results as the centralized approach, and also demonstrates the utility of fusing data in a distributed manner for interdependent critical infrastructures.

Integrating situational awareness with distributed monitoring is an appealing concept for networked critical infrastructures. This is important because the ability to leverage and share sensor data can significantly enhance system resilience and robustness. For example, in the case of a cyber attack, data from intrusion detection systems coupled with data from standard field sensors can be combined to obtain more accurate belief assessments.

The decentralized approach offers the same advantages as the traditional transferable belief model in terms of its ability to deal with uncertainty. It is important to note, however, that the same disadvantages exist (e.g., exponential growth in the computational complexity with respect to the number of hypotheses). Nevertheless, both approaches are useful for modeling beliefs in interdependent infrastructures for the purpose of enhancing situational awareness.

  which has cardinality |Γ(Ω)| = 2 |Ω| . The power set contains all possible subsets of Ω, including the empty set γ 1 = ∅ and the universal set γ 2 |Ω| = Ω.

Table 1 .

 1 BBA assignments for a telecommunications network.

	Set	Node 1 Node 2 Node 3 m12 m123
	∅	0.0	0.0	0.0	0.44 0.770
	{a}	0.1	0.5	0.7	0.11 0.095
	{b}	0.8	0.4	0.2	0.44 0.134
	{a,b}	0.1	0.1	0.1	0.01 0.010

Table 3 .

 3 Centralized algorithm output with incremental aggregation.

	Set	NN 12 NN 123 NN 1234 NN 12345 C-TBM
	∅	0.36	0.468	0.5304	0.6334	0.6334
	{a}	0.18	0.108	0.0648	0.0451	0.0451
	{b}	0.34	0.346	0.3676	0.3070	0.3070
	{c}	0.04	0.046	0.0308	0.0125	0.0125
	{a,b}	0.06	0.024	0.0048	0.0014	0.0014
	{a,c}	0.0	0.0	0.0	0.0	0.0
	{b,c}	0.0	0.006	0.0012	0.0004	0.0004
	{a,b,c}	0.20	0.002	0.0004	0.0001	0.0001

Table 4 .

 4 Temporal edge selection.

	Time t=1 t=2 t=3 t=4 t=5 t=6 t=7
	Edge	e12	e23	e34	e45	e34	e23	e12

Table 5 .

 5 Distributed algorithm output.

	Set	s1 ⊕ s2 s2 ⊕ s3 s3 ⊕ s4 s4 ⊕ s5 s3 ⊕ s4 s2 ⊕ s3 s1 ⊕ s2
	∅	0.36	0.468	0.5304	0.6334	0.6334	0.6334	0.6334
	{a}	0.18	0.108	0.0648	0.0451	0.0451	0.0451	0.0451
	{b}	0.34	0.346	0.3676	0.3070	0.3070	0.3070	0.3070
	{c}	0.04	0.046	0.0308	0.0125	0.0125	0.0125	0.0125
	{a,b}	0.06	0.024	0.0048	0.0014	0.0014	0.0014	0.0014
	{a,c}	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	{b,c}	0.0	0.006	0.0012	0.0004	0.0004	0.0004	0.0004
	{a,b,c}	0.20	0.002	0.0004	0.0001	0.0001	0.0001	0.0001

Acknowledgement

This research was partially supported by the European Commission through the FP7 Cockpit CI Project.