Hamlet Medina Ruiz

Lana Iwaza

Michel Kieffer

Béatrice Pesquet-Popescu

Khaldoun Al Agha

Evaluation of Multicasting Schemes based on Joint Multiple Description and Network Coding

This paper considers a multicast scenario and compares the average reception quality obtained when combining multiple description coding (MDC) and network coding (NC). Plain (single description) network coding (NC-SDC) serves as reference.

In the considered scenario, a single source is multicast to several receivers with various channel conditions. Contrary to a NC-SDC scheme, unable to recover the coded packets when not enough combinations of packets have been received, NC of MDC packets allows a more progressive quality improvement with the number of received packets, and a reduction of the effect of the quantization noise when MDC is performed via frame expansion before quantization.

Considering a probability distribution for the bit transition probability during transmission to any user in the multicast group, the expected signal-to-noise ratio is evaluated. Performance comparisons are made for various error distributions, field sizes, and MDC methods (via frame expansion and correlating transform).

I. INTRODUCTION

The quality of reception in any multimedia transmission in cellular networks or between neighbors in ad-hoc wireless networks is highly dependent on the channel conditions between the receiver and its base station, or between neighbors.

Efficient coding schemes have to be designed in such a way that they adapt to the various channel conditions of the different users: users with better channel conditions should get a content with an improved quality, without penalizing those with worse conditions, the maximum of which should still be able to get the transmitted content with a decent quality.

Scalable or Layered Coding (LC) techniques [START_REF] Ghanbari | Two-layer coding of video signals for vbr networks[END_REF] is one way to implement such adaptive systems. In LC, transmitted data are partitioned into several quality layers with corresponding priority. The decoding quality enhances with the number of layers received. The combination of this technique with a priority encoding transmission (PET) can increase its efficiency. However, with LC, receiving a low priority layer is useless as long as all higher priority layers have not been received.

A way to overcome this drawback is to use Multiple Description Coding (MDC) [START_REF] Goyal | Multiple description coding: Compression meets the network[END_REF]. Data are partitioned into equally important descriptions, and the quality of reception gracefully improves with the number of received descriptions. The price to be paid is some additional redundancy introduced during coding. Routing of MDC data may be optimized using rainbow network coding [START_REF] Shao | Rainbow network flow with network coding[END_REF]. The reconstruction quality is maximized by optimizing the combinations of descriptions associated to each link. Nevertheless, this approach is centralized and requires a perfect knowledge of the network topology.

Network coding (NC) [START_REF] Ahlswede | Network information flow[END_REF], in multicast scenarios, may achieve the max-flow capacity between the source and each destination node [START_REF] Li | On achieving optimal throughput with network coding[END_REF]. Practical random NC schemes have been proposed, see, e.g., [START_REF] Ho | A random linear network coding approach to multicast[END_REF]. NC is quite robust to packet losses, as long as enough coded packets are received. When not enough packets are obtained, no reconstruction is possible. Joint source-network coding approaches have been proposed recently in order to get the best out of MDC and NC. Concatenated MDC and NC have been first proposed in [START_REF] Walsh | A concatenated network coding scheme for multimedia transmission[END_REF], where a network code and an outer PET code are concatenated to allow users with bad channel conditions to be able to decode the most important packets even if the reception of all packets is not guaranteed. In [START_REF] Park | Transmission of correlated information sources with network coding[END_REF], approximate decoding is performed when not enough packets are received, by exploiting existing correlations between transmitted packets. In [START_REF] Iwaza | Joint decoding of multiple-description network-coded data[END_REF], a coding scheme is proposed combining NC and MDC, where descriptions are generated either via frame expansion (NC-MDC-F), or via a correlating transform (NC-MDC-T). A good robustness against packet losses is observed in both cases. While NC-MDC-T is tolerant to more erasures, NC-MDC-F provides a better signal-to-noise ratio (SNR) when the number of lost packets is small thanks to the reduction of a part of the quantization noise. Moreover, NC-MDC-F is able in some cases to provide a reconstruction even if the number of losses is higher than the number of redundant packets introduced during frame expansion. The price to be paid is an increased computational complexity, since a mixed integer quadratic program has to be solved.

This paper evaluates the performance of both NC-MDC approaches and compares them to a scenario where only NC is performed (NC-SDC). These schemes are briefly recalled in Section II. We assume that the same number of potentially corrupted packets are available to each receiver. These packets are assumed to have passed through a binary symmetric channel (BSC) with transition probability ε distributed according to some probability density function (pdf) f (ε). This allows us to characterize for the three coding schemes the average SNR among receivers for various f (ε), as detailed in Section III. Simulation results are described in Section IV, before drawing some conclusions and introducing future research directions in Section V.

II. CODING SCHEMES

Assume that a source S has to transmit some realization x ∈ R k of a random vector to N receiver nodes. Three types of coding schemes are considered. The resulting models of the transmission chain between the source and a typical receiver are now presented. These models may correspond to various communication networks such as a (possibly relay-assisted) cellular network or a wireless ad-hoc network.

A. Scenario 1: single description network coding (NC-SDC)

In this scenario, each entry of the source vector x is quantized using a q-level uniform scalar quantizer to get a vector of quantized indexes z ∈ F k q , where F q is the Galois field with q elements. Each quantized index is put in a separate packet. Then n k random linear combinations of packets over F k q are transmitted. Further combinations may be performed within the network. A receiver then gets m n packets p ∈ F m q , see Figure 1. The effect of the network and the various network coding operations is modeled by the network coding matrix A linking z and p as follows: The estimate x of x at the receiver side requires that A is full rank k, in which case z is obtained via Gaussian elimination and x via inverse quantization of z. If A is of rank less than k, no estimate better than the mean value of x may be obtained.

p = Az. (1)

B. Scenario 2: MDC via correlating transform (NC-MDC-T)

Again, x is quantized to get quantization indexes y ∈ F k q . Then, a full-rank k correlating transform T ∈ F n×k q is applied to y to get a vector z = T y ∈ F n q . Entries of z are again put in separate packets which are network coded, and p ∈ F m q represents the received packets, see (1) and Figure 2. Since T ∈ F n×k q is of full rank k, there exists a parity-check matrix H ∈ F (n-k)×n q of rank n -k such that Hz = 0 for all z such that there exists y ∈ F k q with z = T y. This property may be used at receiver side to estimate x from p. Provided that the matrix

A = A H (2)
is of rank n, z may again be obtained via Gaussian elimination, and x via inverse quantization, see [START_REF] Park | Transmission of correlated information sources with network coding[END_REF], [START_REF] Iwaza | Joint decoding of multiple-description network-coded data[END_REF]. When A is of rank less than n, again, no better estimate than the mean value of x may be obtained.

C. Scenario 3: MDC via frame expansion (NC-MDC-F)

Here, the redundancy is introduced via a frame expansion matrix F ∈ R n×k , producing the vector y = F x ∈ R n with n k. Each entry of y is then quantized with a q-level uniform scalar quantizer to get a vector of indexes z ∈ F n q . The entries of z are then processed as in Scenarios 1 and 2. At the receiver side, p and the fact that x has been expanded may be exploited to obtain x. In [START_REF] Iwaza | Joint decoding of multiple-description network-coded data[END_REF], this problem is cast in the framework of mixed integer quadratic programming. The problem may then be modeled, e.g., with AMPL [START_REF] Fourer | The AMPL Book[END_REF] and solved, e.g., using CPLEX [START_REF]ILOG CPLEX 11.0 User's Manual[END_REF].

Remark 1: In the previous scenarios, packets contain only a single quantized sample. In what follows, we assume that the source waits for realizations of the source vector before starting transmission. Packets of samples are then transmitted, which improves the payload length over total packet length ratio.

III. EXPECTED SIGNAL-TO-NOISE RATIO FOR A GIVEN

RECEIVER

When considering a single receiver, with known channel characteristics, there is no need to perform any joint coding scheme, a separate optimization of the data compression and transmission being optimal, see [START_REF] Shannon | A mathematical thoery of communication[END_REF].

In the considered multicast scenario, the answer is not so clear. The aim of this paper is thus to evaluate the performance of the three previously-introduced coding schemes as a function of the size q of F q and for different channel models.

A. Hypotheses

The number of network-coded packets (size of the generation) is g. For NC-SDC one has g = k whereas for NC-MDC-T and NC-MDC-F, g = n. It is assumed that each of the N receivers obtains n noisy packets from which an estimate of x has to be evaluated. The packets are assumed to have passed through a binary symmetric channel (BSC) with transition probability ε distributed according to some probability density function (pdf) f (ε), identical for all the receivers. This BSC model allows us to account for the packet length, which is not the same in the three considered scenarios.

The fact that all receivers have obtained the same number of noisy packets is quite realistic in a cellular network. It is more questionable in an ad-hoc network, where the min-cut is likely not to be the same between the source and all users.

B. Average signal-to-noise ratio

The average SNR observed by a given receiver depends on the number of received packets as well as the size of the Galois field q. For a given user, whatever the considered scenario, one has to determine the probability mass function (pmf) P Γ (γ | g, n, q) of the number γ independent packets received, i.e., the pmf of the rank of A, given that g packets were combined, n noisy combinations were received, and that the size of the Galois field is q.

If i is the number of error-free packets received, one gets

P Γ (γ | g, n, q) = n i=γ P (γ, i | g, n, q) = n i=γ P (γ | i, g, n, q) P (i | g, n, q) = n i=γ P R (γ | i, g, q) P C (i | g, n, q) .( 3
)
In (3), P C (i | g, n, q) is the probability of getting i noisefree (Clean) packets (this probability depends on f (ε)) and P R (γ | i, g, q) is the probability of having γ informative packets among the i received packets without errors (this is the probability that the rank of A is γ). To evaluate P C (i | g, n, q), one has to introduce f (ε) as follows

P C (i | g, n, q) = P (i, ε | g, n, q) dε = P (i | g, n, q, ε) f (ε | g, n, q) dε = P (i | g, n, q, ε) f (ε) dε, (4)
since f (ε) does not depend neither on n, nor on q. Then

P (i | g, n, q, ε) = n i (1 -ε) iL 1 -(1 -ε) L n-i (5)
indicates the probability of receiving i packets without errors among n received packets sent for a realization ε of the BSC transition probability. In ( 5), L = (g + + seq + crc ) log 2 q represents the total length in bits of a packet. seq and crc stand for the number of symbols in F q used to represent the sequence number (SEQ) and the CRC, respectively, and is the payload. SEQ and CRC form part of the header of networkcoded packets [START_REF] Chou | Practical network coding[END_REF]. The CRC is used to protect the data, the coefficients involved in the linear combination, and SEQ. Now, P R (γ | i, g, q) represents the pmf of the number γ independent linear combinations among the i received, and which is given by

P R (γ | i, g, q) = µ (i, g, γ, q) q ig (6)
where µ (n 1 , k 1 , r 1 , q) [START_REF] Van Lint | A Course in Combinatorics[END_REF] is the number of matrices in

F n1×k1 q with rank equal to r 1 µ (n 1 , k 1 , r 1 , q) = k 1 r 1 q r1 i=0 (-1) (r1-1) r 1 i q q n1i+( r 1 -i 2 )
(7) and n1 k1 q is the Gaussian coefficient [START_REF] Van Lint | A Course in Combinatorics[END_REF] given by

n 1 k 1 q =    1 k 1 = 0 (q n 1 -1)(q n 1 -1 -1)...(q n 1 -k 1 +1 -1) (q k 1 -1)(q k 1 -1 -1)...(q-1) k 1 > 0 ( 8
)
Assume now that the average SNR for Scenario s = 1, . . . , 3 when receiving γ independent coded packets in F q is given by SNR s (γ, q), the average SNR accounting for the transition probability distribution f (ε) is deduced from ( 4) and ( 6) and is expressed as

SNR s (q) = n γ=0 SNR s (q, γ)P Γ (γ | g, n, q) . ( 9
)
When considering NC-SDC or NC-MDC-T, a closed-form expression for SNR s (γ, q) may be obtained. In the case of NC-MDC-F, such expression is more complex to obtain. The following section characterizes the average SNR for various f (ε) and sizes of the Galois field using [START_REF] Iwaza | Joint decoding of multiple-description network-coded data[END_REF].

IV. SIMULATION RESULTS

Simulations are performed with k = 6 and n = 9. The source generates vectors of k independent and identically distributed Gaussian samples with zero-mean and variance σ 2 = 1. Quantization with q ∈ {7, 17, 31, 61} outputs are considered, leading to quantization indexes in F q . The network matrix A is chosen at random in the corresponding Galois field to simulate the effect of network coding. Simulation results are averaged over 1000 realizations of the source and of the network matrix A.

A. Performance as a function of the number of lost packets

First, the average SNR resulting from the reconstruction of the transmitted message is drawn as a function of the number of lost packets for various sizes of the Galois field used in the quantization and NC operations. These curves correspond to the representation of γ SNR s (q, γ)P R (γ | i, g, q) with the number of lost packets being n -i. The potential rank deficiency is thus taken into account, but not the effect of transmission errors.

Figures 4 and5 represent the SNR for NC-MDC-F, NC-MDC-T, and NC-SDC as a function of the field size. For NC-SDC, n network-coded packets containing k independent packets are transmitted (the redundancy is almost the same as in the two first scenarios, only the size of the header is somewhat reduced due to the reduced number of independent packets). As stated in the introduction, Figures 4 and 5 Fig. 4. SNR as a function of the number of lost packets and q, for q = 7 and q = 17. NC-MDC-F (diamond), NC-MDC-T (square), and NC-SDC (circle)

show that NC-MDC-F can mitigate part of the quantization noise, providing better results than the NC-MDC-T when the BSC transition probability is small. NC-MDC-T provides an increased robustness to packet losses compared to NC-MDC-F. Fig. 5. SNR as a function of the number of lost packets and q, for q = 31 and q = 61. NC-MDC-F (diamond), NC-MDC-T (square), and NC-SDC (circle) Figures 4 and 5 also show a degradation of the SNR due to the rank deficiency, even if enough packets have been received. When not enough packets are received, the SNR is null.

B. Average SNR as a function of ε

In this case the average SNR is evaluated as a function of the channel transition probability ε, with q = 31. Figure 6 represents the SNR as a function of ε when ε is between 10 -6 and 10 -1 for packets with payloads of length = 100 symbols. NC-MDC-F outperforms NC-SDC for values of ε in [0, 1.5 × 10 -4 ], but for larger values of ε the average SNR with NC-SDC is larger.

C. Average SNR for various distributions of the transition probability

In this case we study the received signal quality for two different pdfs of the transition probability ε. The uniform distribution f U (ε) = 1 a 0 ≤ ε ≤ a 0 elsewhere [START_REF] Fourer | The AMPL Book[END_REF] and the exponential distribution

f E (ε) = ln(10)10 -ε 10 -a -10 -b a ≤ ε ≤ b 0 elsewhere (11)
are considered. Using [START_REF] Fourer | The AMPL Book[END_REF], (4) becomes

P C (i | g, n, q) = 1 a n i a 0 (1 -ε) iL 1 -(1 -ε) L n-i dε.
(12) Using ( 11), (4) becomes

P C (i | g, n, q) = ln(10) 10 -a -10 -b n i b a 10 -ε (1 -ε) iL 1 -(1 -ε) L n-i
dε.

(13) The average SNR is computed using (9) after substituting [START_REF] Shannon | A mathematical thoery of communication[END_REF] or ( 13) in [START_REF] Shao | Rainbow network flow with network coding[END_REF].

With NC-MDC-F and NC-MDC-T, the number of mixed packets is g = n and with NC, it is g = k. Figure 7 represents the expected SNR for the presented techniques for various field sizes q = {7, 17, 31, 61} using the uniform pdf for the transition probability as a function of a, the upper bound of the support of the uniform pdf.

Similarly, Figure 8 represent the expected SNR for the exponential pdf for the transition probability with a = 10 -10 as a function of b, from 10 -6 to 10 -1 .

In both cases, when there are not too many users with bad channel characteristics, NC-MDC-F outperforms NC-MDC-T, which itself outperforms NC. When the proportion of users with bad channel gets larger (a > 3.10 -4 for the uniform pdf and b > 3.10 -4 for the exponential pdf), NC-SDC becomes better.

V. CONCLUSIONS AND PERSPECTIVES

This paper considers the effects of the field size and channel conditions in the average signal quality received, in the scenario where a single source multicasts data to a set of users experiencing various channel conditions, using NC or a combination of MDC and NC.

Packets are assumed to be received at the output of a BSC which transition probability ε is taken as random according to some pdf f (ε). The performance of the three schemes largely depends on f (ε). When a user is likely to have good channel conditions, NC-MDC-F provides the best results. When his channel conditions may be quite bad, then plain NC-SDC becomes better.

In future work, a distribution of the number of noisy packets obtained by the various receivers may be introduced to have a more realistic description of packet reception in wireless ad-hoc networks. One may also consider the distribution of channel SNR as a function of the user location when considering a cellular network. The effect of relays has also to be better taken into account.

Fig. 1 .

 1 Fig. 1. Block diagram of the NC-SDC scenario

Fig. 2 .

 2 Fig. 2. Block diagram of the MDC scheme via correlating transform (NC-MDC-T)

Fig. 3 .

 3 Fig. 3. Block diagram of the MDC scheme via frame expansion (NC-MDC-F)

  -T, q=31 NC-MDC-F, q=31 NC, q=31 NC-MDC-T, q=61 NC-MDC-F, q=61 NC, q=61

Fig. 6 .

 6 Fig. 6. Average SNR as a function of the transition probability ε

Fig. 7 .Fig. 8 .

 78 Fig. 7. Average SNR for the Uniform Distribution as a function of the upper bound of the support of the pdf of transition probability f (ε)

ACKNOWLEDGMENTS

This work has been partly supported by DIM-LSC NC2, DIM-LSC SWAN. Michel Kieffer was partly supported by the Institut Universitaire de France.


Thierry Gautier

Fabien Lementec

Vincent Faucher

X-Kaapi Fortran programming interface

Keywords: parallel computing, X-Kaapi, Fortran environment de programmation parallèle, X-Kaapi, Fortran

This report defines the X-Kaapi Fortran programming interface.

1 Software installation X-Kaapi is both a programming model and a runtime for high performance parallelism targeting multicore and distributed architectures. It relies on the work stealing paradigm. X-Kaapi was developed in the MOAIS IN-RIA project by Thierry Gautier, Fabien Le Mentec, Vincent Danjean and Christophe Laferrière in the early stage of the library.

In this report, only the programming model based on the C API is presented. The runtime library comes also with a full set of complementary programming interfaces: C, C++ and STL-like interfaces. The C++ and STL interfaces, at a higher level than the C interface, may be directly used for developing parallel programs or libraries.

Supported Platforms

X-Kaapi targets essentially SMP and NUMA platforms. The runtime should run on any system providing:

• a GNU toolchain (4.3),

• the pthread library,

• Unix based environment.

It has been extensively tested on the following operating systems:

• GNU-Linux with x86_64 architectures,

• MacOSX/Intel processor.

There is no version for Windows yet.

X-Kaapi Contacts

If you wish to contact the XKaapi team, please visite the web site at: http://kaapi.gforge.inria.fr

RT n°429

2 Initialization and termination

Description

KAAPIF_INIT initializes the runtime. It must be called once per program before using any of the other routines. If successful, there must be a corresponding KAAPIF_FINALIZE at the end of the program.

Parameters

• FLAGS : if not zero, start only the main thread to avoid disturbing the execution until tasks are actually scheduled. The other threads are suspended waiting for a parallel region to be entered (refer to KAAPIF_BEGIN_PARALLEL).

Return value

Refer to the Error codes section.

Example

Description

KAAPIF_GET_VERSION retrieves a string describing the runtime version. Currently, the string is the git last commit hash in base64.

Parameters

• VERSION : An array large enough to hold the version.

Return value

None.

Example

Description

Concurrency related routines.

Return value

KAAPIF_GET_CONCURRENCY returns the number of parallel thread available to the X-Kaapi runtime.

KAAPIF_GET_THREAD_NUM returns the current thread identifier. Note it should only be called in the context of a X-Kaapi thread.

For both functions, a negative value means an error occured. Refer to the Error codes section.

Example

Description

KAAPIF_SET_GRAINS sets the adaptive loop grains. Grains are used to amortize the act of extracting work for both the parallel (ie. during a steal) and sequential (ie. during a pop) executions. Guessing those grains is problem specific. The general intuition is that they should be inversely propotionnal to a single task processing time (given a task processing time is constant across execution). The parallel grain should be greater than the sequential grain, since a steal operation requires a bigger task to be amortized.

Note that this routine sets global variables used by a subsequent calls to KAAPIF_FOREACH family functions.

To avoid reentrancy issues, KAAPIF_SET_GRAINS should be called just before the corresponding KAAPIF_FOREACH function.

Parameters

• PAR_GRAIN : below this size, a task cannot be split for subsequent parallel execution. default to 32.

• SEQ_GRAIN : the size used by the sequential execution to extract work from its local workqueue. default to 16.

Return value

Refer to the Error codes section.

Example

Description

Capture the current time. Used to measure the time spent in a code region.

Parameters

None.

Return value

The current time, in microseconds.

Example

Refer to examples/kaapif/foreach

Description

Those routines run a parallel loop over the range [FIRST, LAST ] (note this is an inclusive interval). The loop body is defined by BODY whose arguments are given in parameters. It must have the following prototype: SUBROUTINE BODY( I , J , TID , . . . )

• [I, J] the subrange to process (note that interval is inclusive)

• TID the thread identifier

Parameters

• FIRST, LAST : the iteration range indices, inclusive.

• NARGS : the argument count • BODY : the function body to be called at each iteration • ...: the arguments passed to BODY. For KAAPIF_FOREACH_WITH_FORMAT, refer to the KAAPIF_SPAWN documentation.

Return value

Refer to the Error codes section.

RT n°429

Example

Description

Create a new computation task implemented by the function BODY.

BODY is called with the user specified arguments, there is no argument added by XKAAPI:

SUBROUTINE BODY(ARG0, ARG1, . . . )
Each task parameter is described by 4 successive arguments including:

• the argument VALUE,

• the parameter TYPE,

• the element COUNT,

• the access MODE.

TYPE is one of the following:

• KAAPIF_TYPE_CHAR=0, • KAAPIF_TYPE_INT=1, • KAAPIF_TYPE_REAL=2, • KAAPIF_TYPE_DOUBLE=3.
If a parameter is an array, COUNT must be set to the array size. For a scalar value, it must be set to 1.

RT n°429

MODE is one of the following:

• KAAPIF_MODE_R=0 for a read access,

• KAAPIF_MODE_W=1 for a write access,

• KAAPIF_MODE_RW=2 for a read write access,

• KAAPIF_MODE_V=3 for a parameter passed by value.

Parameters

• NARGS : the argument count.

• BODY : the task body.

• ...: the VALUE, TYPE, COUNT, MODE tuple list.

Return value

Refer to the Error codes section.

Example

Description

KAAPIF_BEGIN_PARALLEL and KAAPIF_END_PARALLEL mark the start and the end of a parallel region. Regions are used to wakeup and suspend the X-Kaapi system threads so they avoid disturbing the application when idle. This is important if another parallel library is being used. Wether threads are suspendable or not is controlled according by the KAAPIF_INIT parameter.

Parameters

• FLAGS : if zero, an implicit synchronization is inserted before leaving the region.

Return value

Refer to the Error codes section.

Example

Refer to examples/kaapif/dfg

Description

Synchronize the sequential with the parallel execution flow. When this routine returns, every computation task has been executed and memory is consistent for the processor executing the sequential flow.

Return value

Refer to the Error codes section.

Example

  -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr


T Le-Duc

G Meunier

O Chadebec

J.-M Guichon

A New Integral Formulation for Eddy Current Computation in Thin Conductive Shells

Keywords: Integral equation method, quasi-static problems, surface impedance condition, thin conductive shells

In order to compute eddy current distributions in thin conductive nonmagnetic shells, a new integral formulation is proposed. The method is based on a surface impedance condition which takes into account the field variation through depth due to skin effect. It is general and enables the modeling of various problems whatever their skin-depth and avoiding the meshing of the air region.

I. INTRODUCTION

I NTEGRAL formulations dedicated to the computation of eddy currents in thin shells for the case of a skin depth much greater than the thickness has already been proposed by many authors [START_REF] Salon | An intergro-differential approach to eddy currents in thin plates[END_REF]- [START_REF] Koch | Modeling thin conductive sheets using shell elements in magnetoquasistatic field simulations[END_REF]. In this case, the eddy current distribution is supposed to be uniform across the thickness and surface elements are used avoiding the meshing of the shell's thickness. Such formulations are currently well-known and have shown good accuracy with a few numbers of elements in comparison with finite element methods where the air region needs to be meshed.

However, the computation of eddy currents in thin shells in the general case ( or for instance) is still a difficulty and has only been studied by few authors.

In [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], a general shell element formulation has been proposed. Based on a pseudo-analytical solution (solution of the 1D problem) and a nodal approximation of the scalar magnetic potential, the field variation across the thickness of the shell (i.e., the skin depth) is taken into account with a quite good accuracy. In this formulation, the air region is modeled by boundary integral equations and the coupling is realized at the interface on both sides of the shell.

In [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF], a very similar shell element for modeling thin conducting regions has been proposed. Like in [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], this element takes into account the field variation through depth due to skin effect but in this paper the outside regions are modeled with finite elements method. However, the approach presented in [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF] leads to the mesh of the air region. If shielding devices are modeled, this task is complex and can lead to inaccuracies.

This paper presents a new integral formulation which allows the modeling of nonmagnetic conductive thin shells in the case general ( or or ). Based on a simple surface discretization, the number of unknowns considerably reduces, the depth of the shell being not meshed. Moreover, like in [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], it does not require the meshing of the air region. An advantage of our new formulation in comparison with existing integral ones is its compactness and its generality. Indeed, it can be applied easily to every shell geometry, either delimited by closed surface or not. Moreover, we will see that the formulation is less consuming in terms of degrees of freedom than [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF].

In the second part of the paper, a surface impedance condition for the shell will be presented. In the third part, this surface condition will be coupled with an integro-differential approach. Finally, two numerical examples will be proposed in the last part. Results obtained with our formulation will be compared by those given by the Finite Element Method (FEM).

II. SURFACE IMPEDANCE CONDITION FOR THE SHELL

We consider a nonmagnetic shell with a thickness with a skin depth (Fig. 1). Both thin shell and surrounding regions are nonmagnetic , the conductivity is . In our approach, the tangential component variation of the magnetic field across the thickness of the shell is approximated by the analytical solution of the problem obtained for an infinite plane. For such problem, we get [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]:

(1)
where, and and are the tangential field values on both sides of the shell. The volume eddy curent density being tangential, thus by applying Ampere's law, we get:

(2) By introducing (1) in (2), we get:

(3) Expression of electric field on side "1" of the shell is:

(4)
Equations ( 4) and ( 3) are combined to get:

(5) 0018-9464/$31.00 © 2012 IEEE with In a second step, Faraday's law is projected normally to the shell: [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF] Applying Galerkin method to (6) on side 1 of the shell surface leads to:

j! 0 w 1 B 1 n d0 = 0 0 w(curlE) 1 n d0 (7)
where is a set of nodal surface weighting functions, the shell being limited to . Applying property in the second term of ( 7), we get following expression:

j! 0 w 1 B 1 n d0 = 0 (grad s w 2 E) 1 n d0 0 0 curl(w 1 E) 1 nd0: (8)
We switch terms in the mixed product of the second integral of ( 8) and apply Stokes theorem to the third one to get:

j! 0 w 1 B1n d0= 0 (E 2 n) 1 grad s w 1 d00 w 1 E 1 d (9)
where is the line region delimiting the surface. Using ( 5) and ( 9), we finally have:

0 grad s w 1 (H 1s 0 H 2s )d0 +j! 0 w 1 B 1 1 n 1 1 d0 + w 1 E 1 d = 0: (10)
On the side 2 of the shell, we get a similar equation:

(11) Equations ( 10) and (11) represent the electromagnetic behavior of the conductive shell itself. In order to take into account the surrounding region, it remains to couple the two formulations. In [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF], the authors choose to use the Finite Element Method (FEM), thus they need to mesh the air region. In [START_REF] Krähenbühl | Thin layers in electrical engineering. Example of shell models in analyzing eddy-currents by boundary and finite element methods[END_REF], a Boundary Equation Method (BEM) is preferred but the authors have to manage the coupling on both sides of the shell and have to determine if both sides of the shell are interfaced with a single region or with two different ones. This leads to the development of quite sophisticated geometrical analysis tools if the method has to be applied to general geometry. Moreover, four unknowns per node are needed. We have preferred a new approach, based on the coupling of ( 10) and ( 11) with an integral volume equation. The obtained formulation presents an interesting compactness treating different air volumes as the single physical region. Moreover two unknowns per node are needed (two potentials on both sides of the shell).

III. VOLUME INTEGRAL FORMULATION

Let us now assume that the shell is placed in an inductor field . On the side "1" of the shell, the total magnetic field is the sum of and , the reaction of the eddy current in thin shell. Thanks to the law of Biot and Savart, we can write:

(12)
where is the vector linking the integration point to the point where the field is expressed (Fig. 2). Additional expressions can be given: A reduced magnetic scalar potential is now introduced for each side of the shell. Since there is a jump of the tangential component of the magnetic field through the element, a double layer node surface element is needed (Fig. 3) [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF]. Let us notice that if two potential values are associated to each node, only surface elements are needed in order to discretized the shell. Expressions of tangential magnetic fields on both sides and outside of the shell are:

(18) (19) where is the inductor field limited to its tangential component (same value on both sides of the shell, the thickness being small) and stand for magnetic scalar potentials on both sides. Using (17), ( 18) and ( 19), we obtain the final formulation corresponding to the side "1" of the shell:

(20)

The equation on the side "2" of shell is obtained by permuting the subscripts "1" and "2". Both equations are solved simultaneously in order to get both potential and distributions. It should be pointed out that on line, delimiting the external border of the shell, the potential jump is set to zeros (continuity of the potential ensuring that the current is tangential to the edge) so both integral terms vanish so are not taken into account in the final resolution.

These equations have to be discretized. The easiest way is to mesh the average surface , situated halfway between boundaries of the shells (Fig. 1), into triangular elements associated with a uniform tangential component of the eddy current (meaning 1-order shape functions for the potential). The obtained algebraic linear system has size (two complex magnetic scalar potentials per node; the mesh is being composed of nodes). Let us notice that shape functions with higher orders can be used.

It must be pointed out that this formulation looks like a mix between finite element method and integral volume one. In (20), both first terms leads to sparse matrices (similar to which given by a FEM standard integration process) but the last one in- and is associated to an integral method. The matrix system obtained is so fully populated and thus the integration time and the memory needed could quickly become prohibitive. Compression algorithms must be used if problems with large numbers of elements are addressed. The use of the well-known «Fast Multipole Method» (FMM) [START_REF] Buchau | Parallelization of a fast multipole boundary element method with cluster openMP[END_REF], [START_REF] Banucu | Design and optimization of a contactless magnetically levitated actuator for 4-axis-machining[END_REF] to compress law of Biot and Savart is well-known to be very efficient in such configuration. Moreover, integrals (15) and ( 16) can be singular if we consider the influence of the element on itself. The integral is computed thank the use of an adaptive quadrature.

IV. NUMERICAL EXAMPLES

The new formulation has been implemented for a 3-D geometry and has been validated with two different academic examples. A comparison with FEM method using FLUX software [START_REF] Buchau | Parallelization of a fast multipole boundary element method with cluster openMP[END_REF] is made.

A. Hollow Sphere

Our formulation is based on the infinite plane assumption for the analytical computation of skin effect across the shell thickness. Even if this assumption has already demonstrated its efficiency, some inaccuracies can appear especially closed to shells edges or corners. Our first validation case deals with geometry without this kind of configuration: A conductive hollow sphere ( m, mm, S/m) placed in a uniform axial magnetic field (A/m) (Fig. 4) is considered. Three configurations have been tested: -case 1:

Hz, skin depth greater than the thickness ( mm, ) -case 2:

Hz, skin depth is comparable to the thickness ( mm, ) -case 3:

Hz, skin depth smaller than the thickness ( mm, ) For each case, eddy current losses (Joule losses) in the sphere have been calculated and compared with results given by axisymmetric FEM (Table I). The shell has been meshed with volume elements. In a second step, we have computed the magnetic field in the air region close to the device (calculated on the path AB in Fig. 3) with both methods for two cases: Hz and Hz (Fig. 5). We can see small differences between computed values with our formulation and the reference (axisymmetric FEM). These differences can be explained by the infinite plane assumption and certainly by the use of a 1-order triangular coarse mesh. A good solution would be to refine it but matrix compression algorithms have not been implemented yet.

B. Thin Conductive Disk

The second test case is a thin conductive disk ( m, m, S/m) placed in a magnetic field created by a loop fed with A (Fig. 6). This example is modeled by three different methods. The first one is the axisymmetric FEM. The second one is a shell element formulation implemented in 3-D FEM code [START_REF] Guérin | A shell element for computing 3D eddy current-Application to transformers[END_REF], [10]. The last one is the considered integral method. We still focus on the computed eddy current distribution and Joule losses in the disk at different frequencies (see Table II, Figs. 7 and8).

If we consider that axisymmetric FEM method as our reference, the integral method leads to an error of 2.7% for the Joule losses at Hz (Fig. 8). The lower the skin depth is, the lower the error is. These errors can be explained by the modeling of the geometry with a strong edge effect where the infinite plane assumption leads to inaccuracies. Let us notice that our integral formulation seems to lead to more accurate results than the same implementation of the shell element but with the air region treated with a 3-D FEM. Maybe a more precise nu- merical study should be made in order to better understand this difference.

V. CONCLUSION

In this paper, we have presented an original integral formulation using shell elements in order to model thin conductive nonmagnetic regions. The formulation is general and various skin effects across thickness ( or or are taken into account. The coupling of this formulation with compression algorithm seems to be promising in order to model real shielding applications. In further work, similar formulations enabling the modeling of thin conductive magnetic shells will be focused.

Fig. 1 .

 1 Fig. 1. Thin conductive and nonmagnetic region and associated notations.

Fig. 2 .

 2 Fig. 2. Magnetic field created by eddy current distribution and associated notations.

Fig. 3 .

 3 Fig. 3. Shell element representation with associated potential jump.

Fig. 4 .

 4 Fig. 4. Hollow sphere, notations and path where the external magnetic field is computed.

  volving

Fig. 5 .

 5 Fig. 5. Magnetic field in the path AB calculated by 2 methods with f = 50 Hz and f = 4000 Hz.

Fig. 6 .

 6 Fig. 6. Thin conductive disk, exciting coil and notations. Disk characteristics: R = 1 m, e = 20E 0 3 m, = 6E7 S/m. Coils characteristics: R = 0:5 m, located 1 m high above the disk, I = 1 A.

Fig. 7 .

 7 Fig. 7. Surface distribution of current (A/m) in the disk.

Fig. 8 .

 8 Fig. 8. Difference in % with the axisymmetric FEM solution for the Joule losses computation versus e=.

TABLE I JOULE

 I LOSSES OF THE HOLLOW SPHERE

TABLE II JOULE

 II LOSSES OF THE DISK (W)


Romain Vergne
email: romain.vergne@gmail.com

Pascal Barla

Roland Fleming

Xavier Granier

Surface Flows for Image-based Shading Design

Keywords: CR Categories: I.3.3 [Computing Methodologies]: Computer Graphics-Picture/Image Generation Shading, differential analysis, image deformation

come   L'archive ouverte pluridisciplinaire

Surface Flows for Image-based Shading Design

Romain Vergne 1 * Pascal Barla 2 Roland W. Fleming 1 Xavier Granier 2 1 Justus-Liebig-Universität Gießen 1 Introduction

It has long been known that texture gradients and shading patterns are critically important monocular cues to 3D shape (e.g., [START_REF] Gibson | The Perception of the Visual World[END_REF][START_REF] Koenderink | Photometric invariants related to solid shape[END_REF][START_REF] Todd | Perception of Surface Curvature and Direction of Illuminant from Patterns of Shading[END_REF][START_REF] Todd | Perception of threedimensional form from patterns of optical texture[END_REF]). However, how the brain extracts shape from these cues is still poorly understood. Recently, some authors have suggested that the visual system extracts certain shape and material properties directly from the spatial patterns of intensity gradients in the image (e.g., [START_REF] Breton | Shadows and shading flow fields[END_REF][START_REF] Li | Perception of three-dimensional shape from texture is based on patterns of oriented energy[END_REF][START_REF] Huggins | Folds and cuts: How shading flows into edges[END_REF]). In particular, [START_REF] Fleming | Specular reflections and the perception of shape[END_REF] found that for typical surfaces, the local shading orientation signals are organized into smoothly-varying distortion patterns, which are systematically related to 3D surface shape. Since then, computer vision researchers have shown such flows can be used to identify specular surfaces and reconstruct their shape [START_REF] Delpozo | Detecting specular surfaces on natural images[END_REF][START_REF] Vasilyev | Shape from specular flow: Is one flow enough?[END_REF][START_REF] Tappen | Recovering shape from a single image of a mirrored surface from curvature constraints[END_REF]]. More recently, Fleming et al. [2011a] found that similar measurements also play a key role in perceiving shape from texture patterns. However, for textures, the image distortion measurements indicate surface orientation (i.e., first-order properties) rather than curvatures (i.e., second-order properties) as they do for shading and specular reflections. A number of researchers have also suggested that different distortion patterns-and the intensity gradients within them-strongly affect not only perceived surface shape, but also perceived surface material (e.g., [START_REF] Motoyoshi | Luminance re-mapping for the control of apparent material[END_REF][START_REF] Anderson | Image statistics do not explain the perception of gloss and lightness[END_REF]Kim and Anderson 2010;Fleming et al. 2011b]).

The main objective of this paper is to take these observations one step further: by deforming images in 2D to produce appropriate distortion patterns, it should be possible to reliably create the illusion of 3D shapes and materials, without simulating the physics of light transport. In spirit, this approach is similar to Image Based Material Editing [START_REF] Khan | Image-based material editing[END_REF]], except that instead of altering the appearance of objects in existing images, we aim to allow the user also to create novel images. Like [START_REF] Ramanarayanan | Visual equivalence: towards a new standard for image fidelity[END_REF] we exploit the fact that the human visual system is highly sensitive to certain types of distortion-those that indicate object properties-while being almost completely oblivious to large deviations from physical accuracy in other respects. This approach, which we call image-based shading design, raises two intricately linked issues. First, we must identify what kind of deformation operators lead to perceptually plausible shape and material depiction.

For this purpose we assume that normal and depth images are given as input and relied upon to produce appropriate deformations. This is a reasonable assumption since a variety of solutions exist for creating normal and depth maps from 3D scenes, color images and even from scratch (e.g., [START_REF] Johnston | Lumo: illumination for cel animation[END_REF][START_REF] Toler-Franklin | Illustration of complex real-world objects using images with normals[END_REF][START_REF] Wu | ShapePalettes: interactive normal transfer via sketching[END_REF]Quixel 2011]). Second, we must give artists control over the form and magnitude of the deformation, and the content of the images being deformed, to alter the resulting shading patterns.

We address these issues by making the following contributions:

• Using a differential analysis of reflected radiance in the picture, we identify the two different types of surface flows at the origin of the distortion patterns observed in shaded images.

Related work

Drawing is the technique that gives artists the most immediate control for creating and manipulating appearance directly in pictorial space. However, the creation of complex (e.g., photoreal) pictures demands exceptional drawing skills and a huge amount of time. To address this issue, a number of computer-assisted drawing techniques have been introduced. For instance, when working with vector graphics, the production of smooth-shaded images is greatly helped by tools such as Gradient Meshes (e.g., [START_REF] Sun | Image vectorization using optimized gradient meshes[END_REF]) or Diffusion Curves [START_REF] Orzan | Diffusion curves: a vector representation for smooth-shaded images[END_REF], that provide options for using photographs as references. Although these have been primarily used for the creation of color gradients, they might be extended to the production of texture patterns as exemplified in the work of [START_REF] Jeschke | Estimating color and texture parameters for vector graphics[END_REF]. Similar approaches have been proposed for raster graphics drawing applications, such as Gradient-domain Painting [START_REF] Mccann | Real-time gradientdomain painting[END_REF]. Although these methods considerably speed up the image creation process, they still require advanced skills and become time-consuming with increasing image complexity. This may be made easier by using a photograph as a reference for shading and textures, but then the choice of materials and illumination is tied to the ones found in the input photograph.

Rendering is the natural alternative when one wants to freely vary material and illumination parameters. However, controlling appearance directly in the picture plane requires specific solutions. Inverse rendering techniques [START_REF] Okabe | Illumination Brush: Interactive Design of All-Frequency Lighting[END_REF][START_REF] Pellacini | Lighting with paint[END_REF][START_REF] Pellacini | envylight: an interface for editing natural illumination[END_REF] rely on known material characteristics and user scribbles to infer plausible illumination environments through a constrained optimization process. Conversely, shaded appearance may be designed through a painting interface by inferring material characteristics under known illumination conditions, using optimization as well. It has been demonstrated for both realistic [START_REF] Colbert | BRDF-Figure 13: Material design. Top: logo styles are created in a matter of seconds with the brush tool. Middle: photorealistic effects obtained by deforming the background image. Bottom: fabrics appearance produced with carefully positioned highlights. Shop: Creating Physically Correct Bidirectional Reflectance Distribution Functions[END_REF]] and stylized shading [START_REF] Anjyo | Tweakable light and shade for cartoon animation[END_REF][START_REF] Pacanowski | Sketch and Paint-based Interface for Highlight Modeling[END_REF]. Although they considerably reduce the amount of trial and error required to reach a satisfying result, both kinds of approaches still offer a limited control on final appearance and are restricted either to physically-realistic renderings or specific stylized appearances. [START_REF] Ritschel | Interactive reflection editing[END_REF]2010] provide solutions to depart from the laws of physics by deforming shading effects such as reflections, shadows or textures with user-controlled handles. However, their methods are still bound to the specification of a 3D scene (objects, lights and materials). Therefore, shading effects cannot be created as with drawing applications, but only modified.

Image-based techniques represent intermediate solutions in terms of artistic control and flexibility. Although they may start from 3D-like information such as normal or depth images, they do not require a full simulation of light transport in a 3D scene. For instance, images of textured objects may be created either by segmenting a normal image into patches onto which the texture is synthesized and aligned [START_REF] Fang | Textureshop: texture synthesis as a photograph editing tool[END_REF], or by applying parallax mapping to a depth image integrated from normals [Winnemöller et al. 2009]. Many image-based techniques focus on editing existing images to create new content. For example, [START_REF] Oh | Image-based modeling and photo editing[END_REF]] provided a set of intuitive tools for segmenting, reconstructing and editing depth and color maps from single uncalibrated images enabling a wide range of radical changes, including relighting, removal and copying of content in 3D and viewing from new viewpoints. Image-based material editing [START_REF] Khan | Image-based material editing[END_REF]] is another example: depth and normal images are first estimated from a color image of an object; then its material is modified by means of image processing techniques. The approach provides compelling results but remains limited in the kinds of manipulations it affords: reflections and refractions cannot be positioned freely on the object and one has to work with information found in the image (mainly the background). The method of [START_REF] Yeung | Matting and compositing of transparent and refractive objects[END_REF] provides a more flexible artistic control whereby users can sketch distortion patterns due to refraction, but it is restricted to the matting of transparent objects. The Lit Sphere technique [START_REF] Sloan | The lit sphere: a model for capturing npr shading from art[END_REF]] may also be considered as an image-based shading design technique, although it is strongly limited in terms of artistic control, since it relies on a fixed global mapping from screen-space normal coordinates to the image of a lit sphere.

Overview

Our approach belongs to the image-based category, and like most existing methods it takes normal and depth images as input. Contrary to previous work though, it is not limited to a single shading effect and provides direct artistic control in the picture.

To reach this goal, we first demonstrate in Section 4 that two kinds of surface flows are at the origin of distortion patterns found in rendered images. This is done by conducting a differential analysis of the reflected radiance equation in the vein of [START_REF] Ramamoorthi | A first-order analysis of lighting, shading, and shadows[END_REF], with the important difference that it is performed directly in the picture. This is crucial to our approach, since it explains distortion patterns observed in images of shaded objects in terms of depth and normal variations. We then describe in Section 5 how the gained insights are used in deformation operators to trade-off artistic control versus perceptual plausibility, and we illustrate their usefulness through three types of shading tools. Combining these tools, we produce a variety of shading appearances as shown in Section 6. We discuss limitations and potential extensions to our approach in Section 7.

Surface flows

The main claim of this paper is that it is possible to elicit a powerful illusion of complex shapes or materials without accurately simulating the physics of light transport; one only needs to convey the proper visual cues to the human visual system. But which cues exactly? As mentioned in the introduction, a growing body of evidence suggests that the visual system is particularly sensitive to distortion patterns in the image, which are mostly measurable through spatial variations of luminance. This naturally leads us to perform a differential analysis of the rendering process, with the goal in mind to identify how surface shape participates in the pictorial appearance of a 3D object.

Definitions

We focus on the analysis of reflected radiance Lr, defined at each point (x, y) of the picture plane by:

Lr(x, y) = Ω + n ρn(p, v, ω) Li(p, ω) (n • ω) dω, ( 1
)
where ω is an incoming light direction, n(x, y) is the outward surface normal at the point p(x, y) on the surface, v(x, y) is the view direction and Ω + n is the hemisphere oriented around n. Material characteristics are provided through ρn, the (potentially spatiallyvarying) Bidirectional Reflectance Distribution Function [START_REF] Nicodemus | Geometric Considerations and Nomenclature for Reflectance[END_REF], or BRDF for short. Direct and indirect illumination are encapsulated inside the incoming radiance Li.

Unfortunately, Equation 1 is not directly amenable to differentiation since the domain of integration Ω + n varies itself in the picture plane. This is fixed with a change of variable ← Rnω, where Rn is the smallest rotation that aligns the normal n with the z axis. In the following, we consider only isotropic materials, hence Rn is welldefined. As detailed in the Appendix, Equation 1 then becomes:

Lr(x, y) = Ω + z ρz(p, Rnv, ) Li(p, R -1 n ) (z • ) d , (2)
where Ω + z is the (constant) hemisphere oriented around the z axis, and ρz is the BRDF defined around z. Note that this formulation is still independent of the choice of material, since one only has to replace n by z in the chosen BRDF model.

Differentiation

With this new formulation, we are now able to compute derivatives of Lr as combinations of derivatives of the terms found inside the integral. For ease of notation, we will omit function variables and write Equation 2 as Lr = Ω + z ρzLi dF , where dF = (z • )d is independent of the position in the picture. The derivative of the reflected radiance Lr along the horizontal direction x in the picture plane is obtained by applying the product rule inside the integral:

∂xLr = Ω + z Li ∂xρz dF + Ω + z ρz ∂xLi dF, (3)
where ∂x stands for ∂/∂x. Material variations ∂xρz and illumination variations ∂xLi are obtained by applying the chain rule:

∂xρz = (∂xp • ∇pρz) + (∂xRnv • ∇ν ρz), ( 4
) ∂xLi = (∂xp • ∇pLi) + (∂xR -1 n • ∇ l Li), (5)
where ∇ν ρz is the derivative of ρz in the viewing direction, and ∇ l Li is the derivative of Li in the incident lighting direction. For the sake of simplicity, we have considered that the projection is locally orthographic, and hence ∂xv = 0.

Variations due to surface shape appear in both material and illumination derivatives, and are of two different kinds: they represent either first-order variations (∂xp) or second-order variations (∂xR ±1 n ). Although the meaning of the former is easily grasped (and depends only on depth variations with our local orthographic approximation), the understanding of the latter is indirect since it involves derivatives of the rotation Rn and its inverse. As detailed in the Appendix, we rewrite the derivatives of Rn and R -1 n in terms of cross products that directly involve normal variations:

∂xRnv = Rn∂xn × n × v, (6) ∂xR -1 n = -∂xn × n × R -1 n , (7)
where ∂xn × n corresponds to the axis of the minimal rotation from n to n + ∂xn, and its magnitude represents the corresponding rotation angle. Aside from a more direct interpretation, this formulation affords better segregation between first-and second-order surface structures, as explained in the next section.

Decomposition

The analysis conducted so far has identified two types of variations due to surface shape: first-order (∂xp) and second-order (∂xn×n). We now show how these terms are factored out to yield a decomposition of luminance variations in image space.

We begin by inserting Equations 6 and 7 inside Equations 4 and 5 respectively. The first-order term ∂xp is then directly taken out of the integral. After re-ordering the scalar triple products in Equations 4 and 5, the second-order term ∂xn × n is taken out of the integral as well. This leads to the following compact formulation (see supplemental material for a detailed derivation):

∂xLr = (∂xp • Vp) + (∂xn × n • Vn) . ( 8
)
Vp and Vn give spatial and angular variations respectively:

Vp = Ω + z ∇pρz Li textures + ∇pLi ρz spatial illum. dF, (9)
Vn = Ω + z (v × R -1 n ∇ν ρz)Li reflections -(R -1 n × ∇ l Li)ρz angular illum.
dF. (10) Equation 8 constitutes the main result of this section. It shows that first-and second-order surface variations have distinct effects on luminance variations, and that this distinction is invariant to the choice of materials or lighting environments. So far, we have only considered derivatives in the x direction, but the same analysis holds for the y direction. Taken together, this provides us with a pair of Jacobians [∂xp, ∂yp] T and [∂xn × n, ∂yn × n] T , which we call surface flows to emphasize the distortion effects they impose on other types of 3D variations.

The first-order flow is related to 3D spatial variations, as made explicit by Equation 9: ∇pρz represents variations of BRDF properties across the object (commonly called "texture"), while ∇pLi represents spatial illumination variations due to shadowing or interreflections. The second-order flow is related to 3D angular variations as made explicit by Equation 10: ∇ν ρz represents variations of the BRDF for a varying viewpoint (hence conveying reflections), while ∇ l Li represents angular illumination variations due to the environment. Both Vp and Vn deserve further investigation and should be connected to perceptual findings that have recently been made thanks to experiments on shape-from-texture, shape-fromshadows and shape-from-reflections. However, in this paper, we mainly focus on surface flows since we are interested in the production of believable distortion patterns. We will thus replace Vp and Vn altogether with data coming from images and leave their studies to future work. This solution enables the creation of powerful design operators, as demonstrated in Section 5.

Before we move on to shading design, we need an intuitive way to visualize surface flows in an image. A common solution to represent a Jacobian consists in converting it to a structure tensor [Bigün and Granlund 1987], which corresponds here to a 2 × 2 covariance matrix, from which directions and magnitudes of maximal and minimal change in the image are extracted with an eigen decomposition. The first-and second-order structure tensors Tp and Tn are defined in terms of the directional derivatives of p and n:

Tp = ∂xp • ∂xp ∂xp • ∂yp ∂xp • ∂yp ∂yp • ∂yp ; Tn = ∂xn • ∂xn ∂xn • ∂yn ∂xn • ∂yn ∂yn • ∂yn .
Two important observations must be made. First, the cross product with n does not appear in Tn, as explained in the Appendix. It suggests that this cross product has little effect on the pattern of deformations induced by normal variations. Second, Tp and Tn do not directly correspond to the fundamental forms of Riemannian geometry, since they characterize variations in image space as opposed to variations in tangent space. For a more detailed discussion on this topic, we refer the reader to [START_REF] Fleming | Shape from sheen[END_REF]].

We visualize both types of tensors with a combination of lineintegral convolution (LIC) and a color code as shown in Figure 2.

Visual comparison between compression and stretching patterns in textured or shaded images on the one hand (a), and the representation of their corresponding surface flows on the other hand (b), reveals a strong correlation as expected.

Shading Design

Having identified the surface flows responsible for distortion patterns found in shading and texture, we are now left with a second challenge: how to convey such distortions while still providing direct artistic control? For this purpose, we introduce in Section 5.1 two novel image deformation operators, one for each surface flow.

We then show in Section 5.2 how they might be incorporated into three example drawing tools, each having its own way of conveying additional texture or material properties. Note that the equations presented throughout Section 5 are heuristics inspired by Section 4 and intended to approximate surface flows rather than reproduce them perfectly.

Deformation operators

Let I be the image we wish to deform. We look for a pair of deformation operators which, once applied to I, produce distortion patterns similar to those of first-and second-order surface flows. One may think that Tp and Tn could be used for this purpose, but this is not straightforward: structure tensors only provide non-oriented distortion fields (by construction), and it is not clear how they could be adapted for artistic control since they provide only local information.

Instead, our approach consists in finding user-controllable deformation operators Dp and Dn whose derivatives mimic surface flows. Formally, for a point q = (x, y) in the image, its color after deformation is obtained through composition with one of these operators: I p|n (q) = I(D p|n (q)), where p|n is a shorthand notation for "p or n". The derivative in the x direction of the image deformation process is then given by:

∂xI p|n = (∂xD p|n • ∂xI(D p|n )). (11)
The similarities between Equations 11 and 8 are important. On one side, variations due to deformations ∂xDp or ∂xDn are meant to take the role of ∂xp or ∂xn respectively. On the other side, variations due to the input image ∂xI are designed to mimic 3D spatial (Vp) or angular (Vn) variations. Since our focus is on surface flows, these latter variations will be replaced by variations found in (possibly filtered) images, as shown in Section 5.2.

On top of reproducing 3D-like surface variations, we aim to provide artists with controls over deformation operators. First and foremost, we want to enable direct manipulation through a few anchor points qi, where the deformation is required to have no effect (i.e., D p|n (qi) = qi). Moreover, although the direction of deformation is entirely determined by input normal and depth maps in our implementation, its magnitude should remain controllable with a few parameters. With this trade-off in mind, we have empirically designed the following operators.

First-order deformations Considering a single anchor point q0 for simplicity, we want to locally deform I in a way that preserves its overall structure, in particular its orientation. Consequently, we restrict deformation to a 2D offset of q proportional to (q -q0), with an offset factor fp that depends on the depth difference between q and q0. We use depth instead of full positions for two reasons: depth images are more easily obtained; and depth differences are signed which provides for more intuitive artistic controls.

In the simple case of a single anchor point, the offset vector is given by fp(z -z0)(q -q0), where depth values z and z0 lie in the [0, 1] range and fp(u) = 2u(e γ -1) (e γ + 1) -u(e γ -1) .

As shown in Figure 4(a), fp is a function of relative depth controlled via a single parameter γ ∈ R + . This formula has been designed empirically to provide more natural distortions whenever a distant anchor point is used. Compared to a simple linear function, it prevents extensive distortions on front surfaces (i.e. when z -z0 < 0) while still preserving the overall magnitude of deformation (Figure 5).

Generalizing this idea to multiple anchor points, we obtain the following formula for first-order deformations:

Dp(q) = q + i wi(q)fp(z -zi)(q -qi), (13)
where the wi are used to balance the influence of anchor points. When applied to the whole picture, Dp either inflates (fp < 0) or deflates (fp > 0) the image I. It closely mimics first-order surface variations, as shown in the top rows of Figures 2 and3, where we compare Tp to structure tensors obtained from ∂xDp and ∂yDp.

Second-order deformations Considering once again a single anchor point q0, the second-order operator is obtained with a sim- ilar approach. It consists in offseting a point q in the direction (n -n0) where n and n0 are the normals of q and q0 projected onto the picture plane. The offset factor fn now depends on relative normal variations, and the reason we use their projection is to get a 2D offset direction in the image plane. Considering only q0, the offset vector is given by fn(||n -n0||)(n -n0), with

fn(u) = αu e -β + u(1 -e -β ) . ( 14
)
As shown in Figure 4(b), fn is a function of relative normal variations controlled via two parameters: α ∈ R controls the signed magnitude of the deformation (positive and negative values mimic reflection and refraction patterns respectively), while β ∈ R + controls its rate. The function behaves linearly when β = 0 and nonlinearly when β > 0, hence producing more distortions around forshortened areas, which is useful for some refractive patterns (see Figure 6(e-f)).

As before, generalizing this idea to multiple anchor points, we obtain the following formula for second-order deformations:

Dn(q) = q + i wi(q)fn(||n -ni||)(n -ni). ( 15
)
When applied to the whole picture, Dn distorts the image I in ways that closely mimic second-order surface variations, as seen in the bottom rows of Figures 2 and3. Although our approach is heuristic, in practice the resulting distortions approximate surface flows closely, while maintaining intuitive artistic control.

Figure 6 demonstrates the effects of deformation parameters, here applied to a noise input image with a single anchor point (see also the supplemental video). When used with multiple anchor points, both operators may be seen as weighted sums of offset vectors, one for each anchor point. Provided that we use an interpolating weighting scheme (∀i = j, wi(qi) = 1 and wi(qj) = 0), deformations will have no effect at anchor points. Indeed, f p|n (0) = 0 and hence D p|n (qi) = qi. The choice of weighting scheme depends on selected tools, as described in the next Section.

Drawing tools

Our shading design system also relies on variations found in the input image I to convey surface material (see Equation 11). Operations affecting I, like blurring or color changes, are thus intended to approximate variations due to Vp and Vn. We now introduce shading tools for the control of final image appearance. The brush tool allows users to paint shading in exactly the same way they would paint brush strokes in raster graphics software. It consists in repeatedly applying a brush footprint (I) with a single anchor point in its center; the only difference is that the footprint is deformed before being applied. When used with Dp, it permits the creation of decal textures as shown at the top of Figure 7 (we provide simple scaling and rotation controls of the footprint). When applied with Dn, it becomes easy to create plausible highlights.

Their glossiness can be controlled by varying the blur applied to their footprints. As shown in Figure 1(c) and at the bottom of Figure 7, the tool is effective when multiple footprints are accumulated, hence producing non-physical, yet convincing highlights.

The gradient tool directly takes its inspiration from the classic linear color gradient, except that the gradient image (I) is deformed prior to its application. It makes use of two anchor points q0 and q1, one for each extremity of the gradient, and naturally requires two weighting functions w0 and w1 = 1 -w0. The weight w0 for a point q is obtained by projecting it perpendicularly onto the [q0, q1] segment, and performing a linear interpolation (clamped at extremities). When applied with Dn, different material properties are mimicked depending on the content of I, as illustrated in Figure 8 where we show diffuse, cartoon or metallic appearances.

The image tool allows users to distort a full image (I) to mimic realistic textures or reflections, place highlights or cast shadows.

Examples are shown in Figure 1(b) and Figure 3(b) where we have used a single anchor point. The manipulation of multiple anchor points is particularly useful for controlling the location of particular points of the input image, as shown in Figure 9. We have chosen Shepard interpolants [1968] for the weight functions, but any other interpolation scheme could have been used instead.

Results and comparisons

As shown in Figure 1(d), the final image is obtained by blending the results coming from multiple shading tools. Our implementation makes use of a nodal network that permits any kind of combination, in a way similar to systems found in Nuke or MentalMill for Image-based shading Figure 10 illustrates second-order deformations using normal images coming from different sources. When applied to captured normals, our image tool produces a compelling shading effect with a single drag-and-drop. We have also experimented with a simple enhancement functionality inspired by the Light Warping technique [START_REF] Vergne | Improving Shape Depiction under Arbitrary Rendering[END_REF]; we amplify the mag-nitude of deformations proportionally to their view-centered curvature measure, which yields equivalent results. Our brush tool is the best choice when working with drawn normals: reflections are accurately positioned to create cartoon results. Our method is naturally adapted to edit material appearance in existing images as in [START_REF] Khan | Image-based material editing[END_REF], using estimated normals. The key difference is that our approach provides a more direct and versatile control over apparent shape, materials and illumination effects, for instance to place and move highlights at specific locations on the object.

Texture mapping Figure 11 compares our first-order deformation operator with the texture draping technique of Winnemoeller et al. [2009]. The first important difference is that our approach works directly from a depth image as advocated by our differential analysis, whereas theirs integrates the normal image to get a height field.

Using true depths produces more convincing depictions as best seen around occluding boundaries. The second difference is that our operator is controlled by anchor points located at arbitrary positions in the picture, while theirs relies on parallax mapping, and is thus controlled indirectly by the specification of a view vector.

Complex effects are obtained by combining multiple deformations. For instance in figure 12(a), we add shadows to Figure 1(d) using the brush tool with first-order deformations. In Figure 12(b), we apply a normal texture to create an additional surface normal layer, which is then used as a base for further operations. The translucent object in Figure 12(c) is obtained by combining transparent and specular deformations, blended using Schlick's approximation [START_REF] Schlick | Graphics Gems V. Morgan Kauffman, ch[END_REF]] of the Fresnel term. We also perturbed the input normal map with a noise texture for deforming the structure of refractions, as suggested in [START_REF] Motoyoshi | Highlight-shading relationship as a cue for the perception of translucent and transparent materials[END_REF]].

Material design A variety of material effects is achievable using our image-based shading design approach. This is illustrated in Figure 13, where we have used simple shapes to let the reader better appreciate shading details. Stylized material appearances such as the ones found in logos are trivially obtained with either of our tools, saving a lot of time and training to artists. More photorealistic materials are produced by blurring and deforming real-world photographs. For complex fabric materials, we found it sufficient to layer brush strokes of varying intensity at different locations.

Ease of use Our shading design system is a prototype that is only intended to illustrate our perceptually-motivated approach. Even though it is beyond the scope of the present paper, a proper user study will be useful in the future. Nevertheless, our first experiments are encouraging. Let us first point out that many of the images were created using only a single anchor point, for which the extent of artistic direction of appearance is highly limited. Figure 14(a) provides an extreme example, where we compare our approach with the photorealistic rendering of a mirror object obtained with PBRT [START_REF] Pharr | Physically Based Rendering, Second Edition: From Theory To Implementation[END_REF]. In our result, the original angular environment map is deformed using a single anchor point located at the center of the picture. Slightly more complex examples are provided in Figure 14(b-c): the refractive object is obtained by deforming the same environment map with an opposite magnitude and positioned so that it aligns with the background; whereas the plastic object uses the same configuration as the mirror case, but blurs the environment prior to deformation.

Our second observation is that deformation operators vary appearance within the constraints of a given shape; hence there is a broad "sweet spot" since we exploit key image features that drive shape perception. However, as illustrated in Figure 15, even if deformations remain plausible, unnatural images are obtained when the expected congruence between shading effects is not respected [Kim et al. 2011] (e.g., placing highlights in shadowed regions).

Discussion

Contributions The primary benefits of our approach lie in the speed and directness with which images can be built and edited. Users are free to drag shading, highlights and texture around to edit appearance directly, while preserving apparent shape. Appearance can be honed gradually by adding and tweaking layers, or changed radically in a couple of seconds by swapping between different source textures. Because the results are built out of other images, it is easy to acquire or generate novel source textures. The other advantage is that because the method explicitly relies on the flows that determine perceived shape and material properties, the appearance of these properties can be emphasized, exaggerated or modified directly in the image.

Limitations A minor problem appears when using too many anchor points: shape perception may then be affected since the resulting deformations will differ from plausible variations; on the other side, it will come closer to user desires. Another issue with manually controlled anchor points is animation: when full geometry is given (rather than just normal maps), and objects are rotated and moved, anchor points should move accordingly to yield convincing shading and texture motion. Materials are also controlled with an artistic approach, using simple image processing operators. If plausible renderings can be obtained with this technique, taking the second-order flow into account when applying these operators could lead to more accurate results when it comes to photorealistic appearances. For instance, the current deformation cannot imitate anisotropic materials. With extreme deformations, our system will also tend to produce aliasing. However, we know precisely the extent of deformations at each pixel, which is ideal for designing anti-aliasing strategies.

Future work Although we have only applied deformations to raster images, there is no a priori reason why it would not work with vector graphics. Connecting it with a system that directly outputs vector normal and depth data will likely result in a fine solution for drawing applications (e.g., [START_REF] Shao | Crossshade: Shading concept sketching using crosssection curves[END_REF]). Beyond practical matters, we hope the differential analysis we have conducted in the image plane will foster interest in both the vision and graphics communities. Although our deformation operators closely approximate surface flows, they still deserve further investigations not only to understand their connections with real images, but also to provide new effects, such as anisotropy. In particular, a detailed study of Vp and Vn would provide more insights into how light and material properties interact with shape. They connect variations in a 3D scene to observable pictorial variations, and hence have potential applications not only in rendering and compositing, but also for the study of the human visual system.

Figure 1 :

 1 Figure 1: Our system produces convincing shading results by working purely from image deformations. Here, a surface flow is extracted from a normal image (a). It is first used to deform an entire image that serves as a reflective base coat (b). Then it is employed to deform brush strokes that represent highlights (c). The two components are added together to convey a stylized shading appearance (d).

Figure 2 :

 2 Figure 2: 3D rendering: conventional texture (top) and environment (bottom) mappings (a) produce variations related to first-and second-order surface flows. They are visualized in (b) where hue and saturation indicate direction and magnitude of the flow respectively.

Figure 3 :

 3 Figure 3: Image deformation: our shading design system mimics texture (top) and shading (bottom) effects using image processing alone. It takes depth and normal images as input (a), and uses them to deform images (b) in ways that closely approximate surface flows (c). It provides a convincing, yet artistically controllable illusion of 3D shape conveyed through texture or shading cues.

Figure 4 :

 4 Figure 4: Offset functions: (a) Depth-based offset functions with γ = {1, 2, 3}. The non-linear behavior is chosen to mimic the non-linearity of depth in 3D scenes. Note the invariant fp(0) = 0; (b) Normal-based offset functions with β = {1, 2, 3} for α = 3 (dashed curve) and α = 6 (solid curve). Non-linearity is controlled by the β parameter. Note again the invariant fn(0) = 0.

Figure 5 :

 5 Figure 5: Linear vs non-linear Linear vs non-linear Linear vs non-linear Similar results are obtained with a foreground anchor point (a-b); while a linear function overstretched the texture with a background anchor point (c-d).

Figure 6 :

 6 Figure 6: Effect of deformation operators. In (a), a noise image is anchored to a single point in the picture. In (b-c), it is deformed using Dp p p: progressively increasing γ makes the noise appear as if it were textured on the object. In (d-f) the image is deformed using Dn n n: patterns resembling reflections appear with positive values of α, while negative values imitate refraction patterns; distortions around foreshortened areas are controlled via β.

Figure 7 :

 7 Figure 7: Brush tool. Deformed textures (top row) or shading patterns (bottom row) are applied at arbitrary sizes (red contours) and locations (blue dots). Textures are oriented by the user to position decals. Shading patterns are blurred to convey different materials.

Figure 8 :

 8 Figure 8: Gradient tool. Smooth shading patterns are created by deforming a linear gradient (red curve). Two anchor points (blue and red dots) control its extremities. The gradient image controls material appearance: from diffuse, to cartoon, to metallic.

Figure 9 :

 9 Figure 9: Image tool. A refraction pattern is manipulated using anchor points: (a) the image is positioned (blue dot); (b) the tree refraction is moved downward (green dot); (c) the horizon is aligned with the background. Color insets visualize weight functions.

Figure 10 :

 10 Figure 10: Image-based shading. Image-based shading. Image-based shading. Left: Image tool applied to a captured normal map (RGBN image [Toler-Franklin et al. 2007]), with details enhanced. Middle: Brush tool applied to a drawn normal map. Right: Image-based material editing with our approach: highlights are directly manipulated.

Figure 11 :

 11 Figure 11: Texture mapping: Texture mapping: Texture mapping: compared to [Winnemöller et al. 2009] (a), our approach (b) better conveys discontinuities at occluding contours and produces more plausible deformations.

Figure 12 :

 12 Figure 12: Complex results Complex results Complex results obtained by combining multiple deformations: (a) adding shadows; (b) layering normals; and (c) imitating translucency (see supplemental materials).

Figure 14 :

 14 Figure 14: Comparison with PBRT. Photo-realistic renderings (top) are imitated using the image tool with a single anchor point (bottom) and a few brush strokes for three different materials.

Figure 15 :

 15 Figure 15: Failure case. Varying distortion parameters (here magnitude) often produces plausible shading (a-b), although unconventional anchor point placement leads to unnatural images (c).

Acknowledgments

We would like to thank Adrien Bousseau, Peter Vangorp, Romain Pacanowski as well as anonymous reviewers for their insightful comments on the paper. We are grateful to the AIM@Shape library for their 3D models, and Paul Debevec for his environment maps that we have used in this paper. This research has been supported by the BMBF-NSF joint program on Computational Neuroscience (01GQ1111) and by the ALTA project (ANR-11-BS02-006).

Appendix

Change of variable The BRDF as expressed in Equation 1is dependent on the local normal n. Let's write this dependence more explicitly using the rotation Rn: ρn(p, v, ω) = ρz(p, Rnv, Rnω).

Applying the change of variable ← Rnω to Eq. 1, we obtain:

Noting that Ω + z = RnΩ + n , |Rn| = 1 since Rn is a rotation, and

Derivative of a rotation Let's write the rotation that brings the normal n into alignment with z more explicitly as Rn→z. With the same notation, the first-order approximation of the rotation operator for a neighbor unit normal n = n + ∂xn is given by R n →z = Rn→zR n →n ≈ Rn→z + ∂xRn→z, which yields for an arbitrary vector ν: ∂xRn→zν = Rn→z(R n →n -I)ν.

Let's represent R n →n as a rotation of axis ū and angle θ. To a first-order approximation, we have cos θ = 1, sin θ = ||u||, and ū = u/||u||, with u = n × n = ∂xn × n. Using Rodrigues' formula, the rotation applied to an arbitrary vector v is given by:

Setting v = Rn→zν and plugging this back inside the derivative formula, we obtain ∂xRn→zν = Rn→z∂xn × n × ν. The derivative of R -1 n→z is obtained by taking a rotation axis ū that points in the opposite direction, yielding ∂xR -1 n→z ν = -∂xn×n×R -1 n→z ν.

Second-order structure tensor The structure tensor is a 2 by 2 matrix where each coefficient is obtained as a dot product involving partial derivatives. For the second-order shape term, it also involves cross products with the surface normal n. Lagrange's identity permits expressing each coefficient in terms of dot products only, as exemplified below with off-diagonal coefficients:

This expression is further simplified by noting that the normal is a unit vector (n • n = 1), and thus is always orthogonal to its derivatives (∂tn•n = 0). This yields (∂xn×n)•(∂yn×n) = ∂xn• ∂yn. The derivation of diagonal terms follows the same reasoning.