Bernard Picinbono
email: bernard.picinbono@lss.supelec.fr

Life Fellow

Generating Singular Signals by Filtering Bernoulli Correlated Inputs

Keywords: Statistical signal analysis, Signal and noise modeling, Non-Gaussian signals and noise, Markov processes., Prediction

Singular random signals can be obtained at the output of some linear filters when the input is a Bernoulli white noise. It is shown that this whiteness assumption can be relaxed and some examples of colored Bernoulli signals generating singularity by filtering are presented. Computer experiments are realized in order to verify these results.

I. INTRODUCTION

Singular random signals have been recently presented and analyzed in [START_REF] Picinbono | Singular random signals[END_REF]. A signal x k is said to be singular if the random variables (RV) x k are singular for any k. This means that the marginal distribution function (DF) F (x) of x k is continuous but has a derivative equal to zero almost everywhere. Usually these RVs are considered as mathematical curiosities and one of the main purposes of [START_REF] Picinbono | Singular random signals[END_REF] was to show that, on the contrary, they are very common. In particular the output x k of a linear filter generated by the input signal w k is singular if the filter satisfies some simple conditions concerning its poles and if w k is discrete-valued and strictly white, or is a sequence of independent and identically distributed (IID) random variables (see p. 254 of [START_REF] Picinbono | Random Signals and Systems[END_REF]). This whiteness assumption is introduced in almost all the studies on this question [START_REF] Picinbono | Singular random signals[END_REF], [START_REF] Lukacs | Characteristic Functions[END_REF], [START_REF] Wittke | Infinite series of interference variables with Cantor-type distributions[END_REF], [START_REF] Naraghi-Pour | On the continuity of the stationary state distribution of DPCM[END_REF], [START_REF] Smith | Error probabilities on fading channels with intersymbol interference and noise[END_REF]. Our purpose is to relax this assumption and to show that some correlated inputs can also generate singular signals.

II. THEORETICAL ANALYSIS

The starting point of the analysis is the Lebesgue decomposition theorem. It says that the DF F (x) of a random variable (RV) X can be expressed in an unique way as F (x) = a 1 F c (x) + a 2 F d (x) + a 3 F s (x), a i ≥ 0, a 1 + a 2 + a 3 = 1.

(

) 1
In this equation the three functions F i (x) are DFs and F c (x), F d (x), and F s (x) are the continuous, discrete and singular components of F (x) respectively. If a 1 = 1, the RV is continuous and its PDF is the derivative of F c (x). If a 2 = 1, X is a discrete RV. If, finally, a 3 = 1, the RV X is singular.

The spectrum S F of a RV X is the set of the points of increase of its DF F (x). Its Lebesgue measure L(S F ) is called the spectral measure (SM) of X. If the SM is zero, the continuous part in (1) is zero, or a 1 = 0. Then the RV is either discrete (a 3 = 0), or singular (a 2 = 0 ), or a mixture of a discrete and a singular parts. As a consequence in order to show that a RV is singular it suffices to show that its SM is zero and that there is no discrete part in the Lebesgue decomposition of its DF.

Let now F be a causal filter and h k its impulse response. The output x k generated by the input

w k is x k = ∞ i=0 h i w k-i .
In order to simplify the presentation we shall study only symmetric Bernoulli input signals w k . This means that w k takes only the values ±1 with equal probabilities. The extension to signals with more than two discrete values introduces only complexities of notations. The symmetry assumption of w k implies that x k is also symmetric, which means that x k and -x k have the same DF which is the DF of the RV

X = ∞ k=0 h k w k . (2)
Let F S be the set of filters satisfying the following conditions:

(a) They are dynamical, which means that they are causal, stable, and that their transfer function H(z) is a rational function.

(b) They have an infinite impulse response (IIR), which means that there is at least one nonzero pole.

(c) Their poles are inside the circle of singularity introduced in [START_REF] Picinbono | Singular random signals[END_REF]. It has the same center as the unit circle and its radius is 1/2.

The exponential filter k k = a k belongs obviously to this set if 0 < |a| < 1/2. The main result used below is contained in the following theorem partially indicated in [START_REF] Picinbono | Singular random signals[END_REF].

Theorem 1. Let w k be the input of a filter F belonging to F S and x k the corresponding output. If w k is a Bernoulli symmetric signal then the SM of x k is zero.

Proof. The filter F is defined by its transfer function H(z) or its impulse response h k . Let F be the filter with the impulse response

g k = 2 k h k . It is obvious that its transfer function is G(z) = H(z/2).
Then if the poles of F are z i , those of F are 2z i . The assumption that the poles of F are inside the circle of singularity implies therefore that the poles of F are inside the unit circle, or that F is stable and then ∞ k=0 |g k | < +∞. Let X N and R N be the finite sum and the rest defined by

X N = N k=0 h k w k and R N = ∞ k=N +1 h k w k . Since w k = ±1 we have |R N | ≤ ρ N with ρ N = ∞ k=N +1 |h k |.
As X N can take at the maximum 2 N +1 values, the SM S of X satisfies S ≤ 2 N +1 ρ N . This is valid for all N . Then S ≤ 2 lim N →∞ (2 N ρ N ). But we have

2 N ρ N = 2 N ∞ k=N +1 |h k | < ∞ k=N +1 2 k |h k | = ∞ k=N +1 |g k |, (3)
and the limit for N → ∞ is 0 because ∞ k=0 |g k | < +∞. This implies that the SM of X is 0, or a 1 = 0. It is shown in [START_REF] Picinbono | Singular random signals[END_REF] that if the input signal is white, which means that the w k are IID symmetric Bernoulli RVs, then there is no discrete component in the DF of X, or a 2 = 0, and this implies that X is singular.

Our purpose is to investigate whether, (and if so, under what conditions), the whiteness assumption can be relaxed. We assume then that the w k s are still symmetric Bernoulli, but not necessarily independent. Let X N be the partial sum defined by

X N = N k=0 h k w k . ( 4
)
It takes at the maximum 2 N +1 distinct values v N i . From these values we can construct a tree of representation analyzed in [START_REF] Picinbono | Singular random signals[END_REF]. To each value v N i at the step N we associate a node V N i of this tree and when passing from N to N + 1 each node V N i generates two distinct nodes defined by . The nodes V 0 1 and V 0 2 defined by -v 0 1 = v 0 2 = |h 0 | are generated from an origin node V which does not correspond DRAFT August 15, 2006 to a value of X N . A node V N i is said to be single if it is generated by only one node V N -1 j(i) . In the opposite case it is multiple. If all the nodes are single the number of distinct nodes V N i is 2 N +1 , and we assume that this is satisfied in what follows. This assumption is verified for the exponential filter where h k = a k when 0 < |a| < 1/2 and also for all the filters such that there is no crossing between the branches of the tree, as analyzed in [START_REF] Picinbono | Singular random signals[END_REF].

This assumption of single nodes depends uniquely on the filter. The set of filters F having this property is called F D . Most of the following results are related to filters belonging to the set

F SD = F S ∩ F D . The exponential filter h k = a k belongs obviously to this set if 0 < |a| < 1/2.
The assumption of single nodes means in particular that for any N and i there is a unique path going

from V to V N i . Let i N k (i), 0 ≤ k ≤ N
, be the indices j defining the nodes V k j of this path. These nodes can then be written V k i N k (i) and the index i N k (i) characterizes the unique node of the tree at the step k located on the path going from V to V N i . The problem is to calculate the probabilities

p N (i) = P [x N = v N i ] , 1 ≤ i ≤ 2 N +1 . (5)
When the w k s are IID this probability is 1/2 N +1 . When they are no longer independent, its calculation is much more complicated.

For this we introduce the conditional probability

p N (i, j) = p N [x N = v N i |x N -1 = v N -1 j ], (6)
called also transition probability in the tree. It has two fundamental properties for the analysis that follows.

The first comes from its normalization specified for all j by 2 N +1 i=1 p(i, j) = 1. However a node V N -1 j of the tree of construction generates only two nodes V N i characterized by the indices i + (j) and i -(j) and, according to (4), corresponding to the values

v N i± (j) = v N -1 j ± h N .
As a consequence for a given j there are only two terms in the previous sum and we have

p N [i + (j), j] + p N [i -(j), j)] = 1. (7)
The second starts from the fact that for any node V N i there is only one node V N -1 j(i) at the step N -1 of the tree generating V N i at the step N . Thus p N (i, j) is zero except when j = j(i), and the only nonzero values of p N (i, j) are

q N (i) = p N [i, j(i)] (8)
for N > 0 and q 0 (i) = 1/2.

It results from [START_REF] Smith | Error probabilities on fading channels with intersymbol interference and noise[END_REF] and from the unicity of the path between V and

V N i that p N [(x N = v N i ) ∩ (x N -1 = v N -1 j )] = p N -1 (j) . p N (i, j) δ[j -j(i)], (9)
where δ[k] is the Kronecker delta symbol equal to 1 if k = 0 and 0 otherwise. By a summation on j, which contains only one term, we obtain

p N (i) = p N -1 [j(i)] q N (i). (10)
By repeating this at all the nodes of the unique path between V and V N i characterized by the indices

i N k (i) we obtain p N (i) = N k=0 q k [i N k (i)]. (11)
When the RVs w k are IID we have of course q N [i N k (i)] = 1/2, and we find again that the values

v N i have equal probabilities 1/2 N +1 .
The probabilities p N (i) of (11) are normalized, or i p N (i) = 1, where the sum is extended to all the indices i from 1 to 2 N +1 . This property is valid for N = 0 because q 0 (i) = 1/2. Suppose that it is valid at the step N -1. Since each node V N -1 j generates only two nodes V N i-(j) and V N i+(j) the result comes from (7).

The relation (11) is the basis of the following analysis of the singularity. Indeed if the SM of X is zero and if all the p N (i) tend to 0 when N → ∞, there is no value v ∞ i with a finite probability, and this means that the RV X cannot have a discrete component and then is singular. This can be specified by the following theorem.

Theorem 2. Let X be the RV defined by (2) where h k is the impulse response of a filter belonging to F SD and w k a sequence of Bernoulli RVs. If the transition probabilities q N (i) defined by (8) satisfy

0 < q N (i) < B < 1, (12)
then the RV X is singular, or a 1 = a 2 = 0.

Proof. As F belongs to F SD the SM of X is zero, according to Theorem 1. It remains to show that a 2 = 0. This is a direct consequence of (11) and (12) because p N (i) < (1/2)B N , which tends to zero when N → ∞.

Comments. Note that this theorem introduces only a sufficient condition of singularity. Note also that it can be applied for white input because in this case q

N (i) = 1/2. The condition F ∈ F SD is satisfied DRAFT August 15, 2006
by a large class of filters and some are presented in [START_REF] Picinbono | Singular random signals[END_REF]. However the question of characterizing all the dynamical filters belonging to F SD remains open. As a matter of fact it is possible to extend this theorem to the case where some nodes are not single, but this introduces other conditions that cannot be presented in this paper.

III. EXAMPLES

Application of condition (12) requires the calculation of the probabilities q N (i) defined by (8) and related to the transitions probabilities p N (i, j) defined by [START_REF] Smith | Error probabilities on fading channels with intersymbol interference and noise[END_REF]. For this it is interesting to express these probabilities in terms of those of the RVs w k . For this note first that for a given v N -1 j , or a given node

V N -1 j
, there is only one path coming from

V to V N -1 j , or one sequence of values of w k , 0 < k ≤ N -1.
Let us denote this sequence as S N -1 (j). The transition probability p N (i, j) is simply related to the conditional probabilities of the w k by

p N [i ±1 (j), j] = P [w N = ±1|S N -1 (j)], (13)
where i + (j) and i -(j) are the two indices i defining the two nodes V N i generated by the node V N -1 j .

Note that the probabilities (13) are now transition probabilities in time.

These expressions are especially interesting in the case where the input w k is a Bernoulli symmetric Markovian signal. Let us remind that a symmetric Bernoulli Markovian signal of order P is a signal taking the values ±1 with equal probabilities and which can be expressed as The Markov assumption means that the sequence S N -1 (j) in (13) can be replaced by a sequence using only the past of order P . This can be written

w k = f (w k-1 , w k-2 , ..., w k-P , b k ), (14)
p N [i ± (j), j] = P [w N = ±1|w N -1 , ..., w N -P ]. (15)
Since the RVs w k take only the values ±1, there are only at most 2 P distinct values of p N [i ± (j), j]. If we assume that the signal is not predictable there is no value of p N equal to 1. Then the bound of (12)

August 15, 2006 DRAFT is simply the maximum value of the finite number of transition probabilities p N (i, j). Then Theorem 2 can be applied and the RV X is singular.

It is interesting to present examples of signals for which Theorem 1 can be applied but not Theorem 2. In this case the spectral measure is zero, or a 1 = 0, but the RV X is not purely singular, or a 3 = 1.

We shall first present a case where it is purely discrete, or a 2 = 1.

This especially appears when the symmetric Bernoulli signal w k is predictable. A predictable signal of order P is defined by ( 14) where b k = 0. This means that at each k, w k can be estimated without error in terms of its past P values. As a consequence w k is also completely defined from its P first values, or

w k = h(w 0 , w 1 , ..., w P -1 ). ( 16
)
The simplest example of such a Bernoulli signal is w k = w 0 (-1) k , where w 0 is a Bernoulli symmetric RV. Here P = 1.

It is clear that if w k is predictable, then the RV X of ( 2) is discrete, or a 2 = 1. Indeed since the RVs w i take only 2 values, X takes at the maximum 2 P distinct positive values which means that it is a discrete RV.

In the case of the signal w k = w 0 (-1) k used with the exponential filter defined by h k = a k , we have X = w 0 /(1 + a) and it takes only the two values ±1/(1 + a) with equal probabilities.

It is clear that conditions of Theorem 2 are not satisfied for predictable signals because most of the transitions probabilities p N (i, j) are equal to 1 or 0 in such a way that there is only a finite number of paths going from V to infinity and giving values of X with nonzero probabilities.

Finally it is of interest to note that conditions of Theorem 1 can lead to a RV X with a DF mixture of discrete and singular part, or a 1 = 0 and 0 < a 2 < 1. We shall present in the next section an example of such a situation.

IV. EXPERIMENTS

Our purpose is now to verify the previous theoretical results by some computer experiments. The principle of these experiments is the same as the one used in [START_REF] Picinbono | Singular random signals[END_REF]. In order to analyze the singularity of a signal x k we realize histograms at different scales of a very large number of its values. The singularity is illustrated by the fact that the histograms introduce a fractal structure. This means that they have the same form whatever the scale of analysis. This shows, with the limited precision of any experiment, the lack of derivative.

In order to verify results of the previous section we shall introduce a procedure to generate by computer the values 0 or 1 with probabilities 1 -p and p respectively. Similarly let v k be a Bernoulli strictly white noise taking the values ±1 with the same probabilities 1/2 and independent of u k . Consider now the signal

w k = u k w k-1 + ūk v k , , k ≥ 1, (17)
where ūk = 1 -u k . It has the form (14) where

P = 1 and b k = [u k , v k ].
It is obvious that if w k-1 takes only the values ±1 with probabilities 1/2, w k has the same property.

It can be shown that, whatever the values of w 1 , w k tends to have this property and that the covariance function γ k of w k is p |k| , so it is an exponential covariance function. For p = 0, u k = 0, and w k is the white noise v k with a zero covariance function. On the other hand if p = 1, u k = 1 and w k = w 0 , which introduces constant covariance function or long range memory.

In the experiments presented below we use an exponential filter often mentioned in the previous analysis and also used in [START_REF] Picinbono | Singular random signals[END_REF]. Its input-output relationship can be written in a recursive form with the recursion

x k = ax k-1 + w k
, where w k is given by (17). In Fig. 1 we present histograms obtained when the input w k is white, or when p = 0. The experiment is realized with 10 7 samples and a = 1/3. The fractal structure appears very clearly. In Fig. 2 we present the same histograms realized with correlated inputs defined by the linear Markovian model (17) with p = 0.5. The fractal structure remains valid but the main difference is the fact that the number of samples recorded in each cell of the histograms have much stronger variations than in the case of white inputs. This can be explained theoretically but the analysis cannot be presented here.

Note that when the input is white the histograms are symmetric with respect to their centers. This is no longer true when Markovian inputs are used, except for the histogram centered at 0 because X is symmetric and its histograms are symmetric with respect to 0. This is why only histograms of positive values of x k are represented.

We shall now consider experiments in which the DF of the output x k is a mixture of discrete and singular parts. For this consider again the signal u k = u 0 (-1) k and a white signal v k taking the values ±1 with equal probabilities and independent of u 0 . Suppose now that w k is equal either to u k or to v k for all k with the probabilities α and 1 -α respectively. This implies that its covariance function is

γ k = α(-1) k + (1 -α)δ[k]. The DF of the RVs X is F (x) = αF d (x) + (1 -α)F s (x).
Applying this signal again at the input to the exponential filter h k = a k with a < 1/2 yields an output signal taking the values ±1/(1 + a) with probability α and being singular with probability 1 -α.

Results of computer experiments on this signal are presented in Figures 3 and4. In these experiments a = 1/3, as in Figs. 1 and2, and α = 0.5. This implies that x k takes the values ±0.75 with probability 0.25, which means that its DF is discontinuous for x = ±0.75 and the amplitude of the discontinuity is 0.25. The first histogram of Fig. 3 exhibits these two discontinuities characterized by two lines at the points ±0.75. One of these lines disappears in the second histogram of this figure and there is no longer a line in the last two histograms. The fractal structure introducing a singular component appears clearly in these histograms. In order to verify whether or not the discontinuity of the DF is an experimental artefact, we present in Fig. 4 four histograms isolating the point 0.75 with cells of decreasing widths. In the first histogram we observe still a residual contribution of the singular part of the DF. There is however a line at 0.75 and the fact that its amplitude is constant indicates clearly that there is effectively a discrete component. As the experiment uses 10 6 samples, the discontinuity corresponding to the probability 0.25 must be of the order of 2.5.10 5 , which clearly appears in the four histograms. This shows that the signal

x k is effectively a mixture of a discrete-valued signal and a singular signal.

V. CONCLUSION

When the input of a linear filter satisfying some conditions concerning essentially the location of its poles with respect to the circle of singularity is a Bernoulli white noise, the output can be singular.

This assumption of whiteness can however be partially relaxed while maintaining the singularity of the output. Some sufficient conditions ensuring the singularity with colored inputs have been established.

These conditions are obviously satisfied not only by white noise but also by a large class of correlated signals. It is especially the case of Markovian signals of finite order. The theoretical analysis also shows that the output generated by colored inputs can be a mixture of a discrete and a singular distribution. In particular if the input signal is predictable it was shown that the output of the filter is simply a discrete random variable. Computer experiments in order to verify the theoretical results have been realized and discussed. For this purpose a specific model of linear Markovian signal of order one was introduced and the experimental results are in complete agreement with the theory. Finally a model of Bernoulli input ensuring that the output contains a discrete and a singular part was introduced and here also the experimental results are in perfect agreement with the theory.

Figures captions

Fig.

  where b k is a vectorial white noise. It is clear that the conditional DF of w k at time k, conditionally to the whole past, depends only on w k-1 , w k-2 , ..., w k-P . The function f (.) and the noise b k must satisfy conditions ensuring that if the w l s, k -P ≤ l ≤ k -1 are symmetric Bernoulli signals, w k has the same property. We shall present in the next Section examples of such signals that can be realized in computer experiments.

1 .

 1 Fig. 1. Histograms of x k at different scales. a = 1/3; p = 0.

Fig. 2 .

 2 Fig. 2. Histograms of x at different scales. a = 1/3; p = 0.5.

Fig. 3 .

 3 Fig. 3. Mixture discrete-singular, α = 0.5. Fractal part of the histogram.

Fig. 4 .

 4 Fig. 4. Mixture discrete-singular, α = 0.5. Analysis of the neighborhood of 0.75.

August 15, 2006 DRAFT


Aymen Ben Amor

M Fakhfakh

C Renard

Stefan Dilhaire

Laetitia Vincent

Stéphane Grauby

Thermal conductivity reduction measurement on Ge and Si NWs: heterostructuration

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



Nanostructured material Bulk

𝑆: Seebeck coefficient (V/K), σ: electrical conductivity ( -1 m -1 ), 𝜅 : thermal conductivity (Wm -1 K -1 ).

Typical value for bulk material: 1 at best, Objective: at least 3.

Purpose of this study

), 2 ( 2 ) 3 (    T I R V Tip Tip   𝛼 𝑇𝑖𝑝 : tip temperature coefficient (K -1 )

Experimental protocole

Thermal conductivity of heterostructured nanowires (NWs)

Conclusions and prospects

 Diameter and heterostructuration lead to a thermal conductivity reduction in Si and Ge nanowires,  In the future: -Confirmation on new samples, in particular at 1100°C (more heterostructured domains with larger size), -Development of new experimental set-ups to measure 𝜎 and S in order to estimate 𝑍𝑇.

a

  Université Bordeaux, LOMA, UMR CNRS 5798, 351 cours de la libération, 33405 Talence cedex, France b Université Paris-Sud, C2N, Bât 220, rue André Ampère, 91405 Orsay cedex, France Thermal conductivity reduction measurement on Ge and Si NWs: heterostructuration and diameter size effect  Measuring the thermal conductivity 𝜅 of heterostructured nanowires,  Evaluating the influence of diameter and heterostructuration on 𝜅. * e-mail:aymen.ben-amor@etu.u-bordeaux.fr  Thermoelectric materials: conversion of thermal energy into electrical energy and vice versa (nanoheater, nanocooler, power generator)  Thermoelectric efficiency evaluated by the factor of merit:  Interest of nanostructuration: dimension reduction (< phonon mean free path)  phonon boundary scattering and group velocity reduction  𝜅 reduction  𝑍𝑇 increase  Double challenge:-find the best material with high σ and low 𝜅 : phonon glass and electron crystal  semiconducting materials, -measuring 𝜅 at nanometric scale  𝐴𝐹𝑀 with thermoresistive tip in active regime: 3-Scanning Thermal Microscopy (𝑺𝑻𝒉𝑴).



  Amor a* , M. Fakhfakh b , C. Renard b , S. Dilhaire a , L.Vincent b , S. Grauby a 𝜅 limited by phonon-phonon And phonon-impurity interactions Bordeaux -Juillet 2017  3-Scanning Thermal Microscopy (𝑆𝑇ℎ𝑀) : High 𝜅  high heat flux  low 𝑇(2𝜔)  low 𝑉(3𝜔)  Simultaneously topographical and 3𝜔 -𝑆𝑇ℎ𝑀 thermal conductivity contrast images  Only method to measure the thermal conductivity of individual nanowires (NWs) embedded in their matrix [

 Heterostructured germanium (Ge) NWs [3] Equivalent thermal resistance R eq (K/W)


Lou Gevaux

Cyprien Adnet

Pierre Séroul

Raphael Clerc

Alain Trémeau

Jean Luc Perrot

Mathieu Hébert

Three-dimensional hyperspectral imaging: a new method for human face acquisition

Hyperspectral imaging is an emerging non-invasive method for the optical characterization of human skin, allowing detailed surface measurement over a large area. By providing the spectral reflectance in each pixel, it enables not only color simulation under various lighting conditions, but also the estimation of skin structure and composition. These parameters, which can be correlated to a person's health, are deduced from the spectral reflectance in each pixel thanks to optical models inversion. Such techniques are already available in 2D images for flat skin areas, but extending them to 3D is crucial to address large scale and complex shapes as in the human face. The requirements for accurate acquisition are a short acquisition time for in vivo applications and uniform lighting conditions to avoid shadowing. The proposed method combines wide field hyperspectral imaging and 3D surface acquisition, using a single camera, with an acquisition time of less than 5 seconds. Complete color consistency can be achieved by computationally correcting irradiance non-uniformities using 3D shape information.

INTRODUCTION

Skin is a complex layered organ composed of several types of cells, fibers, glands, vessels, pigments and blood, which together generate a heterogeneous distribution of light absorption and scattering [START_REF] Igarashi | The appearance of human skin: A survey[END_REF]. The measurement of skin optical properties can offer crucial information about a person's health, leading to a growing demand for optical devices in the fields of cosmetology and dermatology. In recent years, research in the measurement of skin optical properties has focused on non-invasive in vivo methods as ex vivo methods are known to alter the properties of the skin, and as invasive analysis may cause patient discomfort. Several optical methods have emerged, such as diffuse reflectance spectroscopy [START_REF] Zonios | Skin melanin, hemoglobin, and light scattering properties can be quantitatively assessed in vivo using diffuse reflectance spectroscopy[END_REF], modulated imaging [START_REF] Cuccia | Quantitation and mapping of tissue optical properties using modulated imaging[END_REF], spectral imaging, and color imaging [START_REF] Doi | Estimation of Local Skin Properties from Spectral Images and its Application to Appearance Reproduction[END_REF] [START_REF] Seroul | Model-based skin pigment cartography by high-resolution hyperspectral imaging[END_REF]. Data retrieved from these methods have been analyzed using various models such as the Monte Carlo method [START_REF] Nishidate | Estimation of melanin and hemoglobin in skin tissue using multiple regression analysis aided by Monte Carlo simulation[END_REF], the radiative transfer equation [START_REF] Magnain | Skin color modeling using the radiative transfer equation solved by the auxiliary function (a) (b) (c) IS&T International Symposium on Electronic Imaging 2018 Material Appearance 2018 method: inverse problem[END_REF], the diffusion approximation [START_REF] Cuccia | Quantitation and mapping of tissue optical properties using modulated imaging[END_REF], flux models such as the Kubelka Munk model [START_REF] Seroul | Model-based skin pigment cartography by high-resolution hyperspectral imaging[END_REF], and independent component analysis [START_REF] Tsumura | Independentcomponent analysis of skin color image[END_REF].

Hyperspectral imaging [START_REF] Lu | Medical hyperspectral imaging: a review[END_REF] is a non-invasive and contactless method that provides more information than standard color imaging and may reveal details invisible to the human eye. Moreover, it can provide information related to color appearance under various illuminants as well as physicochemical composition. By analyzing the acquired spectral reflectance distribution over a large area using an optical model combined with a skin structure model, various skin composition properties can be estimated pixel per pixel, then displayed in the form of images or maps. Since the capacity of this technique to reveal fine details is directly related to the image resolution, an imaging system with high spatial and spectral resolution is preferable. However, as the computation time needed to analyze the spectral information increases proportionally to the number of pixels, the skin model and the optical model chosen should be as simple as possible, to make the time needed for data analysis acceptable. In a previous work [START_REF] Seroul | Model-based skin pigment cartography by high-resolution hyperspectral imaging[END_REF], a Kubelka-Munk inspired model applied onto a two-layer skin structure was shown to give satisfactory results by computing up to five skin composition maps of 1120 × 900 pixels in less than one hour. The hyperspectral images were captured by using the SpectraCam® camera [START_REF] Seroul | Model-based skin pigment cartography by high-resolution hyperspectral imaging[END_REF] [START_REF] Nkengne | SpectraCam®: A new polarized hyperspectral imaging system for repeatable and reproducible in vivo skin quantification of melanin, total hemoglobin, and oxygen saturation[END_REF] developed by Newtone Technologies. This technique, however, is limited to small planar skin areas, as the effect of shadows on complex three-dimensional objects is not taken into account. The combination of this technique with 3D images of the observed skin area thus appears as a solution to surpass the aforementioned limitations and address large, non-flat areas such as a human face.

The few existing near field 3D-spectral imaging systems do not seem really appropriate to high spectral and spatial measurement on living organs. Either they are multispectral systems (e.g. [START_REF] Paquit | 3D and multispectral imaging for subcutaneous veins detection[END_REF] [12] [START_REF] Zhang | A novel 3D multispectral vision system based on filter wheel cameras[END_REF]) providing too low spectral resolutions for optical analysis, or in the case of those 3D-spectral imaging systems that do provide high spectral resolutions (e.g. [START_REF] Kim | 3D Graphics Techniques for Capturing and Inspecting Hyperspectral Appearance[END_REF]), the acquisition process is not adapted to living subjects. This is the reason why we designed the 3D-spectral measurement setup dedicated to skin presented in this paper, which gives a satisfactory tradeoff between high resolution and short acquisition time. In addition to this device, solutions for ensuring good color consistency and reducing shadowing effects are investigated. We focus here on the human face, the part of the body that presents the biggest challenge for 3D scanning, as well as having the highest interest for cosmetology, to create a robust device that is capable of scanning every part of the human body.

The paper is structured as follows: the measurement setup is presented first in Section 2. In the following, the 3D acquisition (Section 3), the hyperspectral acquisition (Section 4) and the optical treatment of the measured data (Section 5) are introduced, before presenting experimental results in Section 6. The question of irradiance correction is briefly discussed in Section 7, before reaching conclusions.

EXPERIMENTAL SETUP

The proposed 3D hyperspectral acquisition system is an extension of the existing hyperspectral camera, SpectraCam® [START_REF] Seroul | Model-based skin pigment cartography by high-resolution hyperspectral imaging[END_REF] [10], incorporating the acquisition of depth map information and allowing wide field acquisition. It is composed of a monochrome CMOS camera of 2048 × 2048 pixels and a Liquid Cristal Tunable Filter (LCTF).

The setup comprises two lighting systems, used sequentially: one for the depth map acquisition, one for the hyperspectral acquisition. The first lighting system is a portable digital LED projector used to project fringes during the 3D scanning step of the acquisition process. The angle between the camera and the projector is approximately 10 degrees, to enable a full face acquisition while minimizing the occluded areas. The second lighting system comprises two LED lighting units designed with blue and white LEDs to provide sufficient irradiance on the skin over the whole visible spectrum.

Linear polarizing filters are added on each lighting unit. They are oriented in order to create a cross polarized (CP) configuration with the LCTF polarization direction, this latter playing the role of an analyzer for the reflected light. The CP configuration is used to remove the gloss component of the reflected light: when light interacts with a diffusing medium, a portion of the flux is reflected at the air-medium interface according to the Fresnel law, without any change in polarization. This specular reflection of light, related to the visual attribute called gloss, is removed by the analyzer in a CP configuration. The remaining portion of flux enters the medium where it is scattered, progressively losing its polarization state, and finally emerging from the medium as an unpolarized light that can partially cross the analyzer. In the rest of the paper, the specularly reflected light will be discarded and only the unpolarized light component will be considered.

The camera, LCTF and lighting are controlled with a micro controller in order to shorten acquisition time. The acquisition is driven with a dedicated software written in C++.

A plane is also needed for reference and calibration. This plane is positioned perpendicularly to the camera's optical axis; together they define the reference coordinates system (X, Y, Z) shown in Figure 1, where Z corresponds to the optical axis. The camera and projector pupils are located in the same plane Z = L. The reference plane is mounted on a translation system which allows precise translations along the Z axis.

Hyperspectral images over the visible spectrum (410 nm -700 nm) and data for 3D reconstruction are sequentially acquired within five seconds, a duration that does not totally obviate the risk of the subject moving, but remains fairly acceptable.

THREE-DIMENTIONAL SCANNING

Various optical and imaging technologies provide non-contact methods for collecting the geometrical shape of an object [START_REF] Remondino | Image-based 3D modelling: A review[END_REF]. These methods can be categorized into those that are passive and active.

Passive methods don't require any controlled illumination and rely on ambient light. Stereoscopic imaging, inspired by binocular vision, is one of the most commonly employed 3D passive techniques: by capturing a scene from at least two points of view, the 3D geometry can be estimated using the triangulation principle. A given point has to be located in both views as a first step of the 3D reconstruction. Therefore, binocular and multi-view stereovision consist in solving a correspondence problem.

Active methods are based on the same triangulation problem. In contrast with stereovision, one camera is replaced with a controllable illumination source that projects a spot, a line or a structured light onto the object. Single point laser scanning emerged in the 1970s, allowing a 3D acquisition without any ambiguity. A similar concept using a sweeping line remains a popular solution for 3D acquisition of static scenes or objects that can be placed on a turntable. However, as the laser scanning process is slow and limited to static scenes, structured light projection techniques are now widely used [16] [17].

The use of a digital projector to project structured light into a full field-of-view replaces the temporal sweeping of the laser beam. This method is faster but not ambiguity-free. Many suggestions have been proposed to overcome the correspondence problem: for example, coded pattern projection (known geometries of spots, color coded fringes, binary coded fringes, etc.) or sinusoidal fringe projection (Fig. 2) with phase unwrapping algorithms [START_REF] Geng | Structured-light 3D surface imaging: a tutorial[END_REF].

The scanning system that we developed is based on a fringe projection method, offering an excellent compromise between speed, quality and cost. This method has already proven its effectiveness in biomedical applications such as 3D intra-oral dental measurements, human body shape measurements, body deformation measurements and in vivo measurements [START_REF] Wang | 3D absolute shape measurement of live rabbit hearts with a superfast twofrequency phase-shifting technique[END_REF]. Structured light projection also allows the use of a single camera for both hyperspectral and 3D acquisition, which ensures a good pixel correspondence between the 3D shape and the texture.

Triangulation principle

Sinusoid fringes are projected onto an object and deformed in function of the object geometry. The measurement of the phase deformation yields a depth map of the scene from the triangulation principle [START_REF] Lanman | Build your own 3D scanner: 3D photography for beginners[END_REF].

Once the phase deformation is known, the triangulation principle is applied to retrieve the depth information between the object and the reference plane. Let us consider a perfect projector whose pupil is centered in a point E, and a perfect camera whose pupil is centered in a point D (Fig. 3). Let us also consider a reference plane where the optical axes of the camera and the projector intercept in points G and C, respectively. Let us now consider a point F on the object, separated from the reference plane by a length h called the depth. Point F is illuminated by a ray from the projector that would intercept the reference plane in a point A were it not for the presence of the object; point F is also viewed by the same pixel (x,y) in the camera as the point B on the reference plane. According to the pinhole model, the depth h satisfies the following equation: By denoting the distance DG between the reference plane and the camera's pupil as L; the distance DE between the camera and projector pupils as b; and the fringe frequency on the reference plane as f, the phase difference between points A and B is:

AB DE DG h h   (1)
obj ref     
where obj  and ref  are respectively the phase at point F and at point B.

From equation ( 1), it is possible to write the relationship between the depth and the phase difference:

2 1 L h fb     (2)

Phase shift principle

In order to get information on each pixel (x,y), the phase shift method is used. N fringe images, defined by the following equation, are projected:

2 ( 1) ' ''cos 2 , n n I I I fx N            1, n N  (3)
where I' is the average intensity and I" is the intensity modulation.

The projected images represent identical fringes, each of which are successively shifted by 2π/N. On the projected images, the phase, 2πfx, depends on the pixel location x. On the images Jn, (n=1,… N) recorded by the camera, the phase is modified and can be computed using:

    2 ( 1) 1 2 ( 1) 1 sin arctan [ , ] cos N n n N n N n n N n J J                 (4)
The selection of the N value depends on the application requirements. When N increases, the acquisition time also increases but noise is reduced [START_REF] Li | Gamma-distorted fringe image modeling and accurate gamma correction for fast phase measuring profilometry[END_REF]. The observed noise, a consequence of the non-linearity of the projector-camera grey scale, is specific to fringe projection 3D acquisition: it is periodic, with a period of N times the period of the projected fringes. A N value around 5 or 6 is a good compromise between reduced noise (which is already much lower than with N = 3) and acquisition time.

Phase unwrapping algorithm

The phase calculated using ( 4) is discontinuous since its values are restricted in the interval [-π, π]. The continuous phase is retrieved through a phase unwrapping algorithm. Phase unwrapping techniques can be classified into two categories: temporal phase unwrapping [START_REF] Zuo | Temporal phase unwrapping algorithms for fringe projection profilometry: A comparative review[END_REF] and spatial phase unwrapping. Temporal phase unwrapping relies on the projection of a large number of images (binary coded fringes, multiple frequency acquisition…) to acquire the necessary information to solve the unwrapping problem. Spatial phase unwrapping uses the continuity between pixels to retrieve the original phase. Therefore, an intuitive method to correct spatial discontinuities is to assume that the phase difference between two neighboring pixels is always less than π: the phase has to be unwrapped each time the difference between two neighboring pixels is higher than π. However, this method is not robust to noise and errors propagate, necessitating more complex solutions. Since a temporal method would increase the acquisition time, the chosen solution is Ghilia's "2D unweighted phase unwrapping" [START_REF] Ghiglia | Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods[END_REF] based on the discrete cosine transform (DCT), which is a spatial method sufficiently robust to noise for facial acquisition (Fig. 4).

Calibrations

The acquisition process comprises three calibration steps for 3D scanning: a radiometric calibration which has to be done prior to generating the projected fringes, a depth calibration which allows the conversion of the phase onto a depth map, and an XY-plane calibration which allows the conversion of the depth map into (x,y,z) coordinates expressed in millimeters.

Radiometric calibration

The projected fringes must have a perfectly sinusoidal profile for optimal 3D reconstruction quality. However, the sinusoidal profile defined in the digital image that is projected is often transformed into another profile due to the radiometric nonlinearity of the projector and the camera. The radiometric calibration consists in measuring and correcting the grey scale non-linearity of the projector-camera system. Various methods of radiometric calibration have been proposed [START_REF] Li | Gamma-distorted fringe image modeling and accurate gamma correction for fast phase measuring profilometry[END_REF]; the method chosen consists in creating corrected fringes to be projected. These fringes are not sinusoidal but designed in order to correspond to the chosen shape when seen through the camera after projection. To computationally correct the projected fringes, 20 grey levels are projected onto the reference plane. The best polynomial of degree two matching with the experimental curve is identified using a least-squares algorithm.

Finally, a lookup table is created to link each targeted grey level of an 8-bit scale with the value that has to be input before projection. When projected onto skin or any turbid media, light diffusion also affects the shape of the projected fringes. Since estimating this phenomena is complex, the 3D scanning is done, when possible, using a wavelength for which the media is as opaque as possible [START_REF] Wang | 3D absolute shape measurement of live rabbit hearts with a superfast twofrequency phase-shifting technique[END_REF]. As the opacity of skin is highest in the short visible wavelengths, we have selected the blue channels for our imaging systems.

Depth calibration

The conversion of phases into depth information from the triangulation principle described by (2) requires knowledge of the system geometrical parameters L, b and f. These parameters, in addition to being dependent on the pixel location for the fringe frequency f, are difficult to measure accurately. Therefore, a calibration must be performed. The simplest calibration is linear, based on the approximate assumption that the phase difference is directly proportional to the depth when L is much larger than h. A non-linear calibration which does not rely on this approximation is also possible. The implemented non-linear calibration, described by Jia [START_REF] Jia | Comparison of linear and nonlinear calibration methods for phase-measuring profilometry[END_REF], requires the acquisition of several phase maps on planes located at several positions (h1, h2 …) along the Z axis (Figure 5) to find two coefficients M and N such as:

( , ) ( , ) ( , ) ( , ) ( , ) xy h x y M x y N x y x y     (5)
A minimum of three acquisitions is required to identify the coefficients, but more measurements can be performed to give a better estimate of M and N by applying a least-squares algorithm.

XY-plane calibration

The depth calibration yields the coefficients M and N in mm -1 that are used in equation ( 5) to compute the depth map. The XYplane calibration step consists in converting the depth map into a point cloud that represents the object surface in the real world coordinates system, which can also be described as the "2.5D" information.

Several methods [START_REF] Zhang | Novel method for structured light system calibration[END_REF] [25] have been developed to fully characterize the camera and projector intrinsic parameters (focal length, pixel sizes, and optical distortion) in order to perform an accurate XY-plane conversion. The method that we implemented is basic and assumes that the zone of interest is not affected by optical distortion and that the acquired depth range is sufficiently small to have the same pixel size for the whole object. A ruler is imaged in order to obtain the pixel-to-millimeter ratio.

HYPERSPECTRAL IMAGING

Hyperspectral imaging is an extension of classical RGB imaging where the number of channels is much higher, corresponding to a number of narrow spectral bands continuously covering the visible spectrum (and possibly UV and IR domains). Hyperspectral image data are also called "hypercubes" since they are three dimensional: two dimensions for the spatial domain, one dimension for the spectral domain. The information captured in each pixel is proportional to the spectral radiance coming from one direction in the camera's field. By dividing the spectral radiance measured from the object by the spectral radiance measured from a perfectly white diffuser that is illuminated in the same way as the object, one obtains in each pixel its spectral reflectance factor. This latter is often considered as a spectral reflectance, which is theoretically true only when the object is a Lambertian reflector.

Many hyperspectral imaging systems have been developed based on various technologies and physical principles, with different specificities in terms of spatial resolution, spectral resolution, surface shape constraints and acquisition time [START_REF] Lu | Medical hyperspectral imaging: a review[END_REF]. The staring approach, which is the method implemented in the proposed system, records the full spatial information for each single spectral band and requires a temporal scan of the spectral dimension. It relies either on active illumination or on adding a filtering system in front of the sensor, which can be a filtering wheel, a tunable filter or an interferometer. Active illumination methods seem promising [START_REF] Goel | HyperCam: hyperspectral imaging for ubiquitous computing applications[END_REF] in terms of cost and acquisition time, given the availability of a wide range of LEDs with different emission spectra. However, since we are focusing on the face, consistent lighting over a large area is important, which seems difficult to achieve with active illumination methods. Therefore, we prefer a method using filters.

The chosen method is based on a Liquid Cristal Tunable Filter (LCTF) which offers the possibility of a rapid filtering scan within the visible spectrum with a sufficiently high spectral resolution for skin study. This kind of filter is based on the association of several Lyot filters [START_REF] Gat | Imaging spectroscopy using tunable filters: a review[END_REF], and is equivalent to high quality interference filters with a band width of 10 nm and a precision of 1 nm, between 400 and 720 nm.

Spectral calibration

Once the hyperspectral acquisition is performed, the luminance Lobj measured on each pixel for each wavelength depends on the surface reflectance properties of the object Robject, the incident irradiance Elighting, the sensor sensitivity Ssensor, the camera transmittance Tcamera and the background noise n:

obj obj lighting sensor camera L R E S T n      (6)
In the case of hyperspectral imaging, the quantities in equation ( 6) are functions of wavelength and pixel location. Hyperspectral calibration aims at retrieving the surface reflectance information from the measured luminance. Elighting, Ssensor, Tcamera, and n are unknown quantities that are calibrated through the acquisition of black and white reference Lambertian samples of respective albedo ρblack and ρwhite. The calibration results in the object reflectance for each pixel and each wavelength after applying the formula:

  The calibrated image (Fig. 6.c) seems darker than the usually perceived visual appearance of skin but the resulting color is closer to the actual skin reflectance.

This calibration method gives satisfactory results when applied on flat objects. However, for 3D objects, the incident irradiance on the reference samples is not the same as that on the object, with the difference resulting in an error in the calibration output. For nonplanar objects, the actual surface reflectance can be assessed through a more complex photometric calculation that requires the surface orientation and location, and knowledge of the incident lighting intensity in space. As the developed system provides both the spectral luminance and the 3D surface information, a "3D irradiance calibration", mentioned later, can be implemented for a consistent color measurement.

OPTICAL ANALYSIS

The optical analysis relies on two models: one for the skin structure that represents the skin tissues as a stack of two layers, and one for the light-skin interaction which predicts the reflectance of the skin according to the absorption and scattering properties of the two layers.

Skin model

Our model considers the skin as a two-layer medium, corresponding roughly to the epidermis and the dermis, whose composition and optical properties are noticeably different. The real skin structure is much more complex, as the epidermis and dermis can be further divided into several sub-layers, each of which have their own specific cell structures and chemical compositions [START_REF] Igarashi | The appearance of human skin: A survey[END_REF]. The upper layer, the epidermis, is composed of the stratum corneum, the stratum lucidum, the stratum granulosum, the stratum spinosum and the stratum basal. Keratinocytes and melanocytes, two kind of cells that produce melanin and keratin fibers, can be found in the epidermis. The second layer, the dermis, is composed of the papillary layer and the reticular layer. They contain collagen fibers, elastin fibers, fibroblasts and erythrocytes (or red blood cells) carrying hemoglobin. They are traversed by blood vessels and capillaries. Beneath the dermis, the hypodermis is an elastic layer that includes a large amount of fat cells. Moreover, the dermalepidermal junction in not flat: the dermis forms upward finger-like projections into the epidermis called dermal papillae [START_REF] Newton | Novel approaches to characterize age-related remodelling of the dermalepidermal junction in 2D, 3D and in vivo[END_REF]. The two-layer model for skin (Fig. 7) is a simplification in which each layer is flat and homogeneous. The epidermis is composed of a baseline and melanin, and the dermis is composed of a baseline, deoxyhemoglobin (Hb), oxyhemoglobin (HbO2) and bilirubin.

Skin, modeled as a two-layer turbid medium, can be characterized by the absorption and scattering coefficients of each layer. The absorption coefficients depend on the spectral absorption properties of chromophores contained in the layer and their volume fraction: according to the Beer-Lambert-Bouguer law, the spectral absorption coefficients of the epidermis Ke(λ) and dermis Kd(λ) can be written according to additive linear laws:

( ) (1 ) ( ) ( ) e m b m m K c K c K       (8)
K c c c K c K c K c K             (9)
We refer to reference [START_REF] Jacques | Optical properties of biological tissues: a review[END_REF] for the spectral absorption coefficient values of the baseline Kb(λ), of melanin Km(λ), of hemoglobin KHb(λ), of oxyhemoglobin KHbO2(λ) and of bilirubin Kbi(λ), as well as for the spectral scattering coefficients of the epidermis Se(λ) and dermis Sd(λ). The volume fractions of melanin cm, of hemoglobin cHb, of oxyhemoglobin cHbO2 and of bilirubin cbi are unknown quantities that can be estimated through an optimization based on the Kubelka-Munk theory.

The Kubelka Munk theory

The hyperspectral reflectance results are analyzed using a twoflux inverse model described in Ref. [START_REF] Seroul | Model-based skin pigment cartography by high-resolution hyperspectral imaging[END_REF], where the light propagation in the skin layers is modeled according to the Kubelka-Munk theory [START_REF] Kubelka | New contributions to the optics of intensely lightscattering materials. Part II: Nonhomogeneous layers[END_REF], whose advantage in comparison with other light scattering models is that it provides analytical formulas for the reflectance and transmittance of a multilayer scattering structure. The epidermis is characterized by its thickness e h , spectral absorption coefficient Ke(λ) and spectral scattering coefficient Se(λ). Its proper reflectance Re and transmittance Te, expressed according to the Kubelka-Munk model [START_REF] Kitware | The Visual Toolkit[END_REF], are (for each wavelength of light): The dermis is assumed to be opaque, therefore infinitely thick from an optical point of view. It is characterized by its spectral absorption coefficient Kd(λ) and spectral scattering coefficient Sd(λ), and its reflectance Rd predicted according to the Kubelka-Munk model is given by:

(a) (b) (c) ( 2 ) d d d d d d d d d K S K K S R a b S       (14)
where parameters ad and bd are functions of Kd(λ) and Sd(λ) defined similarly as in equation [START_REF] Zhang | A novel 3D multispectral vision system based on filter wheel cameras[END_REF].

The total spectral reflectance of the skin, superposition of the dermis and epidermis, is given by Kubelka's formula [START_REF] Kubelka | New contributions to the optics of intensely lightscattering materials. Part II: Nonhomogeneous layers[END_REF]:

Inverse problem and optimization

For each pixel, the volume fraction of each component is estimated through an optimization based on a spectral angle similarity metric, the Spectral Angle Mapping (SAM) [START_REF] Yuhas | Discrimination among semiarid landscape endmembers using the spectral angle mapper (sam) algorithm[END_REF] (also equivalent to a normalized correlation).

The search for chromophore volume fractions requires the calculation of a theoretical reflectance spectrum Rth(λ). The distance d defined by the SAM between Rth(λ) and the measured spectrum Rm(λ) is computed as: Thanks to the use of the SAM metric, results of the skin analysis are very good even on shadowed region in the input data. Indeed, this metric improves the algorithm robustness by reducing its sensitivity to variations due to lighting non uniformities. As an example (illustrated by Figure 9), let's compare three regions of skin. The first is a reference region characterized by the measured spectral luminance S0(λ). The second region has the same chromophore composition as the first but with a lower incident irradiance, which implies the perceived color is darker and the spectral luminance is S1(λ) = kS0(λ) with k < 1. The third region of skin differs in chromophore composition from the first two and is characterized by the measured spectral luminance S2(λ). When comparing the three spectral luminances using the rootmean-square (RMS) metric, the distance between S0(λ) and S1(λ) is high while the distance between S0(λ) and S2(λ) is low due to the invariance of this metric in respect to lightness. This comparison is consistent with perceived color, which can be affected by varying irradiance, but doesn't represent the skin chromophore properties. When comparing the three spectral luminances using the SAM metric, on the other hand, the distance between S0(λ) and S1(λ) is null, while there is a noticeable distance between S0(λ) and S2(λ). This metric is more sensitive to changes in the spectrum shape and insensitive to irradiance non-uniformities, and therefore yields better results for skin analysis optimization.

RESULTS AND PERFORMANCE

We propose in this section to present the experimental results achieved with our method, considering in turn the different steps developed above: the 3D acquisition, the hyperspectral acquisition, and the optical analysis. Lastly, 3D maps of skin parameters will be presented and discussed.

3D meshes

3D scanning based on fringe projection results in a point cloud representing the surface of the object. This point cloud has to be converted into a mesh in order to add a texture on it and display it as a 3D object. Before conversion, the depth map is filtered with a low pass Butterworth filter in order to reduce noise (Fig. 10.b). After noise reduction, the filtered depth map is converted into a .vtk 3D file using a function from The Visual Toolkit (VTK) [START_REF] Kitware | The Visual Toolkit[END_REF] library (see Figure 10.c).

The precision of the acquisition system has been defined by measuring a tilted plane at several inclinations, with and without smoothing. The maximum, average and standard deviations are presented in Table 1, before and after smoothing (Fig. 11). The smoothing step consists in reducing the periodic noise and filtering high frequencies. Table 2 shows the tilt angle between the XY plane and planes A and B, measured from the 3D data and manually. The analysis of the 3D results shows that the noise is not white but is the sum of white noise, periodic noise and deformation (due to errors in the conversion of phase map into depth map). On an opaque plane such as the one used for the measurement, the noise is mainly the result of the projector-camera non linearity. This measured deviation has to be interpreted as a minimum value characterizing the acquisition method. For in vivo skin acquisition, the noise can also be caused by: unintentional movements during acquisition; blurred fringes due to the limited depth of field of the projector [START_REF] Li | Gamma-distorted fringe image modeling and accurate gamma correction for fast phase measuring profilometry[END_REF]; diffuse fringes due to the acquired object not being opaque but turbid. The acquisition setups have to be optimized to minimize noise: the acquisition time has to be short, the projected fringes have to be blue and the camera focus has to be carefully adjusted. The acquisition precision of the current method is sufficient for full face acquisition. However, for more detailed applications, such as measurement of wrinkles, greater precision would be required.

Hyperspectral acquisition results

Images obtained with the hyperspectral camera are hypercubes of 2048 x 2048 pixels x 30 wavelengths covering the visible spectrum (410 to 700 nm) with an acquisition time of approximately 2.5 seconds. From each hypercube, it is possible to calculate the RGB colors perceived under a given illuminant after spectral calibration (Fig. 13.d). Images obtained at particular wavelengths show the specific skin characteristics: melanin stains are visible at low wavelengths (Fig. 13.a), blood vessels appear between 560 nm and 600 nm (Fig. 13.b), and finally, when seen in infrared (Fig. 13.c), the skin appears uniform, more transparent and more diffusing. The hypercube resolution is good enough for analysis, however the LCTF and the camera bring chromatic distortion, resulting in spectral errors on small objects, such as spots and hairs. Chromatic errors can also result from small movements of the subject during the acquisition process. This phenomenon is especially visible on the eyelid, which is susceptible to involuntary micro-movements caused by the lighting used during acquisition.

Optical analysis results

The skin analysis optimization yields chromophore concentrations in each pixel according to the method summarized in Section 5, then displayed under the form of maps that show the relative volume fractions of blood, oxygen and melanin (Fig. 14).

When the optical analysis is performed, an error map showing for each pixel the residual error (defined by Equation ( 16)) between the measured spectrum and the theoretical spectrum is created. The full dynamic error map (Fig. 13.a) shows that the optical model is good at distinguishing skin (low error) from other elements such as hair, fabric, make-up and background (high error). The error map on more details to be seen on the face shows that there is larger error on small elements, such as blood vessels on the eyelid. This is caused by the chromatic aberration of the acquisition system: on small elements or on borders, the pixel spectrum is modified by the chromatic aberration, impacting the pixel-by-pixel chromophore map. For now, we consider that the number of pixels impacted by this error is small enough to be negligible in the final image.

The error is also higher in areas where the incident irradiance varies, as visible on the neck, ears, lips and nostril sides. This shows that the SAM metric (part 5.3) is not sufficient to avoid the influence of varying irradiance. Another hypothesis for explaining errors in regions such as nostril sides and ears is that the geometry of the surface presents convex zones where interreflections [START_REF] Forsyth | Mutual illumination[END_REF] can occur. Interreflections consist of multiple light reflections within a cavity that modifies the perceived spectral luminance as a nonlinear function of the spectral reflectance of the concave surface. These spectral modifications impact the chromophore map reconstruction and results in an error. We intend to investigate this topic in the near future.

Finally, the optical analysis is based on the assumption that skin is Lambertian (ignoring the gloss component) and that the Kubelka Munk model can be applied in the same way, whether the observed surface is normal to the observation direction or not. However, for non Lambertian surfaces, the quantity of light reflected depends on the angle between the normal and the observation direction. An extended model taking into account the angular dependence of the skin reflectance would be necessary to further improve the accuracy of the method.

3D hyperspectral results

When the depth map is converted into a 3D mesh, texture is added: this can represent the reflectance information for a given wavelength, the color information, or the chromophore information calculated with the skin analysis algorithm (Fig. 15). Figure 15.b, c and d represent the results of the entire process, comprising 3D scanning, hyperspectral imaging and skin analysis. The visualization of the face in 3D can supplement the chromophore map information, potentially aiding the study of skin properties that implicate shape.

TOWARD IRRADIANCE CORRECTION

As mentioned in section 5.2, a simple luminance correction based on a flat reference sample is not sufficient to provide accurate spectral measurement. When a hyperspectral acquisition is performed on a full face, the reconstructed color is usually a faithful representation of the actual spectral reflectance only on the forehead

(a) (b) (c) (d) (a) (b) (c)
and the cheeks. To overcome this limitation, a calibration method taking into account the full 3D geometry of the object and lighting must be implemented. The proposed method is based on two assumptions. Firstly, the surface is assumed to be Lambertian. Secondly, light sources are approximated as punctual sources. In the proposed system, each lighting unit is built with ten LEDs, each of which is considered as a punctual light source. The total incident irradiance is the sum of each LED contribution.

Let's consider a small area of the object surface that corresponds to one pixel (x, y) when imaged on the camera sensor. For each pixel, Bouguer's law and Lambert's law result in equation [START_REF] Geng | Structured-light 3D surface imaging: a tutorial[END_REF], which gives the luminance Lx,y of the pixel as a function of: the incident intensity Ik(θx,y,k, φx,y,k) of each LED k in the direction of the pixel; the 3D surface (i.e., the distance to the LED dx,y,k and the angle αx,y,k between the surface normal and the LED axis); and the albedo ρx,y of the surface: The surface albedo, which is the parameter that has to be determined, can be retrieved by performing a lighting calibration using a black and a white flat sample to provide the missing information Ik(θx,y,k, φx,y,k) for each LED. Figure 16 shows the results obtained in a simple configuration based on a single monochromatic source, before (Fig. 16.a) and after irradiance correction (Fig. 16.b). Figure 16.c shows the irradiance horizontal profile on the forehead before and after correction. In this example, the gloss component of skin reflectance is visible as no polarizing filter has been used. The implemented method corrects the irradiance variation observed on the forehead. The curve of the blue line on Figure 16.c is the consequence of the irradiance non uniformities on the forehead but doesn't describe the skin reflectance on this zone. The red line shows that the reflectance profile is more homogeneous after applying the irradiance correction. Local variations are due to skin texture on the forehead.

CONCLUSION AND DISCUSSION

The system developed is, as far as we could see, one of the first 3D hyperspectral imaging systems conceived for living tissues, satisfying good spectral and spatial resolutions and short acquisition time. The combination of acquisition system, shadow correction and optical models is also original and brings a new dimension to skin analysis. Moreover, the proposed system contributes to the general trend of replacing an object with its numerical mockup, allowing virtual analysis and interaction with low constraints.

The use of 3D information for correcting irradiance nonuniformity and improving the measurement of spectral reflectance was the first motivation for implementing a 3D scanning system. However, it has demonstrated high potential for other applications, such as providing the shape in addition to the color, or investigating the Bidirectional Reflectance Distribution Function (BRDF) of objects without shape constraints.

Although the chromophore maps show very interesting results, the implemented method for skin analysis does not currently take into account scattering, although it is an important parameter for skin. The quantification of scattering could widen the range of applications for full face skin analysis, providing a better description of skin properties such as inflammation and aging, which are of particular interest to dermatology and cosmetology.

The acquisition system quality could also be improved by using a brighter projector to shorten acquisition time and by offering 3D scanning from more than one point of view in order to ensure good 3D reconstruction on all parts of a human face.

Figure 2 .

 2 Figure 2. Illustration of 3D scanning based on structured light projection.

Figure 1 .

 1 Figure 1. Schematic representation of the acquisition system.

Figure 3 .

 3 Figure 3. Illustration of the triangulation principle.

Figure 4 .

 4 Figure 4. 3D acquisition intermediate results: recorded fringe image (a), wrapped phase computed from 6 fringe images (b), unwrapped phase on a selected zone of interest (c).

Figure 5 .

 5 Figure 5. Scheme of the depth calibration setup.

Figure 6 .

 6 Figure 6. Acquired and corrected hyperspectral images converted into RGB color (illuminant D65): (a) Non uniform illumination on the white reference plane, (b) uncalibrated image (c) calibrated image

Figure 7 .

 7 Figure 7. Skin model used in the optical analysis of hyperspectral data.

  and

22 (

 22 consists in searching for the parameters that minimize d:

Figure 9 .

 9 Figure 9. Spectral luminances (a) and corresponding RGB colors under D65 illumination (b). Distances to the reference spectrum S0 according to the RMS metric and the SAM metric are: dRMS(S1)=14; dRMS(S2)=2; dSAM(S1)=0; dSAM(S2)=6.10 -2 .

Figure 10 .

 10 Figure 10. 3D hyperspectral acquisition results: (a) depth map before smoothing, (b) depth map after smoothing, (c) 3D mesh without texture.

Figure 11 .

 11 Figure 11. Point cloud display of the measured titled plane before smoothing (a) and after smoothing (b). The color corresponds to the deviation amplitude.

Figure 12 .

 12 Figure 12. Example of hyperspectral image: (a) Image at 420 nm, (b) image at 590 nm, (c) image at 700 nm, (d) RGB color reconstruction under a D65 illuminant.

Figure 13 .

 13 Figure 13. Error map from the skin optical analysis, (a) full dynamic (SAM ∈ [0, 0.5]) and (b) optimized dynamic (SAM ∈ [0, 0.1]) to observe errors on skin.

  Fig. 13.b, in which the display dynamic has been optimized to allow

Figure 14 .

 14 Figure 14. Optical analysis results, volume fraction map of: melanin (a), oxygen (b), blood (c).

Figure 15 .

 15 Figure 15. 3D hyperspectral acquisition results: 3D visualization of the color (a), 3D visualization of the melanin map (b), 3D visualization of the oxygen rate map (c), 3D visualization of the blood map (d).

Figure 16 .

 16 Figure 16. Original image with incident irradiance non-uniformities (a), image after correction (b), profile of the luminance intensity on the forehead for the non-corrected image (blue) and the corrected image (red) (c).

Table 1 : 3D scanning precision on the tilted plane A

 1

	Tilted plane A	Maximum	Average	Standard
	measurement	deviation	deviation	deviation
	Before smoothing	1.39 mm	0.25 mm	0.32 mm
	After smoothing	0.51 mm	0.17 mm	0.20 mm

Table 2 : Tilted plane A and B, angle measurement

 2

		3D acquisition	Manual
	Plane A	17.52°	17.35° ± 0.35°
	Plane B	33.23°	34.02° ± 0.71°

Acknowledgement

This work was supported by the Région Auvergne-Rhône-Alpes within the program ARC 6, and the French National Research Agency (ANR) within the program "Investissements d'Avenir" (ANR-11-IDEX-0007), in the framework of the LABEX MANUTECH-SISE (ANR-10-LABX-0075) of Universite de Lyon.

Author Biography

Lou Gevaux received her B.S. in optical sciences from the University of Saint-Etienne (France, 2014) and an engineering degree from the Institut d'Optique Graduate School (France, 2014). She is currently pursuing a Ph.D. at the University Jean Monnet of Saint-Etienne (France), advised by Professor Mathieu Hébert. In association with Newtone Technologies (France, Lyon) and the Department of Dermatology, Hôpital Nord (Saint-Etienne, France), her research focuses on optical methods for skin analysis.