Mariem Graa
email: mariem.graa@it-sudparis.eu

Nora Cuppens-Boulahia

Frédéric Cuppens
email: nora.cuppens@telecom-bretagne.eu

Ana Cavalli
email: ana.cavalli@it-sudparis.eu

Detection of illegal control flow in Android System: Protecting private data used by Smartphone Apps

Keywords: Android system, Smartphones, Dynamic analysis, Static analysis, Control dependencies, Leakage of sensitive information

. But this technique does not detect control flows. In particular, untrusted applications can circumvent Android system and get privacy sensitive information through control flows.

In this paper, we propose a hybrid approach that combines static and dynamic analysis to propagate taint along control dependencies in Android system. To evaluate the effectiveness of our approach, we analyse 27 free Android applications. We found that 14 of these applications use control flows to transfer sensitive data. We successfully detect that 8 of them leaked private information. Our approach creates 19% performance overhead that is due to the propagation of taint in the control flow. By using our approach, it becomes possible to detect leakage of personal data through control flows.

Introduction

Smartphone operating systems use has been increasing at an accelerated rate in recent years. Android surpassed 80% market share in the third quarter of 2013 [START_REF] Van Der Meulen | Gartner says smartphone sales accounted for 55 percent of overall mobile phone sales in third quarter of 2013[END_REF] and it is the most targeted OS by the cyber criminals with more than 98% of malware applications [START_REF] Enck | Taintdroid: An information-flow tracking system for realtime privacy monitoring on smartphones[END_REF]. This is due to the prevalence of third party app stores (48 billion apps have been installed from the Google Play store in May 2013 [START_REF] News | Bbc google activations and downloads update may 2013[END_REF]). These applications are used to capture, store, manipulate, and access to data of a sensitive nature in mobile phone. An attacker can launch control flow attacks to compromise confidentiality of the Android system and can leak private information without user authorization. In the study presented in the Black Hat conference, Daswani [START_REF] Wilson | Many android apps leaking private information[END_REF] analyzed the live behavior of 10,000 Android applications and showed that more than 800 of them were found to be leaking personal data to an unauthorized server. Therefore, there is a need to provide adequate security mechanisms to control the manipulation of private data by third-party apps. Many mechanisms are used to protect sensitive data in the Android system, such as the dynamic taint analysis that is implemented in TaintDroid [START_REF] Enck | Taintdroid: An information-flow tracking system for realtime privacy monitoring on smartphones[END_REF]. The principle of dynamic taint analysis is to "taint" some of the data in a system and then propagate the taint to data for tracking the information flow in the program. Two types of flows are defined: explicit flows such as x = y, where we observe an explicit transfer of a value from x to y, and implicit flows (control flows) shown in Figure 1 were there is no direct transfer of value from a to b, but when the code is executed, b would obtain the value of a. The dynamic taint analysis mechanism does not detect control flows which can cause an under tainting problem i.e. that some values should be marked as tainted, but are not. The under tainting problem can cause a failure to detect a leak of sensitive information. Thus, malicious applications can bypass the Android system and get privacy sensitive information through control flows. In a previous work [START_REF] Graa | Detecting control flow in smarphones: Combining static and dynamic analyses[END_REF], we have proposed an approach that combines static and dynamic taint analysis to propagate taint along control dependencies and to track implicit flows in the Google Android operating system. Our approach enhances the TaintDroid approach by tracking control flows in the Android system to solve the under-tainting problem. In this paper, we present implementation details and experimental results of the proposed approach. We show effectiveness of our approach to propagate taint in the conditional structures of real Android applications and to detect leakage of sensitive information. This paper is organized as follows: Section 2 presents a motivating example. We discuss related work about static and dynamic taint analysis and we analyze existing solutions to solve the under tainting problem in Section 3. Section 4 describes the proposed approach and the corresponding implementation details. We analyse a number of Android applications to test the effectiveness of our approach and we study our approach taint tracking overhead in Section 5. Finally, Section 6 concludes with an outline of future work.

An attacker can exploit an indirect control dependencies to leak private data. Let us consider the control dependence attack shown in Figure 2. The variable X contains the private data that is the user contact. The attacker tries to get the user contact name by comparing it with symboles of Ascii table in the second loop. He stored the character of private data founded in Y . At the end of the loop, the variable Y contains the correct value of the user contact and it is not tainted because taint is not propagated in the control flow statement. The attacker exploits untainted variable that should be tainted (under tainting problem) to leak private data. Thus, Y is leaked through the network without being detected. Therefore, an attacker can leak a sensitive information by exploiting control flows.

Related Work

Many security mechanisms are used to protect sensitive data in smartphones. TaintDroid [START_REF] Enck | Taintdroid: An information-flow tracking system for realtime privacy monitoring on smartphones[END_REF], an extension of the Android mobile-phone used to control in realtime the manipulation of users personal data by third-party applications. It implements a dynamic taint tracking and analysis system to track the information flow and to detect when sensitive data leaves the system. AppFence [START_REF] Hornyack | These aren't the droids you're looking for: retrofitting android to protect data from imperious applications[END_REF] extends Taintdroid to implement enforcement policies. One limit of TaintDroid and AppFence approaches is that they cannot propagate taint in control dependencies. The methods proposed in [START_REF] Egele | Detecting privacy leaks in ios applications[END_REF][START_REF] Chin | Analyzing inter-application communication in android[END_REF][START_REF] Fuchs | Scandroid: Automated security certification of android applications[END_REF] statically analyze third party application code for detecting data leaks. But, these static analyses approaches cannot capture all runtime configuration. Some approaches combine static and dynamic analysis to solve the under-tainting problem. BitBlaze [START_REF] Song | Bitblaze: A new approach to computer security via binary analysis[END_REF] presents a novel fusion of static and dynamic taint analysis techniques to track all information flow. DTA++ [START_REF] Kang | Dta++: Dynamic taint analysis with targeted control-flow propagation[END_REF] uses Bitblaze and enhances the dynamic taint analysis to limit the under-tainting problem. However DTA++ is evaluated only on benign applications. Trishul [START_REF] Nair | A virtual machine based information flow control system for policy enforcement[END_REF] correctly identifies control flow to detect a leak of sensitive information. Furthermore, these approaches are not implemented in smartphones application. Fenton [START_REF] Fenton | Memoryless subsystem[END_REF] defined a Data Mark Machine, an abstract model, to handle control flows. This model does not take into account the implicit flow when the branch is not executed. Denning [START_REF] Denning | Secure information flow in computer systems[END_REF] enhances the run time mechanism used by Fenton with a compile time mechanism to solve the under-tainting problem. Denning inserts updating instructions whether the branch is taken or not. We draw our inspiration from the Denning approach, but we define formally a set of taint propagation rules to solve the under-tainting problem. In [START_REF] Graa | Detecting control flow in smarphones: Combining static and dynamic analyses[END_REF], we propose a hybrid approach that tracks control flows in smartphones. We define a set of formal propagation rules to solve the under-tainting problem. We prove the correctness and completeness of these rules and we propose a correct and complete algorithm to solve the under tainting problem [START_REF] Graa | Formal Characterization of Illegal Control Flow in Android System[END_REF]. In [START_REF] Graa | Protection against Code Obfuscation Attacks based on control dependencies in Android Systems[END_REF], we show that our approach can resist to code obfuscation attacks based on control dependencies in the Android system using the taint propagation rules. But, we do not evaluate the overhead and effectiveness of our approach and do not test a real Android applications. We provide test exprimental results and we evaluate the overhead and the false positives of our approach. We show that it successfully detects sensitive information leakage by untrusted Android applications.

Approach Overview and Implementation

Our objective is to detect private information leakage by untrusted smartphone applications exploiting implicit flows. We control the manipulation of private data by third party application in realtime.

Our approach consists of two main components: the StaticAnalysis component and the DynamicAnalysis component (see Figure 3). We implement our proposed approach in the TaintDroid operating system. We add a StaticAnalysis component in the Dalvik virtual machine verifier that statically analyzes instructions of third party application Dex code at load time. Also, we modify the Dalvik virtual machine interpreter to integrate the DynamicAnalysis component. We implement the two additional rules using native methods that define the taint propagation.

Static Analysis Component

In this component, we check the instructions of methods to create the control flow graph (CFG). A CFG is composed of basic blocks and edges. The basic blocks represent nodes of the graph. The directed edges represent jumps in the control flow. For each control instruction, we insert a BasicBlock at the end of the basic blocks list. Then, we specify the target of basic blocks. We perform the post dominator analysis (A node v is post-dominated by a node w in the control flow graph G if every path from v to Exit, not including v, contains w) on the control flow graph to determine the flow of the condition dependencies from different blocks. After this, we allocate a BitmapBits for tracking condition dependency. We store the control flow graph using the DOT language of graphviz tool [START_REF] At | Graphviz[END_REF] in the data directory of the smartphone.

Dynamic Analysis Component

The dynamic analysis is performed at run time by instrumenting the Dalvik virtual machine interpreter. We assign a context taint to each basic block. The context taint includes the taint of the condition on which the block depends. We compare arguments in the condition using the following instruction : res cmp = ((s4)GET REGIST ER(vsrc1) cmp (s4)GET REGIST ER(vsrc2)). Based on the comparison result, we verify wether the branch is taken or not. We Combine the taints of different variables of the condition as follows: SET REGIST ER T AIN T (vdst, (GET REGIST ER T AIN T (vsrc1)|GET REGIST ER T AIN T (vsrc2))) to obtain the Context T aint. If res cmp is not null then the branch is not taken. Thus, we adjust the ordinal counter to point to the first instruction of the branch by using the function ADJU ST P C(2). Otherwise, it is the second branch (else) which is not taken then we adjust the ordinal counter to point to the instruction in this branch by using the function ADJU ST P C(br) where br represents the branch pointer. We instrument different instructions in the interpreter to handle conditional statements. For each instruction, we taint the variable to which we associate a value (destination register). In the case of f or and while loops, we process by the same way but we test whether the condition is still true or not in each iteration. We make a special treatment for Switch instructions. We deal with all case statements and all instructions which are defined inside Switch instructions. Note that, we only taint variables and do not modify their values. Once we handle all not taken branches, we restore the ordinal counter to treat the taken branches and we assign taints to modified variables in this branch. We make a special exception handling to avoid leaking information. If the type of exception that occurred is listed in a catch block, the exception is passed to the catch block. So, an edge is added in the CFG from the throw statement to the catch block to indicate that the throw statement will transfer control to the appropriate catch block. If an exception occurs, the current context taint and the exception's taint are stored. The variables assigned in any of the catch blocks will be tainted depending on the exception's taint.

Evaluation

In this section, we analyse a number of Android applications to test the effectiveness of our approach. Then, we study our taint tracking approach overhead using standard benchmarks. We evaluate the false positives that could occur using our approach. We use a Nexus One mobile device running Android OS version 2.3 enhanced to track implicit flows.

Effectiveness

To evaluate the effectiveness of our approach, we analyse 27 free Android applications downloaded from the Android Market [1] that manipulated private data. As shown in Table 1, five applications require permissions for contacts and five applications require permissions for camera at install time.

Most of these applications access to locations and phones identity. Also, our analysis showed that these permissions are acquired by the implicit or explicit consent of the user. For example, in the weather application, when the user selects the option "use my location", she gives permission to the application to use and to send this information to the weather server. We found that 14 of these 25 analyzed Android applications (marked with * in the Table 1) leak private information:

-The IMEI numbers that identify a specific cell phone on a network is one of the information that is transmitted by 11 applications. Nine of them do not present an End User License Agreement (EULA). -Two applications transmitted the device's phone number, the IMSI and the ICC-ID number to their server. -The location information is leaked by 15 third-party applications to advertisement servers. These applications do not require implicit or explicit user consent. Just two applications require an EULA. We use dex2jar tool [2] to translate dex files of different applications to jar files. Then, we use jd-gui [START_REF]Java decompiler[END_REF] to obtain the source code that will be analysed. As shown in Table 2, we found that 14 of tested Android applications listed by types of accessed sensitive data use control flows to transfer private information. Eight of them leaked private data. Sensitive data is used in the if ,f or and while control flow instructions. We verify that variables to which a value is assigned in these instructions and that depend on a condition containing private data are not tainted using TaintDroid. Our approach has succesfully propagated taint in these control instructions and detected leakage of tainted sensitive data that is reported in the alert messages.

Performance

In this part of the paper, we study our taint tracking approach overhead. The static analysis is performed at load and verification time. At load time, our approach adds 33% overhead with respect to the unmodified system. At verification and optimization time, our approach adds 27% overhead with respect to the unmodified system. This time increase is due to the verification of method instructions and the construction of the control flow graphs in the static analysis phase. We install the Caffeine-Mark application [START_REF] Corporation | Caffeinemark 3.0[END_REF] in our Nexus One mobile device to determine the java microbenchmark. Note that the CaffeineMark scores roughly correlate with the number of Java instructions executed per second and do not depend significantly on the amount of memory in the system or on the speed of a computers disk drives or internet connection [START_REF] Corporation | Caffeinemark 3.0[END_REF]. 4 presents the execution time results of a Java microbenchmark. We propagate taint in the conditional branches especially in the loop branches and we add instructions in the processor to solve the under tainting problem. Then, the loop benchmark in our approach presents the greatest overhead. We taint results of arithmetic operations in explicit and control flows. Thus, the arithmetic operations present the greatest overhead. The string benchmark difference between unmodified Android system and our approach is due to the additional memory required in the string objects taint propagation. We observe that the unmodified Android system had an overall score of 3625 Java instructions executed per second. Whereas, our approach had an overall score of 2937 Java instructions executed per second. Therefore, our approach has a 19% overhead with respect to the unmodified system.

False positives

Our analysis and tests indicated that almost of 50% of studied Android applications use control flows and leak sensitive data. Our approach generates 25% of false positives. We detect an IMSI leakage vulnerability when it is really used as a configuration parameter in the phone. Also, we detect that the IMEI is transmitted outside of smartphone but it is the hash of the leaked IMEI. Thus, we can not treat these applications as privacy violations.

Conclusion

In order to detect the leakage of sensitive information by third-party apps exploiting control flows in smartphones, we have proposed a hybrid approach that propagates taint along control dependencies to solve under-tainting problem. We have analysed 27 free Android applications to evaluate the effectiveness of our approach. We found that 14 applications use control flows to transfer sensitive data and 8 leak private information. We showed that our approach generates significant false positives that can be reduced by considering expert rules (ad hoc rules). Also, we can use an access control approach to authorize or not the transmission of the data outside the system. Our approach incurs 19% performance overhead that is due to the propagation of taint in the control flow. To improve performance of our system, we suggest implementing the taint propagation mechanism in Just In Time Compiler (JIT) that provides better performance than the interpreter such as a minimal additional memory usage. By implementing our approach in Android systems, we successfully protect sensitive information and detect most types of software exploits caused by control flows.

Fig. 1 .

 1 Fig. 1. Implicit flow example.

Fig. 2 .

 2 Fig. 2. Control dependence attack

Fig. 3 .

 3 Fig. 3. Our Approach Architecture

Fig. 4 .

 4 Fig. 4. Microbenchmark of java overhead

Figure

  Figure 4 presents the execution time results of a Java microbenchmark. We propagate taint in the conditional branches especially in the

Table 1 .

 1 Third party applications grouped by the requested permissions (L: location, Ca: camera, Co: contacts, P: phone state).

	Third party applications	Permissions
		L Ca Co P
	The Weather Channel * ; Cestos; Solitaire; Babble; Manga	x
	Browser (5)
	Bump; Traffic Jam; Find It * ; Hearts; Blackjack; Alchemy;	x	x
	Horoscope * ; Bubble Burst Free; Wisdom Quotes Lite * ; Paper
	Toss * ; Classic Simon Free; Astrid * (12)
	Layar ; ProBasketBall (5) x x	x
	Wertago * ; Dastelefonbuch * ; RingTones * ; Yellow Pages * ;	x	x x
	Contact Analyser (5)

* ; Knocking * ; Coupons * ; Trapster *

Table 2 .

 2 Third party applications used control flows

	Category	application	Leaked
		Name		Data
		Wertago		x
	Contact and Phone Iden-	Dastelefonbuch x
	tity
		Yellow Pages	x
	Camera	Knocking	x
		ProBasketBall
		The	Weather	x
		Channel
		Cestos
		Classic	Simon
		Free
		Bubble	Burst
		Free
	Location and Phone Iden-	Bump
	tity
		Traffic Jam
		Horoscope	x
		Paper Toss	x
		Find It		x


Xingbo Wang
email: xingbo.wang@hig.no.

Marius Pedersen

Jean-Baptiste Thomas

THE INFLUENCE OF CHROMATIC ABERRATION ON DEMOSAICKING

Keywords: demosaicking, chromatic aberration, colour filter array

The wide deployment of colour imaging devices owes much to the use of colour filter array (CFA). A CFA produces a mosaic image, and normally a subsequent CFA demosaicking algorithm interpolates the mosaic image and estimates the full-resolution colour image. Among various types of optical aberrations from which a mosaic image may suffer, chromatic aberration (CA) influences the spatial and spectral correlation through the artefacts such as blur and mis-registration, which demosaicking also relies on. In this paper we propose a simulation framework aimed at an investigation of the influence of CA on demosaicking. Results show that CA benefits demosaicking to some extent, however CA lowers the quality of resulting images by any means.

INTRODUCTION

The advent of colour filter array (CFA) by Bayer [START_REF] Bayer | Color imaging array[END_REF] permits a single-sensor system capable of capturing a colour image at one exposure. The advantages mentioned above are gained at the expense of reduced spatial image resolution due to the intrinsic property of CFAs that merely one colour channel is sensed by each pixel. Therefore certain post processing becomes necessary to estimate other channels to a certain extent, known as demosaicking.

Optical images are commonly distorted by optical elements in form of various types of optical aberrations, and one of them resulting in colour artefacts is denoted chromatic aberration (CA) [2, p. 257]. It occurs because lenses, typically made of glass or plastic, bear different refractive indices for different wavelengths of light (the dispersion of the lens). The refractive index decreases with increasing wavelength. The main consequence of CA in imaging is that light rays at different wavelengths are focused at different image distances (axial/longitudinal CA) and at different locations in the image (transverse/lateral CA). In theory and reality only rays at a certain wavelength are focused accurately on the image plane, all other rays are focused before or behind the image plane, which leads to blur when observed from the image plane. Similarly, magnification of a lens varies from wavelength to wavelength as well. On the optical axis, there appears as axial CA, since magnification is zero. In a plane perpendicular to the optical axis, the difference of magnification rate turns more prominent as the distance from the optical axis increases, which produces transverse CA. The problem of CA may be corrected optically to some extent by using superior glass substrate whose refraction index varies to a lower degree and combining positive and negative lenses, nevertheless it is not avoidable in practice where cost, compactness, and weight matter.

CA influences the spatial and spectral correlations of images through the artefacts, such as blur and mis-registration, which also affect demosaicking. To the best of our knowledge, an evaluation of the influence of CA on demosaicking has not been performed. Therefore we address this issue in this work by building a modular framework that simulates the CFA imaging process.

The rest of the paper begins with a description of the experimental framework that includes an introduction to the simulation of CA in Section 2.1, followed by an explanation of the demosaicking algorithms under investigation in Section 2.2. Section 2.3 analyses the image quality metrics employed in the experiment, and the experimental setup is detailed in Section 2.4. Results are demonstrated and discussed in Section 3. Conclusions are drawn and future work is proposed in Section 4.

EXPERIMENTAL FRAMEWORK

The framework consists of four modules: CA simulation, CFA mosaicking, demosaicking and quality evaluation. Three of them are detailed in this section, while the conditions of CFA mosaicking are referred to in Section 2.4.

Simulation of CA

CA occurs in the formation of an optical image, therefore the best way of observing and studying CA is through various optical design. However, this is not flexible and convenient. In consequence, we opt to simulate CA with spectral images as an alternative solution. The two types of aberrations are simulated individually and jointly. The simulation of axial and transverse CAs are demonstrated in Fig. 1 and Fig. 3. Keeping in mind the reason for axial CA, we suppose the rays at 550 nm are focused on the sensor, such that green channel at 550 nm remains in focus and thus little changed, which is consistent with the peak in the CIE luminosity function. Then the rays at other wavelengths will be focused at a range of planes before and behind the intended image plane. This is seen as a blur. CA is an optical aberration closely bound to the optical design of a lens, and the degree of blur depends primarily on the position of the image plane for a given lens at given aperture size. It is often analysed in optical design and fabrication by means of ray tracing. However it varies from one optical design to another, and there is not a parametric model of such aberration to the best of our knowledge. To that end, we rely on the ISETBIO toolbox [START_REF]ISETBIO toolbox[END_REF] that makes use of an ocular CA model [START_REF] Thibos | Effect of ocular chromatic aberration on monocular visual performance[END_REF][START_REF] Thibos | The chromatic eye: a new reduced-eye model of ocular chromatic aberration in humans[END_REF]. With this reduced eye model, the chromatic refractive error, that is the dioptric difference of refraction, can be derived by Eq. 1,

Defocus in diopters

∆R E (λ) = n(λ 1 ) -n(λ 2 ) rn D (1)
where λ 1 and λ 2 refer to the wavelengths in question, and n(λ 1 ) and n(λ 2 ), the refractive indices at these two wavelengths. n D is the refractive index for the sodium D-line (589 nm), which is wavelength for which the model eye is emmetropic. r is the corneal radius of curvature, set to 5.5 mm. The refractive index is obtained by Eq. 2,

n(λ) = a + b λ -c (2)
where a = 1.320535, b = 0.004685, c = 0.214102. Now the defocus obtained is expressed in diopters. To make it suitable for the next step, which is the calculation of the pupil function, it is necessary to convert the defocus expressed in diopters to that expressed in micrometres by Eq. 3.

D m = r 2 D d 16 √ 3 (3)
where r refers to the corneal radius of curvature in millimetres, and D m and D d denote defocus in micrometres and diopters respectively. The pupil function, or sometimes referred to as aperture function, is a complex function of the position in the pupil or aperture (an iris in this case) that represents the amplitude and phase of the wavefront across the pupil, also known as wavefront aberrations. It is an important tool to study optical imaging systems and their performance. Optical system aberrations have historically been described, characterized, and catalogued by power series expansions, where the wave aberration is expressed as a weighted sum of power series terms that are functions of the pupil coordinates. Each term is associated with a particular aberration or mode, such as spherical aberration, coma, astigmatism, field curvature, distortion, and other higher order modes. Zernike polynomials form a complete set of functions or modes that are orthogonal over a circle of unit radius and are convenient for serving as a set of basis functions. Pupil functions are calculated here with the use of Zernike polynomials [START_REF] Maeda | Zernike polynomials and their use in describing the wavefront aberrations of the human eye[END_REF] and pupil function at a given wavelength is demonstrated in Eq. 4.

P F (λ) = Ae -i2π Z λ (4)
where A represents the amplitude, calculated entirely based on the assumed properties of the Stiles-Crawford effect, and the exponential function is actually the phase of the aberration. The Zernike polynomials, denoted by Z, may consist of a series of modes, however, as the only aberration we are interested in this context is defocus (4 th mode), we simply neglect others.

The pupil function is related to the point spread function (PSF) by the Fourier transform [7, p. 131]. The reduced eye model simplifies the aberration and assumes the axial CA is shift-invariant, the generated PSF is static for each wavelength. An example of PSFs at a few wavelengths over the visible spectrum is demonstrated in Fig. 2. It is clear that the support of the PSFs vary in size as the wavelength increases, and achieves the minimum at 550 nm. A PSF is in fact the impulse response of an imaging system in spatial domain, and the Fourier transform of a PSF leads to the optical transfer function that integrates the modulation transfer function (MTF) and the phase transfer function (PhTF). A convolution of a PSF and an image will produce the simulated image distorted by axial CA. The mechanism of transverse CA lies in a varying rate of geometric distortion, that is in fact a radial operation depending on wavelengths. Therefore we first represent the images of each band in polar coordinates. Subsequently the radius can be manipulated in such a manner that the magnification varies in radial direction, following the ray tracing analysis of a real lens design. In practice, the measurement of transverse CA is often provided in form of a series of samples by the ray tracing software, which reflects the wavelength dependent radial distortion, as shown in Table 1. Accordingly, interpolation becomes essential to map and resample the original band images in order to simulate transverse CA.

As mentioned above, the reduced eye model simulates merely axial CA, and transverse CA depends on the results of ray tracing analysis. When combined, the two types of CA may be simulated sequentially as implemented by ISET [START_REF] Farrell | Digital camera simulation[END_REF]. Namely, transverse CA is applied prior to axial CA.

Demosaicking algorithms under investigation

Decades of research on CFA based imaging methods has resulted in a large number of evolving demosaicking algorithms. They can be well classified into two groups according to the intrinsic properties of the images they make use of. To be specific, some algorithms take advantage of spatial corre- lation merely within one channel, whereas others benefit from spectral correlation between colour channels. Here we implemented and experimented with two representative methods, namely bilinear interpolation and gradient-corrected linear interpolation. The former relies on merely spatial correlation, whereas the latter also makes use of, and thus sensitive to, spectral correlation. Indicated by its name, bilinear interpolation is an extension of linear interpolation for interpolating functions of two variables (e.g., x and y) on a regular 2D grid. In fact, bilinear interpolation estimates the unknown values by means of a distance weighted average of its neighbouring pixels. It is applied to each spectral band individually, as a result, it should be sensitive to variation of spatial correlation, however insensitive to alteration of spectral correlation.

Gradient based demosaicking solutions are aimed at reduced artefacts by avoiding interpolating across the edges. Malvar et al. [START_REF] Malvar | High-quality linear interpolation for demosaicing of bayer-patterned color images[END_REF] advance a gradient-corrected bilinear interpolated approach, with a gain parameter to control how much correction is applied. In other words, the results of bilinear interpolation is corrected by a measure of the gradient for the known colour at the pixel location. To determine appropriate values for the gain parameters, a Wiener approach is employed which computes the values leading to minimum meansquare errors, given second order statistics computed from the Kodak data set [START_REF]Kodak lossless true color image suite[END_REF].

Evaluation of demosaicking performance

Evaluation of the demosaicking algorithms can be seen in general as an image quality evaluation task.

In the last decade when demosaicking gains much attention, the evaluation of demosaicking also attracts the interest of the academic community. Longère et al. conduct perceptual quality evaluation on a few demosaicking approaches, and state that perceptual results cannot be easily predicted using an image metric [START_REF] Longère | Perceptual assessment of demosaicing algorithm performance[END_REF]. Even so, there are several attempts for objective quality evaluation. Lu and Tan propose two types of quality measures specifically for demosaicking, one computes the PSNR (peak signal-to-noise ratio) and CIE ∆E * ab for edge and smooth regions respectively, whereas the other one deals in particular with the zipper artefact [START_REF] Lu | Color filter array demosaicking: new method and performance measures[END_REF]. Yang et al. show that low-level features (colours or edges) extracted from demosaicked images are affected by resolution, whatever the demosaicing method used, and therefore propose new criteria designed for low-level image analysis [START_REF] Yang | Quality evaluation of color demosaicing according to image resolution[END_REF].

The most widely used criteria for the evaluation of demosaicing quality are MSE (mean squared error) and PSNR, primarily because they are simple thus easy to implement. However, the PSNR criterion provides a general estimation of the demosaicking quality, but does not really reflect the human judgement. Alternatives such as CIE ∆E * ab and S-CIELAB [START_REF] Zhang | A spatial extension of CIELAB for digital color image reproduction[END_REF] are also widely used, however they require apriori information like the reference white or illuminant [START_REF] Losson | Comparison of color demosaicing methods[END_REF].

In this work, we concern mostly the signal fidelity of each colour channel. Therefore we base our judgement on PSNR. In comparison with PSNR, the Structural SIMilarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] index provides more structural information and tends to be more consistent with subjective image quality assessment, which matches CA that blurs and shifts structures. Therefore we also adopted the SSIM as a method for measuring the similarity between the images before and after demosaicking.

In contrast to straightforward assessment of demosaicking methods or image quality, an evaluation of the influence of CA on demosaicking seems a bit more complicated, as there are two types of source images, i.e., with and without CAs. In consequence, with the two full reference quality metrics, we not only compute the image quality between the source and demosaicked images with and without distortions, but also cross compare the non-distorted source images and distorted demosaicked images.

Experimental setup

The experiment was conducted in the following manner. First, simulated CA was applied to spectral images separately. Then, a mosaic image is formed by filtering the distorted spectral image with a certain CFA mosaic and converting it to sRGB colour space. Next, the mosaic image was interpolated with one of the demosaicking methods. Finally, the demosaicked images were fed into a quality metric channel by channel.

Conditions and parameters in relation to the experiment are presented as follows. and the images were cropped into square matrices by the shorter dimension. Spectral range was set to 400 nm to 700 nm, with an interval of 10 nm. The CFA conforms to the common design of Bayer. Transverse CA was simulated according to real lens data provided by the ISET [START_REF] Farrell | Digital camera simulation[END_REF]. The SSIM index was computed with the empirical Fig. 4. Thumbnails of the images used. Scene 1-5 were selected from the Foster set [START_REF] Foster | Frequency of metamerism in natural scenes[END_REF] and scene 6-11 from the CAVE set [START_REF] Yasuma | Generalized assorted pixel camera: Postcapture control of resolution, dynamic range and spectrum[END_REF].

formula specified on the official webpage to determine the scale for images viewed from a typical distance [START_REF]The ssim index for image quality assessment[END_REF].

RESULTS AND DISCUSSION

Experimental results are depicted below. Lines and markers in red, green and blue represent the results for the corresponding channels, and those in black indicate the average values. Due to the space limitation, some similar results are not shown here. Concerning the performance with and without CAs, Fig. 56show the results in terms of PSNR, while Fig. 7-8 display the results in terms of SSIM index. Similarly, Fig. 9-10 demonstrate the results in terms of PSNR and SSIM respectively for the cross comparison between non-distorted sources and distorted demosaicked images. It can be seen from Fig. 5-6 that combined and axial CA yield higher performance, regardless of the demosaicking Fig. 6. PSNR computed from straight comparisons between distorted source images and distorted demosaicked images for R/G/B channels and the average with gradient corrected linear interpolation. methods used. This is true for all red, green and blue channels. However, in Fig. 5 the best result was achieved by the blue channel under axial and combined aberrations, whereas in Fig. 6 for the green channel under the same aberrations. Such findings are consistent with Fig. 78. This indicates that CA surprisingly benefit demosaicking. In particular, blur benefit both demosaicking methods by higher spatial correlation. In theory, transverse CA in form of mis-registration should reduce correlation, thus decreasing the performance. On the contrary, the figures show opposite results. Perhaps the resampling involved blurs the images to some extent as well. In comparison with intra-channel methods, interchannel interpolation is more sensitive to the content of the green channel. The drop of Image 10 in Fig. 7-8 may be because of the specific characters of the image, i.e., oil painting as the foreground on top of black background. Furthermore, the two distinct models of CA simulation may also affect the results. Suggested by Fig. 9 and Fig. 10, axial CA results in higher performance in terms of cross comparison, indicating that blur influences less the signal fidelity and similarity. In comparison, both of the other two types of CA involve transversal CA and lead to significantly worse results. Obviously CA lowers the overall image quality by any means.

CONCLUSION

In this paper we propose a framework to investigate the influence of chromatic aberration (CA) on demosaicking. Experimental results show that CA benefits demosaicking to some extent, however any type of CA decreases image quality by means of blur and mis-registration.

Certainly the results deserve further study. More demosaicking algorithms and sample images would be helpful. And experiments in subjective assessment of image quality are expected to be conducted, which may reveal better how human observers perceive the issue. Further, the degree of CA may be varied to obtain more data, and ray tracing analysis of real lenses is expected to make the simulation more realistic.

Fig. 1 .

 1 Fig. 1. Workflow for the simulation of axial CA.

Fig. 2 .Fig. 3 .

 23 Fig. 2. A series of PSFs at visible wavelengths.

Fig. 5 .

 5 Fig. 5. PSNR computed from straight comparisons between distorted source images and distorted demosaicked images for R/G/B channels and the average with bilinear interpolation.

Fig. 8 .

 8 Fig. 8. SSIM indices computed from straight comparisons between distorted source images and distorted demosaicked images for R/G/B channels and the average with gradient corrected linear interpolation.

Fig. 9 .

 9 Fig.9. PSNR computed from cross comparisons between non-distorted source images and distorted demosaicked images for R/G/B channels and the average with gradient corrected linear interpolation.

Fig. 10 .

 10 Fig. 10. PSNR computed from cross comparisons between non-distorted source image and distorted demosaicked image for R/G/B channels and the average with gradient corrected linear interpolation.

Table 1 .

 1 Positions of radial samples before and after transverse CA.

	Radial samples	400 nm	Distorted radial samples 475 nm 550 nm 625 nm	700 nm
	0.0518	0.0507	0.0510	0.0511	0.0512	0.0513
	0.1037	0.1014	0.1019	0.1022	0.1024	0.1025
	0.1555	0.1519	0.1527	0.1531	0.1534	0.1537
	0.2073	0.2023	0.2033	0.2039	0.2044	0.2047
	0.2591	0.2525	0.2538	0.2545	0.2551	0.2554
	0.3110	0.3025	0.3040	0.3049	0.3055	0.3060
	0.3628	0.3521	0.3539	0.3549	0.3557	0.3562
	0.4146	0.4014	0.4034	0.4046	0.4054	0.4061
	0.4664	0.4503	0.4525	0.4539	0.4548	0.4555
	0.5183	0.4987	0.5012	0.5028	0.5038	0.5046
	0.5701	0.5467	0.5495	0.5511	0.5523	0.5531
	0.6219	0.5942	0.5972	0.5990	0.6003	0.6012
	0.6737	0.6411	0.6443	0.6463	0.6477	0.6487
	0.7256	0.6874	0.6909	0.6930	0.6945	0.6956
	0.7774	0.7330	0.7368	0.7391	0.7407	0.7418
	0.8292	0.7780	0.7821	0.7846	0.7862	0.7875
	0.8810	0.8224	0.8267	0.8293	0.8311	0.8324
	0.9329	0.8660	0.8705	0.8734	0.8753	0.8767
	0.9847	0.9089	0.9137	0.9167	0.9187	0.9202
	1.0365	0.9511	0.9561	0.9593	0.9614	0.9629


Xingbo Wang
email: xingbo.wang@hig.no

Philip J Green

Jean-Baptiste Thomas

Jon Y Hardeberg

Pierre Gouton

Keywords: colorimetric performance, colour filter array, multispectral imaging, single-sensor

Single-sensor colour imaging systems mostly employ a colour filter array (CFA). This enables the acquisition of a colour image by a single sensor at one exposure at the cost of reduced spatial resolution. The idea of CFA fit itself well with multispectral purposes by incorporating more than three types of filters into the array which results in multispectral filter array (MSFA). In comparison with a CFA, an MSFA trades spatial resolution for spectral resolution. A simulation was performed to evaluate the colorimetric performance of such CFA/MSFA imaging systems and investigate the trade-off between spatial resolution and spectral resolution by comparing CFA and MSFA systems utilising various filter characteristics and demosaicking methods including intraand inter-channel bilinear interpolation as well as discrete wavelet transformed based techniques. In general, 4-band and 8-band MSFAs provide better or comparable performance than the CFA setup in terms of CIEDE2000 and S-CIELAB colour difference. This indicates that MSFA would be favourable for colorimetric purposes.

Evaluation of the colorimetric performance of single-sensor image acquisition systems employing colour and multispectral filter array 1 Introduction

Single-sensor trichromatic imaging systems mostly employ a colour filter array (CFA) in order to sense a portion of the incoming spectra selectively on a pixel-by-pixel basis. An example of CFA that has achieved commercially notable success is known as Bayer filter mosaic consisting of three types of filters, i.e., red, blue and green [START_REF] Bayer | Color imaging array[END_REF], as shown in Figure 1. Thanks to the spatial and spectral inter-pixel correlation an image may possess, the lost information about the incident stimuli can be estimated through demosaicking. Demosaicking is an operation in the image processing chain carried out on the mosaicked image read from the sensor. Consequently each pixel will comprise three components, thereby recovering a full colour image. The success gained by the CFA based single-sensor colour imaging systems has awakened particular interest from the academia and the industries in generalise this concept to the multispectral domain by integrating more than three types of filters into one filter array, which results in the multispectral filter array (MSFA). Two instances of MSFAs can be seen in Figure 2 and3. In general, the development of a MSFA based imaging system involves design of filter transmittances, spatial arrangement of mosaic patterns paired with an associated demosaicking algorithm, and a regression process to recover colorimetric or spectral information.

In recent years, sustained research effort went into designing multispectral filter array (MSFA) [START_REF] Miao | A generic method for generating multispectral filter arrays[END_REF] and developing associated demosaicking algorithms [START_REF] Miao | A generic binary tree-based progressive demosaicking method for multispectral filter array[END_REF][START_REF] Baone | Demosaicking methods for multispectral cameras using mosaic focal plane array technology[END_REF][START_REF] Brauers | A color filter array based multispectral camera[END_REF][START_REF] Monno | Multispectral demosaicking using adaptive kernel upsampling[END_REF][START_REF] Monno | Multispectral demosaicking using guided filter[END_REF][START_REF] Wang | Median filtering in multispectral filter array demosaicking[END_REF][START_REF] Wang | Discrete wavelet transform based multispectral filter array demosaicking[END_REF]. Also widely explored is the influence of filter characteristics on colour/spectrum reproduction [START_REF] Parulski | Color filters and processing alternatives for one-chip cameras[END_REF][START_REF] Imai | Comparison of spectrally narrow-band capture versus wide-band with a priori sample analysis for spectral reflectance estimation[END_REF][START_REF] Wang | A study on the impact of spectral characteristics of filters on multispectral image acquisition[END_REF][START_REF] Wang | Multispectral imaging: narrow or wide band filters[END_REF]. Nevertheless, to our best knowledge, little is known about how filter design and demosaicking algorithms affect the colorimetric performance of a CFA and particularly MSFA image acquisition system.

In comparison with CFAs, MSFAs populate higher number of channels, thus reducing the number of pixels assigned to a certain channel for a given sensor. Obviously this lowers spatial resolution, however MSFA may offer higher spectral resolution. While the former effect generally lowers the colorimetric performance of the system, the latter should improve the accuracy of colour reproduction. It is therefore of particular interest to evaluate the colorimetric performance of such MSFA imaging systems and investigate the trade-off between spatial resolution and spectral resolution by comparing CFA and MSFA systems utilising various filter characteristics and demosaicking methods.

The following sections of the paper are organised as follows. We first present the methods used including the simulated framework illustrated in Section 2.1, the filter design strategy described in Section 2.2, the mosaic generation demonstrated in Section 2.3, the demosaicking algorithms introduced in Section 2.4, the colour and spectral reflectance estimation method for device calibration presented in Section 2.5, the evaluation means explained in Section 2.6 and the experimental conditions listed in Section 2.7. Results are shown in Section 3 that leads to reasonable conclusions drawn in Section 4.

Methods

Simulated framework

A simulated workflow was constructed so as to conduct the research due to practical difficulties in physical implementation of the MSFAs [START_REF] Lapray | Multispectral filter arrays: Recent advances and practical implementation[END_REF]. As shown in Figure 4, the framework consists of a chain of processing that starts from the hyperspectral images used as virtual optical images. A MSFA mounted sensor is merely a combination of the filter array and the image sensor, and the mosaicked optical image is therefore formed in between. As a result, a hyperspectral image can be considered as a spectrally sampled optical image which will then be spectrally filtered and spatially interleaved by the filter array. In this manner, the process of mosaicking is simulated. Next, an ideal image sensor populating the same number of pixels as the images integrates the incident power at each pixel over the spectrum. Quantum efficiency of the sensor is integrated with the spectral transmittance of filters, so that filter characteristics referred to in this work represent sensor sensitivities. And neither optical crosstalk nor optical/electronic noise is considered in this work.

Experimental procedures

Mosaic image, namely the sensor output, is actually a digital representation of a spatially or spectrally subsampled and interleaved trichromatic or multispectral image. Therefore, it needs to be interpolated spatially and/or spectrally, in order to recover the information lost in the mosaicking process. In other words, interpolation yields a colour or multispectral image of full spatial resolution.

At this stage, the recovered image is not colorimetrically meaningful as the digital counts have not been assigned any physical meaning. This is addressed by device colour calibration. It first models the device by associating stimuli with known colorimetric characteristics and the corresponding sensor response, and later estimates original colorimetric information of the unknown stimuli from the corresponding sensor response.

CIE tristimulus values are computed before colour difference between the original and reproduced hyperspectral images is calculated.

Filter design in MSFA systems

Among the derivatives of Bayer mosaic, some possess complementary colour filters in comparison to the primary colour filters utilised in the original patent [START_REF] Bayer | Color imaging array[END_REF].

Literature presents distinct results. It is obvious that complementary colour filters intrinsically bear wider pass-band than their primary counterparts, and it is widely accepted that the former gives rise to better colour reproduction and signal-to-noise ratio in sufficient lighting conditions, whereas the latter offers higher sensitivity as well as resolution [START_REF] Parulski | Color filters and processing alternatives for one-chip cameras[END_REF][START_REF]Panasonic: Primary color and complementary color filters[END_REF]. Our results, nevertheless, show that appropriate pass-bands outperforms some narrower ones [START_REF] Wang | A study on the impact of spectral characteristics of filters on multispectral image acquisition[END_REF][START_REF] Wang | Multispectral imaging: narrow or wide band filters[END_REF]. In addition, our previous research on multispectral demosaicking poses the question of filter design in relation to the inter-band correlation [START_REF] Wang | Median filtering in multispectral filter array demosaicking[END_REF][START_REF] Wang | Discrete wavelet transform based multispectral filter array demosaicking[END_REF]. An instance of 4-band filter set used in this research is depicted in Figure 5. Following the aforementioned findings, we are interested in narrowband and broadband bandpass filters as well as corresponding inverted ones such as bandstop filters. The FWHM (Full Width at Half Maximum) of passband and stopband have been set to 10 nm and 40 nm respectively. In practice, a passband of 10 nm simulates very narrow bandpass filters like LCTF (Liquid Crystal Tunable Filter), a stopband of 10 nm mimics notch filters relying on destructive interference. Similarly, a passband and a stopband of 40 nm resemble the spectral transmittances of thin-film filters.

In addition to the filters mentioned above, we introduced two more types. One is based on the principle that the transmittances of filters should sample the spectrum evenly with their FWHM. The other is in fact the result produced by a filter selection algorithm [START_REF] Hardeberg | Acquisition and Reproduction of Color Images: Colorimetric and Multispectral Approaches[END_REF] that chooses a given number of optimal(or suboptimal) filters from a set of available candidates that are physically practicable, on the assumption that high spectral performance is yielded by the "brightest" filter that transmits the most light combined with other filters which are orthogonal to each other in a vector space. Here we employ a set of transmittance data measured from Wratten filters produced by Kodak, as shown in Figure 6.

Mosaic pattern generation

Filter arrays experimented with in this project were designed with the help of a generic binary tree based generation method of MSFA spatial arrangement starting from a chequerboard pattern introduced by Miao et al. [START_REF] Miao | A generic method for generating multispectral filter arrays[END_REF]. By manipulating the pattern through a combination of decomposition and subsampling steps, the method presented may generate MSFAs that satisfy varied design requirements proposed by the authors including probability of appearance, spectral consistency and uniform distribution. It is shown, through case studies, that most of the CFAs currently used by the industry can be derived as special cases.

MSFA demosaicking

A large number of CFA demosaicking algorithms have been proposed over the decades [START_REF] Menon | Color image demosaicking: An overview[END_REF]. In this work, we used intra-channel bilinear interpolation, channel difference bilinear interpolation [START_REF] Adams | Interactions between color plane interpolation and other image processing functions in electronic photography[END_REF] and discrete wavelet transform (DWT) based demosaicking [START_REF] Driesen | Wavelet-based color filter array demosaicking[END_REF].

Bilinear interpolation makes use of merely intra-pixel correlation and estimates unknown colour components by exploiting the spatial correlation between sampled colours in a certain spectral plane, and function plane by plane.

Channel difference interpolation brings inter-channel correlation into play and interpolates the difference between one colour plane and another [START_REF] Adams | Interactions between color plane interpolation and other image processing functions in electronic photography[END_REF] on the assumption that hue changes smoothly in images.

DWT transforms an image into various frequency bands, and natural images often possess rather similar high-frequency information among these bands, which provides yet another solution to demosaicking problem. This algorithm has been extended to a 4-band MSFA as reported in the literature [START_REF] Wang | Discrete wavelet transform based multispectral filter array demosaicking[END_REF].

Device calibration

Colour calibration is an inverse problem aimed at an estimation of the tristimulus values of the stimuli from the corresponding measurements obtained from sensors. In concrete terms, a colour acquisition process can be described in a linear form as

R = QS (1)
where R refers to the sensor responses, Q corresponds to system responsivities and S represents the incoming stimuli.

Colour calibration aims at an estimation of S from R. Equation 1 is solvable if Q is known and invertible, however it is not true in the case of colorimetric calibration. Nevertheless it can be estimated by means of training where a collection of training stimuli S t and corresponding responses R t are utilised to derive an approximation of Q -1 . In this work, we employed the method of linear least squares [START_REF] Hironaga | Estimating the noise variance in an image acquisition system and its influence on the accuracy of recovered spectral reflectances[END_REF] which attempts to solve (1) by means of pseudoinverse which leads to (2)

S = S t R + t R (2)
where S is an estimation of S and

R + t is a right pseudoinverse of R t : R + t = R T t (R t R T t ) -1 .

Performance evaluation

An evaluation of the colorimetric performance of CFA and MSFA based imaging systems can be solved by means of colour difference formula. However, considering that the targets are digital images rather than uniform colour patches, a metric incorporating some low-level HVS features, such as S-CIELAB, might be suitable as well and may provide more information. CIEDE2000 is the latest colour difference formula developed and recommended by the CIE [START_REF]Colorimetry -Part 6: CIEDE2000 Colour-Difference Formula[END_REF]. It further improves perceptual uniformity in comparison with the CIE94 formula by introducing a few revised compensation terms for lightness, chroma and hue respectively. In addition, there are three corresponding parameters that are usually set to 1:1:1 and can be adjusted according to specific applications. For instance, CIE recommended 2:1:1 for textile industry. In this work, we used 2:2:1 to evaluate image colour difference.

The S-CIELAB metric extends the CIELAB Delta E metric to colour images by adding a spatial pre-processing step to the standard CIE ∆E a * b * metric to account for the spatial-colour sensitivity of the human eye [START_REF] Zhang | A spatial extension of cielab for digital color-image reproduction[END_REF]. It measures how accurate the reproduction of a colour is to the original when viewed by a human observer.

Experimental conditions

The experiments were conducted in such conditions as follows. Spectral range covers the spectrum between 400 nm and 700 nm with 10 nm interval. CIE D65 was used as the illuminant. Among the 48 hyperspectral images used as virtual scenes, 16 are from Foster database [START_REF] Foster | Frequency of metamerism in natural scenes[END_REF] and 32 are from CAVE database [START_REF] Yasuma | Generalized assorted pixel camera: Post-capture control of resolution, dynamic range and spectrum[END_REF]. Three types of MSFA were considered, namely 3-band CFA, 4-band and 8-band. For the least-square regression, 170 spectral reflectances of natural objects [START_REF] Vrhel | Measurement and analysis of object reflectance spectra[END_REF] and the corresponding CIE XYZ tristimulus values were utilised as the training targets. Tristimulus values were calculated with colour-matching functions for the CIE 1931 standard colorimetric observer [START_REF]CIE Colorimetry -Part 1: Standard Colorimetric Observers[END_REF]. For the calculation of S-CIELAB colour difference, the viewing distance was set to 60 cm, and the resolution was set to 95.78 dpi, so as to mimic a 23-inch LCD monitor of 1920 × 1080 pixels and an aspect ratio of 16:9.

Results and discussion

Results are presented in Figure 7 and8. It is of great moment to realise that the colour difference shown here reflect the overall performance of the system consisting of filter characteristics, spatial arrangements, demosaicking methods as well as colour estimation techniques. However, a comparative analysis of the results reveal some clues.

From the results we can observe that increased number of bands in general offer lower or comparable colour difference especially when paired with 10 nm and 40 nm bandpass filters and a selected range of Wratten filters. In particular, the 40 nm bandpass filters result in the lowest colour difference among all of the methods and configurations, whereas the 10 nm bandstop filters yield significantly larger errors.

In general, the DWT based demosaicking outperforms the other two where the widths of passband are significantly broader, whereas bilinear interpolation carried out on channel differences does not perform satisfactorily. This is also related to insufficient high-frequency components in the image databases [START_REF] Wang | Discrete wavelet transform based multispectral filter array demosaicking[END_REF], as inter-channel demosaicking should benefit from inter-channel correlation at high frequencies and broader passbands may boost this correlation.

In most cases S-CIELAB results coincide with CIE DE2000 ones, although the former tends to exaggerate the discrepancy of the results between the CFA and the MSFAs.

Conclusion

A simulation was performed to investigate the colorimetric performance of MSFA based image acquisition systems. In total, 48 virtual scenes were captured by the simulated camera for 3-band, 4-band and 8-band MSFAs respectively, each paired with 6 different types of filter characteristics and 3 demosaicking algorithms. Results were all transformed to CIE XYZ tristimulus values and evaluated with CIEDE2000 (2:2:1) and S-CIELAB colour difference. CIEDE2000 and S-CIELAB results coincide in most cases. In general, the 4-band MSFA provides better or comparable performance in comparison with the 3-band setup except the case of 10 nm bandreject and normal bandpass filters. Similarly the 8-band MSFA delivers higher colour accuracy expect the case of 10 nm and 40 nm bandreject filters. Therefore MSFA is generaly helpful for an application where colorimetric reproduction is required.

Moreover, it is obvious to see that spectral characteristics of a filter set not only make a direct impact on the colour reconstruction, but also influence the spectral correlation of the observed image on which some demosaicking methods depend.

Certainly the validity of the results obtained in this work is limited by the realisticness of the simulation. In a real system, the presence of various types of noise will probably impact the results. However, this laid a foundation for the design of a new MSFA sensor planned in 2015 [START_REF] Lapray | Multispectral filter arrays: Recent advances and practical implementation[END_REF].

Fig. 1 .

 1 Fig. 1. Bayer CFA. Fig. 2. 4-band MSFA. Fig. 3. 8-band MSFA.

Fig. 4 .

 4 Fig. 4. The experimental framework.

Fig. 5 .

 5 Fig. 5. Transmittances of a 4-band filter set.

Fig. 6 .

 6 Fig. 6. Transmittances of a set of Wratten filters.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Average CIEDE2000 colour difference among test images.


Pierre-Henri Conze

François Rousseau

Vincent Noblet

Fabrice Heitz

Riccardo Memeo

Patrick Pessaux

Semi-automatic liver tumor segmentation in dynamic contrast-enhanced CT scans using random forests and supervoxels

Pre-operative locoregional treatments (PLT) delay the tumor progression by necrosis for patients with hepato-cellular carcinoma (HCC). Toward an efficient evaluation of PLT response, we address the estimation of liver tumor necrosis (TN) from CT scans. The TN rate could shortly supplant standard criteria (RECIST, mRECIST, EASL or WHO) since it has recently shown higher correlation to survival rates. To overcome the inter-expert variability induced by visual qualitative assessment, we propose a semi-automatic method that requires weak interaction efforts to segment parenchyma, tumoral active and necrotic tissues. By combining SLIC supervoxels and random decision forest, it involves discriminative multi-phase cluster-wise features extracted from registered dynamic contrast-enhanced CT scans. Quantitative assessment on expert groundtruth annotations confirms the benefits of exploiting multi-phase information from semantic regions to accurately segment HCC liver tumors.

Introduction

Hepato-cellular carcinoma (HCC) is the most common type of liver cancer and the third most frequent cause of cancer-related death. Pre-operative locoregional treatments (PLT) tend to downstage HCC tumors by necrosis. Standard evaluation scores (RECIST, mRECIST, EASL or WHO) used to predict the response to PLT do not provide fully satisfactory results [START_REF] Ronot | Alternative response criteria (Choi, EASL and mRECIST) versus RECIST1.1 in patients with advanced hepatocellular carcinoma treated with Sorafenib[END_REF]. A more efficient HCC patient follow-up is reached through tumor necrosis (TN) rate which provides more significant correlation with survival rates.

To overcome inter-expert variability induced by visual qualitative assessment, we present a computed-aided diagnosis method for TN rate computation. Assessing TN remains an open issue due to a wide variability of shape, size, location, contour aspect and intensity of HCC tumors. It requires the segmentation of healthy liver parenchyma as well as tumoral active and necrotic areas (Fig. 4).

For this task, dynamic contrast-enhanced images provide discriminative information since HCC is characterized by arterial enhancement followed by venous washout in response to contrast agent injection [START_REF] Ronot | Hepatocellular carcinoma: Diagnostic criteria by imaging techniques[END_REF].

Dynamic contrast-enhanced images have been recently exploited for HCC liver tumor segmentation via extraction of multi-phase voxel-wise features used in level sets [START_REF] Lee | Estimation of necrosis volumes in focal liver lesions based on multi-phase hepatic CT images[END_REF], k-means [START_REF] Raj | Visualization and segmentation of liver tumors using dynamic contrast MRI[END_REF] or graph cuts [START_REF] Fang | Segmentation of liver tumor using efficient global optimal tree metrics graph cuts[END_REF]. These multi-phase features capture the dynamic in response to agent injection and tend to build a full perfusion model from low temporal resolution data (Fig. 3b). However, to reach a better accuracy, user interaction appears necessary since arterial enhancement and venous washout depend on contrast agent kinetic and injection protocol [START_REF] Ronot | Hepatocellular carcinoma: Diagnostic criteria by imaging techniques[END_REF]. In this direction, tumor extraction has been covered in an interactive perspective for single-phase images [START_REF] Shimizu | Ensemble segmentation using Adaboost with application to liver lesion extraction from a CT volume[END_REF][START_REF] Geremia | Spatial decision forests for ms lesion segmentation in multi-channel mr images[END_REF] by relying on supervised ensemble learning with single-phase voxel-wise features including spatial characteristics. Although more appropriate, such strategy is voxel-wise as in [START_REF] Lee | Estimation of necrosis volumes in focal liver lesions based on multi-phase hepatic CT images[END_REF][START_REF] Raj | Visualization and segmentation of liver tumors using dynamic contrast MRI[END_REF][START_REF] Fang | Segmentation of liver tumor using efficient global optimal tree metrics graph cuts[END_REF] and therefore requires a significant amount of interaction while relying on a limited spatial context

In this work, we propose to exploit robust multi-phase cluster-wise features extracted from registered multi-phase contrast-enhanced CT scans by combining clustering and supervised ensemble learning. Performing interactive learning and prediction on semantic regions allow tumor segmentation to take advantage of discriminative dynamic information at an extended spatial extent. Moreover, it ensures weak interactions efforts for practitioners. This method applied to TN rate estimation is a key step toward accurate PLT response assessment.

Methodology

Traditional approach: single-phase voxel-wise random forest

In an interactive setting for image segmentation, the user usually defines manually a set of K labeled voxels S = {v k , c(v k )} k∈{1,...,K} where c(v k ) ∈ {c 1 , ..., c N } is the label at voxel v k with N the total class number. The set S , referred as training data, is used to build a voxel/label mapping model whose aim is to predict the label c(v ) of each test voxel v . In the context of ensemble learning methods, random decision forests [START_REF] Breiman | Random Forests[END_REF] have grown in popularity due to their ability to offer a unified framework for many machine learning tasks [START_REF] Criminisi | Decision forests: A unified framework for classification, regression, density estimation, manifold learning and semi-supervised learning[END_REF].

A random decision forest consists of T independent trees made of both internal nodes which split input data according to binary tests and terminal nodes which reach all together a final data partition. At each internal node, the split sends voxels to left and right children nodes. For this task, the associated binary test focuses on a randomly subset θ(v ) of the visual features θ(v ) extracted for input voxels v and halve the input dataset according to:

h(v , θ) = true, if τ low < θ(v ) < τ up false, otherwise (1)
where θ(v ) is compared to thresholds τ low and τ up . In the context of tumor segmentation, methods relying on random decision forest such as [START_REF] Geremia | Spatial decision forests for ms lesion segmentation in multi-channel mr images[END_REF][START_REF] Cuingnet | Automatic detection and segmentation of kidneys in 3D CT images using random forests[END_REF] usually focus on one single image and therefore involve single-phase visual features θ(v ). During training, each tree takes the training voxel set S as input and optimizes its own internal nodes ({τ low , τ up , θ(v )}) via information gain maximization [START_REF] Criminisi | Decision forests: A unified framework for classification, regression, density estimation, manifold learning and semi-supervised learning[END_REF] to obtain the most discriminative binary tests with respect to S . After this optimization, each leaf node l t of the t th tree receives a partition S lt of the training data S and produces an entire class probability distribution: p lt (c i |S ) ∀i ∈ {1, ..., N }. To predict the label c(v ) of a given test voxel v with associated single-phase features θ(v ), v is injected into each optimized tree which makes it reaches a leaf node l t per tree following split rules. For each label c i , we get:

p(c(v ) = c i ) = 1 T T t=1 p lt (c(v ) = c i |S ) = 1 T T t=1 |{v k , c(v k )} ∈ S lt | c(v k ) = c i | |S lt | (2) The final predition of c(v ) corresponds to the label c i maximizing p(c(v ) = c i ): c(v ) = arg max ci p(c(v ) = c i ) (3)
2.2 Proposed methodology: multi-phase cluster-wise random forest

From voxels to semantic regions. In the traditional approach (Sec. 2.1), voxels are acting without inter-dependencies which may result in a lack of spatial consistency regarding classification results. Features are sometimes explicitly related to spatial context [START_REF] Geremia | Spatial decision forests for ms lesion segmentation in multi-channel mr images[END_REF][START_REF] Cuingnet | Automatic detection and segmentation of kidneys in 3D CT images using random forests[END_REF] but spatial extent remains limited. A-posteriori regularization techniques such as conditional random field (CRF) [START_REF] Bauer | Fully automatic segmentation of brain tumor images using support vector machine classification in combination with hierarchical conditional random field regularization[END_REF] would more accurately introduce spatial constraints but it strongly increases computational complexity. Moreover, such voxel-wise learning-based approaches require a significant amount of interaction for the end-user to get a large enough training set. To reduce interaction efforts and intrinsically introduce spatial consistency, we perform random forest on training and test semantic regions instead of voxels.

Exploiting multi-phase input data. Since practitioners focus HCC diagnosis on the association of both arterial hypervascularity and venous washout [START_REF] Ronot | Hepatocellular carcinoma: Diagnostic criteria by imaging techniques[END_REF], we propose to take full advantage of multi-phase contrast-enhanced CT scans. In practice, for a given examination, a contrast agent is injected to the patient and CT scans are acquired at different phases : before injection (WI) but also after at arterial (AR), early venous (EV) and late venous (LV) phases. Each examination The set of registered multi-phase images (Fig. 1) is S Ĩ = { ĨWI , ĨAR , I EV , ĨLV } where . denotes warped scans. To obtain S Ĩ from S I , we apply a symmetric diffeomorphic non-rigid registration based on the variational formulation of [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]. A liver segmentation mask (Fig. 1) is assumed to be available at EV phase.

Proposed protocol. Our interactive tumor segmentation method translates in multi-phase cluster-wise random forest classification for TN rate estimation. N = 3 classes are considered: parenchyma, active and necrotic areas. First, the liver volume is over-segmented using a 3D extension of the simple linear iterative clustering (SLIC) superpixel algorithm [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]. Starting from all the voxels belonging to the liver mask in I EV , SLIC provides a set of K R 3D compact clusters R = {r i } i∈{1,...,K R } (Fig. 2a). Then, the interactive training cluster selection occurs. Instead of brushing strokes on many voxels, the practitioner has only to select and label a subset of R (Fig. 2b). We claim that such interaction is more suitable for clinical practice due to its simplicity and its possible integration in multi-examination training. It results in a training cluster set S = {r j , c(r j )} j⊂{1,...,K R } combining cluster r j with groundtruth labels c(r j ) ∈ {c 1 , c 2 , . . . , c N }. A random forest is then built based on the training set S via a training procedure for which internal nodes are optimized with respect to multi-phase cluster-wise visual features θ(r j ) assigned to each cluster r j .

Test clusters r (r ∈ R\S ) are then propagated into the forest to get a label prediction c(r ) (Fig. 2c) based on their own multi-phase cluster-wise features.

Finally, the TN rate τ is computed as follows:

τ = r |r |.1 c(r )=c 1 r |r |.[1 c(v )=c 0 +1 c(v )=c 1 ]
where |r | is the number of voxels in cluster r , c 0 and c 1 active and necrosis labels.

Multi-phase cluster-wise features. The accuracy of our method is related to the ability of the features to discriminate the different tissues. We assign to each semantic region, r , 20 multi-phase cluster-wise visual features divided into 3 groups (Fig. 3a). The two first groups introduce spatial characteristics in terms of intensity and gradient magnitude at the cluster spatial extent. In both cases and for each phase, mean and standard deviation values are computed among all voxels of r . Thus, we quantify intrinsic intensity, visual homogeneity, textural information and texture-scale repartition. The third group fully exploits dynamic contrast-enhanced data by combining several multi-phase intensity information averaged among the voxels of r : dynamic features used in [START_REF] Lee | Estimation of necrosis volumes in focal liver lesions based on multi-phase hepatic CT images[END_REF][START_REF] Raj | Visualization and segmentation of liver tumors using dynamic contrast MRI[END_REF][START_REF] Fang | Segmentation of liver tumor using efficient global optimal tree metrics graph cuts[END_REF] including peak enhancement (PE) and area under curve (AUC) (Fig. 3b) to which is added inter-phase intensity differences ∆ EV /AR and ∆ LV /EV since PE only compares WI and AR. The baseline pre-contrast (BL) (Fig. 3b) is taken into account with intensity features. These dynamic features discriminate clusters based on their own dynamic characterized by arterial enhancement and venous washout.

Results

Evaluation on clinical data. The protocol has been tested on data collected from 7 examinations {e 1 , e 2 , . . . , e 7 } performed on patients with HCC. Each examination results in a set of dynamic contrast-enhanced CT scans including I WI , I AR and I EV with additional I LV for e 1 and e 2 . For each e i , 6 equally reparted 2D axial slices have been selected in I EV to cover the tumor spatial extent and labeled by 4 experts in hepato-digestive surgery to reach groundtruth (GT) segmentation masks delimitating parenchyma, active and necrotic tissues. Manual segmentation is performed on temporally averaged images over all phases. It results in a database of 42 slices with associated fused GT masks obtained by fusing all the expert annotations with STAPLE [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF]. This database has been created since, to our knowledge, no such data is freely available.

Each examination e i is processed independently. For each one, we extract among all the 3D SLIC clusters those overlapping the annotated 2D axial slices. Among these clusters, we identify those whose intersection with annotated slices has a predominant GT label (consensus over at least 95% of voxels). This predominant label is assigned to the 3D cluster as its own GT label. Finally, 1/3 of these labeled 3D SLIC clusters are randomly selected (to simulate user interaction) and used to train our multi-phase cluster-wise random forest (MpCl -RF). Once trained, the forest is employed to classify all the remaining 3D SLIC clusters. We provide comparisons with single-phase voxel-wise (SpVx ), single-phase cluster-wise (SpCl ) and multi-phase voxel-wise (MpVx ) random forest (RF). To make comparisons possible, training for Sp/MpVx focus on voxels belonging to clusters selected for training in Sp/MpCl. Each forest contains T = 100 trees.

Features differ for each method. Cluster-wise strategies assign to clusters the features of Fig. 3 for MpCl and a modified set for SpCl since multi-phase features are ousted and spatial intensity and gradient ones only focus on I EV . In voxel-wise, MpVx assigns to voxels intrisic and spatially averaged (3 3 windows) intensity and gradient at each phase as well as PE, AUC, ∆ EV /AR and ∆ LV /EV whereas SpVx only involves intrisic and averaged intensity and gradient in I EV .

SpVx, SpCl, MpVx and MpCl are quantitatively assessed via TN rate error ∆τ = |τ -τ GT | comparing estimated and GT TN rates and DICE coefficients between obtained and fused GT masks for parenchyma, active and necrotic tissues: DICE prcm , DICE activ and DICE necro . To be less sensitive to the variability due to RF random aspects, results for each e i are averaged over 10 realizations.

Discussion. We present in Tab. 1 a comparative assessment of the 4 methods through ∆τ and DICE coefficients averaged over the whole database. Segmentation results are displayed on Fig. 4 with corresponding fused GT mask (one slice per examination). The comparative study reveals better results using MpCl -RF. It obtains the smallest ∆τ with 5.26 and the highest DICE for parenchyma, active and necrotic areas with 74.4, 71.9 and 93.3. In comparison, MpVx -RF reaches equivalent DICE necro but less accurate ∆τ , DICE activ and DICE prcm . According to visual results (Fig. 4), MpVx -RF allows an obvious worst spatial consistency while requiring more interaction efforts compared to MpCl -RF.

The comparison between single and multi-phase approaches confirms that exploiting multi-phase CT images instead of one single scan significantly improves the results with gains of DICE around 8.6 (11.2), 8.1 (6.8) and 3.6 (6.9) resp. for cluster and voxel-based strategies. The significant differences between Sp/MpCl -RF tend to justify the use of robust multi-phase features when making RF working on clusters since a classification error on one single cluster can have a strong impact on the accuracy. In addition, the averaged maximum inter-expert variability computed over the GT masks provided by the 4 experts is 21.2 ± 15.2 in TN rate. This significant variability reinforces the necessity to work toward an efficient computed-aided diagnosis method dedicated to TN rate estimation in order to overcome the bias induced by visual qualitative assessment.

Conclusion

In this work, we addressed the semi-automatic evaluation of pre-operative locoregional treatments (PLT) response for hepato-cellular carcinoma (HCC). Toward this goal, we proposed a method that estimates the tumor necrosis (TN) rate from dynamic contrast-enhanced CT scans. While ensuring weak interaction efforts for practitioners, it accurately segments parenchyma, active and necrotic tissues. Our approach applies random forest on supervoxels and involves robust multi-phase cluster-wise features. Quantitative assessment on clinical data confirms the benefits of exploiting dynamic information extracted from multi-phase images at a cluster spatial extent. Multi-examination learning would deserve further investigation to make the proposed strategy becoming fully automatic. More generally, it could be easily applied to other tumor types, organs and modalities.

Fig. 1 .

 1 Fig. 1. Registered multi-phase input data ( ĨWI, ĨAR, IEV and ĨLV) obtained from dynamic contrast-enhanced CT scans with associated liver segmentation mask.

Fig. 2 .

 2 Fig. 2. Semi-automatic tumor segmentation with multi-phase cluster-wise random forest: by performing user interaction (b) and learning on 3D supervoxels (a), it segments parenchyma (blue), active (red ) and necrotic (green) areas of the 3D liver volume (c).

Fig. 3 .

 3 Fig. 3. Multi-phase cluster-wise features involved into random forest (a) to discriminate the different tissues and illustration of dynamic features (b).

Fig. 4 .

 4 Fig. 4. HCC liver tumor segmentation results via single-phase voxel-wise (SpVx ), single-phase cluster-wise (SpCl ), multi-phase voxel-wise (MpVx ) and the proposed multi-phase cluster-wise (MpCl ) random forest (RF) with groundtruth (GT) masks. Parenchyma, active and necrotic areas are respectively in blue, red and green.

Table 1 .

 1 Quantitative comparisons of single-phase voxel-wise (SpVx ), single-phase cluster-wise (SpCl ), multi-phase voxel-wise (MpVx ) and the proposed multi-phase cluster-wise (MpCl ) random forest (RF) through TN rate error and DICE coefficients averaged over the whole database. Best results are emphasized in bold.

	methods	SpVx -RF SpCl -RF MpVx -RF MpCl -RF
	∆τ	6.40 ± 2.85 9.13 ± 4.78 6.60 ± 3.32 5.26 ± 3.90
	DICEactiv 54.3 ± 17.2 65.8 ± 15.3 65.5 ± 12.4 74.4 ± 12.6
	DICEnecro 65.0 ± 21.6 63.8 ± 25.8 71.8 ± 17.6 71.9 ± 19.5
	DICEprcm 80.5 ± 13.1 89.7 ± 4.90 87.4 ± 9.00 93.3 ± 3.08

Acknowledgments. This work received the financial support from Fondation Arc, www.fondation-arc.org.