Jacopo Brunetti
email: jacopo.brunetti@insa-lyon.fr

Francesco Massi

Walter D'ambrogio

Yves Berthier

Numerical investigation on the mode coupling contact dynamic instabilities

Keywords: Frictional Contact, Mode Coupling Instability, Unstable Induced Vibration

In complex mechanical systems the frictional contact is at the origin of significant changes in the dynamic behavior. Unstable vibrations can reach large amplitude that could compromise the integrity of the systems and are often associated with annoying noise emissions. In this paper the numerical analysis of a lumped system constituted by several degrees of freedom in frictional contact with a slider is presented, where the introduction of friction gives rise to an unstable dynamic behavior. Two different approaches are used to investigate the effects of frictional forces: i) a linear Complex Eigenvalue Analysis (CEA) allows for calculating the complex eigenvalues of the system that can be characterized by a positive real part (i.e. negative apparent modal damping); ii) a non linear model has been developed that accounts for the non linear behavior at the contact interface (stick, slip, detachment) to solve the time history solution and analyze the unstable vibrations.

INTRODUCTION

Complex mechanical systems under frictional contact are always subjected to vibrations induced by the frictional contact. Such vibrations can be either of low amplitude, due to the system response at the broadband noise excited at the contact surface, or characterized by an unstable response of the system dynamics (stick-slip, sprag-slip, mode lock-in) and their consequent large amplitude. The presence of frictional contact can give rise to mode-coupling instabilities that produce harmonic "fiction induced vibrations". Unstable vibrations can reach large amplitude that are often associated with annoying noise emission.

The study of this kind of contact dynamic instability has been subject of many numerical [START_REF] Massi | Structural modifications for squeal noise reduction: numerical and experimental validation[END_REF], [START_REF] Coudeyras | A new treatment for predicting the self-excited vibrations of nonlinear systems with frictional interfaces: The Constrained Harmonic Balance Method, with application to disc brake squeal[END_REF] and experimental [START_REF] Massi | Brake squeal: Linear and nonlinear numerical approaches[END_REF] studies on specific mechanical systems, such as automotive brakes, typically affected by such issue. These works have shown how the stable and unstable behavior of the system can be identified by a Complex Modal Analysis (CEA) and that the frequency found to be unstable both experimentally and numerically is close to one of the unstable modes resulting form the CEA. Theoretic and numerical analyses of simple lumped models considered mainly systems with no more than one unstable mode to observe the effect of some parameters on the complex eigenvalues of the system [START_REF] Sinou | Methods to reduce non-linear mechanical systems for instability computation[END_REF] or on size of the limit cycle of oscillations [START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF]. Often, contact non linearities have been replaced by non linear stiffness in order to apply non linear resolution methods [START_REF] Sinou | Methods to reduce non-linear mechanical systems for instability computation[END_REF]. By the analysis of the contact behavior of complex mechanical systems it can be observed that the variation of contact status strongly affects the limit cycle; periodical variations of contact status can be found locally at the contact [START_REF] Tonazzi | Instability scenarios between elastic media under frictional contact[END_REF], affecting the maximum amplitude of vibration reached by the mechanical system (limit cycle).

In this paper, the numerical analysis of a lumped system constituted by several degrees of freedom in frictional contact with a slider is presented, where the introduction of friction can give rise to an unstable dynamic behavior. The non linear effects of the contact are in this case accounted for by considering the possibility of switching between sliding, sticking and detachment of the bodies in contact. Two approaches are used to investigate the effects of friction forces. The CEA allows for calculating the complex eigenvalues of the system that can be characterized by positive real part (i.e. negative apparent modal damping). The effects of the main parameters on the system stability are investigated. In the second approach, a non linear model has been developed to account for the real contact behavior at the interface to calculate the transient solution and analyze the unstable vibrations. The transient response of the system is directly related to the results from the CEA.

LUMPED MODEL

The analyzed system is composed by 4 masses m i , where 2 masses are in contact with rigid sliders with a constant friction coefficient µ (Fig. 1). The sliders are inclined at the angle θ i and move with speed v p,i .

The two masses belonging to the same vertical line are constrained to have the same horizontal position. External forces F ex,1 and F ex,2 applied to the upper masses m 1 and m 3 ensure a preload at the contact. The contact between the sliders and the lower masses m 2 and m 4 acts as an unilateral constraint (y 2 ≥ 0 e y 4 ≥ 0) and the reaction forces can be expressed as follow:

             N 1 = 0, T 1 ≤ -µ ẋ2 -v p ẋ2 -v p N 1 if y 2 = 0 N 2 = 0, T 2 ≤ -µ ẋ4 -v p ẋ4 -v p N 2 if y 4 = 0 N 1 = 0, T 1 = 0 if y 2 > 0 N 2 = 0, T 2 = 0 if y 4 > 0 (1)
Referring to the DoFs introduced in Fig. 1, the equations of motion of the system can be expressed as follows:

                                               -m 1 ẍ1 -c 1 ẋ1 + c 3 ( ẋ3 -ẋ1 ) -k 1 x 1 + k 3 (x 3 -x 1 ) = F 1 -m 1 ÿ1 -c 2 ( ẏ1 -ẋ2 sin θ 1 -ẏ2 cos θ 1 ) -c 5 ẏ1 + c 6 ( ẏ3 -ẏ1 ) -k 2 (y 1 -x 2 sin θ 1 -y 2 cos θ 1 ) -k 5 y 1 + k 6 (y 3 -y 1 ) = F ex, 1 -m 2 ẍ2 + c 2 ( ẏ1 -ẋ2 sin θ 1 -ẏ2 cos θ 1 ) sin θ 1 +k 2 (y 1 -x 2 sin θ 1 -y 2 cos θ 1 ) sin θ 1 = -F 1 cos θ 1 -T 1 -m 2 ÿ2 + c 2 ( ẏ1 -ẋ2 sin θ 1 -ẏ2 cos θ 1 ) cos θ 1 +k 2 (y 1 -x 2 sin θ 1 -y 2 cos θ 1 ) cos θ 1 = F 1 sin θ 1 -N 1 -m 3 ẍ3 -c 3 ( ẋ3 -ẋ1 ) -k 3 (x 3 -x 1 ) = F 2 -m 3 ÿ3 -c 4 ( ẏ3 -ẋ4 sin θ 2 -ẏ4 cos θ 2 ) -c 6 ( ẏ3 -ẏ1 ) -k 4 (y 3 -x 4 sin θ 2 -y 4 cos θ 2 ) -k 6 (y 3 -y 1 ) = F ex,2 -m 4 ẍ4 + c 4 ( ẏ3 -ẋ4 sin θ 2 -ẏ4 cos θ 2 ) sin θ 2 +k 4 (y 3 -x 4 sin θ 2 -y 4 cos θ 2 ) sin θ 2 = -F 2 cos θ 2 -T 2 -m 4 ÿ4 + c 4 ( ẏ3 -ẋ4 sin θ 2 -ẏ4 cos θ 2 ) cos θ 2 +k 4 (y 3 -x 4 sin θ 2 -y 4 cos θ 2 ) cos θ 2 = F 2 sin θ 2 -N 2 (2)
Internal constraints can be expressed as:

x 1 = x 2 cos θ 1 -y 2 sin θ 1 x 3 = x 4 cos θ 2 -y 4 sin θ 2 (3)
so that the system is reduced to 6 DoFs (y 1 , x 2 , y 2 , y 3 , x 4 , y 4 ). To analyze the nonlinear behavior of the system the following status of contact have been considered:

Sticking ( ẋt = v p ) ;
Negative Sliding ( ẋt -v p < 0): sliding contact with negative relative speed and positive friction force;

Positive Sliding ( ẋt -v p > 0): sliding contact with positive relative speed and negative friction force;

Detachment (x n > 0): no contact interaction between the mass and the slider.

These conditions can be reached independently on both the contacts between masses and sliders,and the system can reach 4x4 possible combinations of the contact status.

The simulations here reported are performed with the following system parameters where masses are expressed in [kg] and stiffness in [N/m]:

m 1 = 37.50 m 2 = 0.25 m 3 = 3.75 m 4 = 1.50 k 1 = k 2 = k 3 = 2.0e+5 k 4 = 2.3e+5 k 5 = k 6 = 2.0e+4. C = βK β = 5e-5 s (4)

COMPLEX EIGENVALUE ANALYSIS

To find the eigenmodes of the system the equilibrium position is accounted for. It is supposed that all the contacts are in sliding condition. In this case the motion is constrained in the direction witch is normal to the slider. The reaction force N generates a tangential force T = µN whose orientation agrees with that of the rigid slider speed.

In sliding condition the system has 4 DoFs. The Damping matrix is introduced as to be proportional to the stiffness matrix. If the friction coefficient is nil the modes of the system can be expressed solving the eigenvalue problem for the mass and stiffness matrices.

λMΨ = KΨ (5)
where λ = ω 2 are the real eigenvalues of the system and Ψ are the eigenvector that represent the modal shapes (Fig. 2). The modes at lower frequency correspond to a relevant motion of the left subsystem in vertical direction and to a in phase motion of the two subsystems in horizontal direction while the modes at higher frequency correspond to a relevant motion of only the right subsystem.

If the friction coefficient is different from zero, it produces an asymmetry on the stiffness and damping matrices. The system modes can be found in this case by the means of a complex modal analysis. The system can be expressed in the state space:

Mẍ + C ẋ + Kx = F M ẋ -M ẋ = 0 (6) 0 M M 0 ẋ ẍ + K C 0 -M x ẋ = F 0 (7) A = 0 M M 0 B = K C 0 -M (8) λAξ = -Bξ (9)
where λ are in this case the complex eigenvalues and ξ are the complex eigenvectors of the system. The eigenvalues can be expressed as a function of the angular frequency ω and the damping factor ζ of the mode:

λ = -ωζ + iω 1 -ζ 2 . ( 10
)
The real part is correlated to the modal damping factor: a positive real part corresponds to modes with apparent negative modal damping factor (unstable modes). When the effects of friction forces increase, eigenvalues have initially different frequency coalesce together (lock-in) and the real parts start to diverge (Hopf bifurcation point). Figure 3 reports the trend of real and imaginary parts of each eigenvalue with respect to the friction coefficient µ. Eigenvalues at lower frequency are the firsts to coalesce for a critical friction coefficient µ I = 0.35. The two modes at higher frequency coalesce for a friction coefficient µ II = 0.63. It can be noticed that positive real part of the lower eigenvalues overcome the real part of the higher eigenvalue for a friction coefficient between 0.35 and 0.70. The opposite occurs for µ > 0.70.

Complex eigenvectors give information about both the relative amplitude of vibration and the initial phase of the vibration of each degree of freedom. Each pair of complex conjugate modes gives the following (real) contribution to the displacement solution:

x(t) = ψ i e λit + ψ i e λ i t = 2 Re ψ i e λit (11)
The orbits of the masses of the system allow for a physical interpretation of the different eigenvectors both in terms of vibration amplitude and in terms of the phase delay between the different DoFs of the system (Fig. 4). Convergent spirals correspond to stable modes while divergent spiral correspond to unstable modes. For unstable modes the counterclockwise verse produces a variation of the normal and tangential force that agrees with the tangential speed variation giving a power absorption of the system. On the contrary, in stable modes the clockwise verse produces a variation of the contact forces that is in opposition with the tangential speed variation and this gives a power dissipation of the system.

TRANSIENT NON LINEAR ANALYSIS

When dealing with linear systems, the solution can be expressed as a linear combination of all the modes of the system. Basing on the real part of the eigenvalues, the system behavior can be stable or unstable. For linear systems, if there are several unstable modes, the mode with greater real part of the eigenvalues grows exponentially faster than the others, and the contribution of the other modes becomes quickly irrelevant. Real systems don't behave linearly because of the contact nonlinearities that can be the change of contact status (stick, slip, detachment) or the nonlinear contact stiffness. When the vibration amplitude increases, the variation of the contact status (detachment or sticking) confines the exponential growth of the vibration amplitude. This is why the linear model is fully representative only of the initial part of this dynamic instability phenomena. To simulate in a reliable way the whole transient analysis in the case of unstable friction induced vibrations, a non linear model is needed.

In this paper transient analyses have been performed to observe the time evolution of the nonlinear system in presence of contact forces considering the nonlinearities introduced by the contact itself. The response of the system has been calculated by a numerical integration of the differential equations of the system starting from an initial condition close to the equilibrium position.

Starting from the contact sliding condition and integrating in time, different scenarios can occurs for each contact pair: i) if the normal contact force reaches the nil value, the mass leaves the slider and pass to the detachment condition; ii) if the relative speed of the mass with respect to the slider reaches the nil value the system can either pass to the sticking status or overcome the nil value with the consequent inversion of relative motion (sign of the friction force).

For each combination of contact status of the contact pairs a different set of differential equations can be written; each time a switch condition is reached, the position in the state plane of the system becomes the initial condition of the new differential set of equations.

Effects of friction coefficient on transient behavior.

With reference to the different behaviors of the system observed by the complex eigenvalue analysis, transient nonlinear simulations have been developed in function of the friction coefficient (see Fig. 3): µ < 0.35: the system is stable; 0.35 < µ < 0.63: the system is unstable because one of the eigenvalue at lower frequency has a positive real part; 0.63 < µ < 0.70: the system is unstable because there are eigenvalues at both lower and higher frequency with positive real part. The real part of the eigenvalue with lower frequency is greater than the other; µ > 0.70: the system is unstable because there are eigenvalues at both lower and higher frequency with positive real part. The real part of the eigenvalue with higher frequency is greater than the other.

In case of stable behavior of the system, the amplitude of vibration decreases with exponential decay factor. If the initial amplitude of vibration is less than the amplitude that brings to the contact status variation, the system behaves linearly and tends exponentially to reach its static equilibrium position.

In case of unstable behavior of the system, the amplitude of vibration, initially induced by a little perturbation, increases exponentially until the first mass reaches one of the possible status variation. Then the exponential increase of the vibration amplitude stops and the system vibrations reaches its limit cycle. In figure 5, a typical behavior of the system when at least one of the modes is unstable is reported. Results reported in figure 5 are referred to a friction coefficient µ = 0.65. The initial perturbation is

ẏ1 (0) = 1e -3 m/s.
The state diagrams of the y 1 and y 3 coordinate highlight the limit cycle of the system. The complexity of the state diagram of the y 3 coordinate is due to the rational ratio 9/29 between the two frequencies f 1 = 12.70 Hz and f 2 = 40.93 Hz.

Results of three different analyses with three different values of the friction coefficient and the same initial perturbation are reported in figure 6. The DFT of the system vibrations in the initial part of the simulation show that when the system behaves linearly, the excited frequencies are exactly the same of the complex modal analysis (Fig. 3). Their values increase slightly with the friction coefficient.

In the cases presented in figure 6 the initial perturbation acts on the degree of freedom y 1 and it excites (initially) mainly the first mode at a frequency of about 12 Hz. Increasing the friction coefficient, the peak amplitude corresponding to the second mode (at about 40 Hz) increases faster than the amplitude of the first mode according to the real part trend of the linear parametric analysis (fig. 3).

Although in the first part of the simulation the system behaves differently for different friction coefficient, due to the different grow rate of the two unstable modes, in the second part of the simulation the limit cycles are very similar in spite of the different value of µ.

The frequencies of vibration calculated by the DFT in this part of simulation are higher than the frequencies calculated during the first (linear) part because of the stiffening of the system due to the presence of sticking motion that acts as a constraint on the direction tangential to the contact.

Furthermore, for increasing friction coefficient the amplitude of vibration of the limit cycle increases slightly. When the friction coefficient is µ = 0.50 the amplitude of displacement oscillation measured on the y 1 direction (see fig. 1) is ∆y 50 1 = 0.0398; for µ = 0.65 the amplitude of vibration is ∆y 65 1 = 0.0412 and for µ = 0.75 the amplitude of vibration is ∆y 75 1 = 0.0417.

Effects of the initial condition.

Figure 7, shows the DFT of the system response for different perturbations (initial conditions), both in amplitude and in shape. Result highlighted that different initial perturbations affect only the initial part of the simulation, when the system behaves linearly. The first nonlinearities to appear can also depend on the shape and amplitude of the initial perturbation that may be more similar to an unstable mode rather than to another. In this case, the mode with a shape that is more similar to the initial condition gives to the system response an initial greater contribution. Sometimes if the contact nonlinearities arise when the dynamic of the system is still affected by the initial perturbation, the relative amplitude of the peaks found (in the linear part) do not agree with the real parts of the respective eigenvalues. On the contrary the steady state reached by the system (DFT in the right graph in Fig. 7) is not at all affected by the initial perturbation, confirming that the harmonic content and the behavior of the system during the limit cycle is directly correlated to the system parameters.

Effects of vertical static load on contact nonlinearities.

Figure 8 shows the results of three simulations with different external static load F ex1 and F ex2 . The rest of the system parameters, boundary conditions and initial conditions are kept constant. Friction coefficient for these simulation is µ = 0.75 and the initial perturbation is ẏ1 (0) = 1e -3 m/s. One of the main effects of the variation of the external load concerns the amplitude of vibration during the limit cycle. The amplitude variation is strictly related to the periodic switches of the different statuses of the contact pairs.

For high values of preload (F Ex = 7000 N and F Ex = 5000 N ), the contact status alternates between the sticking status and the sliding status (leftmost and central graphs in figure 9) and this corresponds to a very small variation in terms of harmonic content. Horizontal portions of plots dx 2 /dt and dx 4 /dt indicate a sticking status for the mass m 2 and/or m 4 : in this case the dx 2 /dt and/or dx 4 /dt are equal to the speed of the rigid slider. On the contrary, other portions denote a sliding status. In the leftmost and in the central graph of figure 9 no detachment occurs and dy 2 /dt and dy 4 /dt are zero.

Reducing the external load (F Ex = 3000 N ) some detachments of mass m 2 appear during the limit cycle and are indicated by portions of the plot of dy 2 /dt that are different from zero (rightmost graph in figure 9).

In this last case, both the amplitude and the width of the peak at the lower frequency are reduced (figure 8) that means a lower energy content of the system around that frequency.

CONCLUSIONS

In this paper, the development and analysis on a lumped parameters non linear model, composed by several degrees of freedom, is presented. The presence of coulombian contact with a constant friction coefficient can make the system unstable. By a parametric Complex Modal Analysis on the linearized model, the effects of the variation of the friction coefficient on the system stability have been highlighted. System configurations with either one or several unstable modes are recovered.

The non linear transient analyses, considering the possibility of the contact status variation at the contact (switching between sliding, sticking or detachment), allow for finding limit cycles of the system by accounting for the contact nonlinearities. The effects of boundary and initial conditions on the initial part (linear) and on the limit cycle (non linear) of the system response have been analyzed.

Results reported in this paper show that when the friction coefficient is in the instability range, variations of the value produce significant variations only in the initial part of the transient simu- lation, where the system behaves linearly and on the beginning of the non linear behavior. In the limit cycle the same harmonic content can be found and only variation of the vibration amplitude are noticed in the case of variation of the friction coefficient. On the contrary, the external static loads applied to the system modify the contact status during the limit cycle and relevant variations of the harmonic contents can be noticed.

2 Figure 1 :

 21 Figure 1: Lumped system considerd for calculations.

Figure 2 :

 2 Figure 2: Modal shapes in case of µ = 0.

Figure 3 :

 3 Figure 3: Real and Imaginary part of the eigenvalues versus the friction coefficient

Figure 4 :

 4 Figure 4: Complex modal shapes for a friction coefficient µ = 0.75

  028 -0.027 -0.026 -0.025 -0.024 -0.023 -0.022 -0.021 -0

Figure 5 :

 5 Figure 5: Transient results. µ = 0.65

Figure 6 :

 6 Figure 6: Discrete Fourier Transform of the initial (linear) part and of the final (steady state) part of the system response signal for different friction coefficients

dy 1 / 3 /

 13 dtm(0)=1e-3m[m/s] dy 1 /dtm(0)=1e-4m[m/s] dtm(0)m=m1me-4m[m/s]

Figure 7 :

 7 Figure 7: Discrete Fourier Transform of the initial part and of the steady state part of the system response for different initial perturbation.

Figure 8 :

 8 Figure 8: Discrete Fourier Transform of the steady state part of the system response for different external static load.

Figure 9 :

 9 Figure 9: Velocity response during the limit cycle. Graphics from left to right are referred to F ex = 7000 N, F ex = 5000 N and F ex = 3000 N


Boris Raymond

Gaël Guennebaud

Pascal Barla

Multi-Scale Rendering of Scratched Materials using a Structured SV-BRDF Model

Keywords:

a) Parallel scratches (BRDF) (b) Parallel scratches (SV-BRDF) (c) Deep layered scratches (SV-BRDF) (d) Thin layered scratches (SV-BRDF) Figure 1: Our approach is tailored to the high angular and spatial frequencies exhibited by scratched materials. In (a) we show only angular variations via our scratch BRDF model, which imitates the appearance of brushed metal. A more realistic look is obtained in (b) by incorporating spatial variations thanks to our SV-BRDF model. Changing the profile and distribution of scratches in (c,d) yields different material appearances. Our model reproduces the characteristic look of scratched metals at multiple scales (insets show 3D camera zooms). This includes directionally-smeared reflections in (a,b), patterned highlights in (b), glint lines in (c,d) and hazy environment reflections in (d).

Introduction

Research on material appearance is pushing forward the development of ever-sophisticated and physically-plausible models. Traditionally in computer graphics, opaque materials are represented using Bidirectional Reflectance Distribution Functions (BRDF). For a homogeneous material, a single BRDF instance might be used to render an entire object. However, real-world materials are seldom, Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm. if ever, perfectly homogeneous: they often exhibit a lot of spatial micro-variations, in either reflectance or geometric details.

In this work, we focus on materials exhibiting many tiny scratches, as found mostly on metals, but also plastics or finished woods. The scratches may be due to polished or deteriorated objects, and are usually distributed all over the surface. As shown in Figure 2, they produce a number of interesting visual effects that are visible even from afar: glint lines around highlights, patterned highlight silhouettes, and directionally-smeared reflections.

As recently demonstrated by [START_REF] Yan | Rendering glints on high-resolution normal-mapped specular surfaces[END_REF], homogeneous BRDF models are not adapted in these cases; hence the modeling and rendering of high-frequency surface details remains a challenging problem. For extreme close-ups, scratches are individually visible in all their details. When getting farther away, reflections from individual scratches may still remain visible even though they are much smaller than a pixel: they lead to glint lines. At a distance, the distribution of scratches still has a visible impact on appear-ance: it modulates highlight silhouettes and smears environment reflections. Scratched materials thus require a BRDF model that is not only spatially-varying, but also multi-scale.

The key insight of our work is to exploit the regular microstructure of scratched materials. Individual scratches are mostly one-dimensional, and many of them are distributed parallel to each other. Complex scratch distributions (Figure 2-left) exhibit a combination of parallel scratches in different orientations, whereas simpler ones (Figure 2-right) only consist of a single set of orderly parallel scratches. In the following, we call a scratch layer a set of similar parallel scratches; a scratch distribution becomes a combination of one or more layers. This is not to be confused with the concept of layered materials (e.g., [Jakob et al. 2014a]): here we are only concerned with the outermost, scratched interface.

Our approach is to model scratched materials by combining the reflectance properties of a myriad of individual scratches, exploiting the regularity of their distribution to tackle the high dimensional complexity of the problem. We make the following contributions:

• We introduce a BRDF model for scratches, which, contrary to most previous approaches, explicitly deals with interreflections inside a scratch, including Fresnel effects.

• We extend it to a multi-scale spatially-varying BRDF model (SV-BRDF) for the rendering of scratched materials, using arbitrary layered scratch distributions.

• We provide controls at the scratch level: users may adjust their profile, Fresnel coefficients, width, height, density and orientation, while observing results interactively.

• We describe how to evaluate our SV-BRDF in arbitrary settings, either using importance or light sampling.

• We demonstrate our model in progressive rendering for editing and with global illumination for high-quality results.

The main appeal of our approach is to provide a lightweight, fast and controllable solution reproducing the main visual effects of scratched materials, as demonstrated in Figure 1. The accompanying video shows that our method reproduces the characteristic motion of reflections on such materials: strongest highlights are enriched with glint lines that remain coherent in space and time.

Previous Work

Appearance modeling is a very wide topic, hence we focus here on approaches that permit to reproduce scratched materials.

Distribution-based models

Most BRDF models rely on the micro-facet theory [START_REF] Torrance | Theory for offspecular reflection from roughened surfaces[END_REF], where the shape of a BRDF is mainly determined by the normal distribution function (NDF) of micro-facets.

Reproducing the appearance of scratched materials requires an anisotropic BRDF, and hence an anisotropic NDF. To this end, [START_REF] Ward | Measuring and modeling anisotropic reflection[END_REF] first proposed to use an elliptical Gaussian NDF. Then, more complex anisotropic distributions including an additional Fresnel term have been introduced (e.g., [START_REF] Ashikhmin | An Anisotropic Phong BRDF Model[END_REF][START_REF] Dong | Predicting appearance from measured microgeometry of metal surfaces[END_REF]). These analytic models are fast to evaluate but they are valid only for very distant views, assuming that the structure of the micro-geometry is never to be distinguished.

More recent work has attempted to extend distribution-based BRDF models to handle multiple scales, starting from a detailed normal map or height map. A first body of methods [START_REF] Tan | Multiresolution reflectance filtering[END_REF][START_REF] Tan | Filtering and rendering of resolution-dependent reflectance models[END_REF][START_REF] Han | Frequency domain normal map filtering[END_REF]] builds a multi-modal NDF from the input normal map. However, the presence of multiple modes involves non-linear interpolation across scales to ensure smooth transitions.

A second family of approaches makes use of a unimodal anisotropic NDF which is updated at increasing scales. For instance, [START_REF] Bruneton | Realtime realistic ocean lighting using seamless transitions from geometry to brdf[END_REF] use trochoid averaging for the special case of ocean surfaces, whereas [START_REF] Olano | Lean mapping[END_REF] and [START_REF] Dupuy | Linear efficient antialiased displacement and reflectance mapping[END_REF] make use of surface slope moments working for arbitrary but static input displacement maps.

Regarding scratched materials, distribution-based models are limited in two respects. First, they only consider single-scattering and include, at best, approximate shadowing and masking terms. However, we anticipate that inter-reflections inside a scratch impact the visual material appearance. Second, a smooth NDF is valid only when a large number of scratches fall below a pixel. Hence costly sub-pixel sampling is required for intermediate scales, otherwise details will be over-smoothed during transitions.

Data-driven models

A different approach to the modeling of materials with spatial variations relies on dense and discrete data that can either be measured or computer-generated. The most general approach is the 6D Bidirectional Texture Function (BTF), which depends on two additional spatial coordinates compared to a 4D BRDF. [START_REF] Ma | Level-of-detail representation of bidirectional texture functions for real-time rendering[END_REF] demonstrate multi-scale renderings of BTFs by means of a Laplacian pyramid. However, owing to their high-dimensionality, BTFs are restricted in terms of spatial and angular resolutions, which is problematic to reproduce the high frequency of scratched materials. Anisotropic BRDF measurements [START_REF] Filip | Effective acquisition of dense anisotropic brdf[END_REF]] already exhibit reduced angular resolutions; measuring high-quality BTFs of scratched materials would be even more difficult.

A number of alternatives to BTFs has been recently proposed. [START_REF] Wu | Characteristic point maps[END_REF][START_REF] Wu | Physicallybased interactive bi-scale material design[END_REF] introduced two of them, one based on characteristic point maps, and the other one on pre-computed matrices.

Both deal with arbitrary micro-geometry. The latter even supports dynamic edition while obtaining renderings for distant views at interactive rates. However, they require costly pre-computations and induce large memory footprints limiting the angular resolution too much to be adapted to scratched materials. [START_REF] Yan | Rendering glints on high-resolution normal-mapped specular surfaces[END_REF] proposed a hierarchical pruning technique to efficiently evaluate the NDF found under each pixel footprint from very high-resolution normal maps. The data-driven NDF is then plugged in a microfacet model. To the best of our knowledge, this is the first method to produce convincing renderings of scratched materials, exhibiting vivid glint lines. However, this approach ignores inter-reflections inside scratches and requires very high resolution maps to faithfully represent all the variations of the micro-geometry. Very recently, [START_REF] Dong | Predicting appearance from measured microgeometry of metal surfaces[END_REF] studied the generation of effective NDFs from measured micro-geometry of metals. Their main contribution is in the prediction of brushed metal BRDFs using data-driven NDFs from profilometer measurements. They also present a compact ellipsoidal NDF model to achieve 2D variations from 1D measured exemplars. This latter approximation smoothes out fine variations and shares the limitations of other distribution-based models (i.e., no inter-reflection, approximate shadowing-masking).

Primitive-based models

A third approach consists in modeling materials based on distributions of small primitives. [START_REF] Poulin | A model for anisotropic reflection[END_REF] pioneered the technique by using distributions of micro-cylinders, yielding anisotropic BRDFs. As with other BRDF models, inter-reflections are ignored and micro-geometry is never visible. Individuallyvisible scratches have been investigated for close-up render- ing [Mérillou et al. 2001;[START_REF] Bosch | A Physically Based Model for Rendering Realistic Scratches[END_REF][START_REF] Bosch | Real-time path-based surface detail[END_REF]. However, these methods do not deal with inter-reflections inside scratches, nor do they consider distributions of scratches for rendering in more distant views.

Cloth-rendering has also motivated the development of material models based on fiber primitives. [START_REF] Pont | Split off-specular reflection and surface scattering from woven materials[END_REF] have first demonstrated how simple sinusoidal fiber profiles lead to complex split off-specular reflections and sophisticated grazing angle effects. However, their analysis is restricted to the plane of incidence and ignores inter-reflections. More recent work [START_REF] Irawan | Specular reflection from woven cloth[END_REF][START_REF] Sadeghi | A practical microcylinder appearance model for cloth rendering[END_REF] provide complete descriptions of BRDF models from fiber properties. They rely on the mirror reflection cone first introduced by [START_REF] Kajiya | Rendering fur with three dimensional textures[END_REF] and ignore inter-reflections between fibers. [START_REF] Irawan | Specular reflection from woven cloth[END_REF] also provide a texture model, even though multi-scale rendering with transitions between their BRDF and texture models is not demonstrated.

Overview

Our approach is tailored to the modeling and multi-scale rendering of scratched materials. We introduce a specific solution to deal with the high angular and spatial frequencies typical of such materials, without restricting the choice of scratch profiles or distributions.

The different steps of our method are illustrated in Figure 3.

Our BRDF model is both primitive-based and data-driven, as described throughout Section 4. The BRDF for a single scratch is simulated using a 2D ray-tracing algorithm (Figure 3a), including interreflections and realistic Fresnel effects. The result is stored in the three chanels of a 2D texture (Figure 3b) thanks to a parametrization inspired by the mirror cone of [START_REF] Kajiya | Rendering fur with three dimensional textures[END_REF].

The SV-BRDF (Figure 3c) is obtained as a linear combination of scratch BRDFs (bottom) driven by the area of scratches under a pixel or ray footprint (top). The rendering of our scratched materials at multiple scales (Figure 3d) is covered in Sections 5 and 6, where we detail evaluation methods using either light or importance sampling, in both direct and global illumination settings.

Scratched material model

Hypothesis We define scratches as extruded 1D profiles that are dug into a base surface, which is considered unaltered in-between scratches. As in previous work, we assume the base surface is locally flat, with scratches in a same layer running straight and parallel to each others. In other words, surface curvatures and scratch direction variations (splay and geodesic curvature) are considered negligible at the scratch scale. These conditions ensure that all material-level inter-reflections will be confined within a single scratch. We explain in Section 5.3 how to relax our hypothesis to handle scratch profiles exhibiting bulges. We further restrict the micro-BRDF inside a scratch to be a Fresnel reflector. However, the micro-BRDF outside a scratch, which corresponds to the base material, may be chosen freely.

As demonstrated in our results, these assumptions are reasonable for the modeling of scratched materials, enabling a wide range of scratch profiles, distributions and appearances. Most importantly, they allow us to devise a compact and efficient SV-BRDF model as explained in the following.

Mirror scratch BRDF

Since inter-reflections are confined to a single scratch, each scratch may safely be considered independently of others. We thus start by building a BRDF model for a single mirror scratch extruded from an arbitrary profile. As illustrated in Figure 4a, an important property of such an ideal scratch is that light directions contributing to reflection lie on a half-cone around the scratch direction u. This property was first observed by [START_REF] Kajiya | Rendering fur with three dimensional textures[END_REF] for fibers. The aperture of this cone depends on viewing elevation. Similarly, we identify the family of cones by mapping any unit vector ω ω ω in the hemisphere to angular coordinates (θ, φ) ∈ [-π 2 , π 2 ] 2 given by:

(θ, φ) = (sin -1 (ω ω ω • u), atan2(ω ω ω • v, ω ω ω • n)),
where n is the surface normal, and v is a vector orthogonal to both n and u, hence forming a surface tangent frame.

Let ρm(θo, φo, θi, φi) denote our scratch BRDF in this parametrization, where the subscript m stands for mirror, and the subscripts o and i refer to the outgoing and incoming directions respectively. The aforementioned half cone of reflection implies that for a given view direction, only lights on the iso-curve θ = -θo will contribute (purple curve in Figure 4b). In other words, ρm will be zero whenever θi = -θo, and the dimension of our BRDF can be reduced from 4 to 3. Moreover, since a scratch consists of an extruded profile, varying the view or light elevations along the scratch direction (i.e., θo or θi) will not change the BRDF behavior along the respective cone. Therefore, ρm is invariant with respect to θ, and our BRDF simplifies to a 2D function:

ρm(θo, φo, θi, φi) = ρm(φo, φi) if θi = -θo, 0 otherwise. ( 1
)
This dimensionality reduction allows us to quickly pre-compute our mirror scratch BRDF model within a compact and discrete M × M 2D grid parametrized by φo and φi (see Figure 4c) that we fill using ray-tracing and density estimation. The unique input of our algorithm is a scratch profile function for which some examples are depicted in Figure 5. For each discrete value of φi, we shoot N ≫ M rays at random positions over a scratch, perform recursive ray-tracing with mirror reflection within the scratch until the ray exits the scratch with a direction φo. We accumulate its relative integrated quantity of radiation M N cos(φo) into the respective grid entry, where cos(φo) stands for the classic cosine term in our parametrization. Thanks to the above properties with respect to θ, all computations are performed in the 2D plane supporting the scratch profile (see Figure 3a for an illustration), that is for θi = -θo = 0.

By construction, our mirror scratch BRDF ensures energy conservation as well as Helmholtz reciprocity up to numerical accuracy. This last property implies a symmetry along the diagonal of the 2D domain. For symmetric scratch profiles, the BRDF also exhibits a symmetry along the anti-diagonal. We exploit these two symmetries to get up to 4 times the number of accumulated rays for free.

Figure 5 depicts the behavior of ρm when changing the input scratch profile. We visualize the impact of the number of interreflections using red, green and blue colors for the first, second and higher-order bounces respectively. Steeper scratch profiles lead to an increased number of inter-reflections, while flatter profiles tend toward perfect mirror reflection (the anti-diagonal where φi = -φo, in red). Measured profiles result in more complex BRDFs; in particular, they are not symmetric across the anti-diagonal.

Full scratch BRDF

The scratch BRDF of the previous section is an idealized model since no real-world material comes anywhere close to a perfect mirror. We thus now explain how to include Fresnel effects inside a scratch to obtain a more physically-plausible scratch BRDF.

In the micro-facet theory, Fresnel effects are easily included via a single multiplicative term. This is because only single scattering is considered, with additional inter-reflections removed by shadowing. In our case, computing the Fresnel effect is necessarily more involved because intensity along a light path must be attenuated at each specular reflection according to Fresnel formula. Since this attenuation depends on the angle θ made between the ray and the micro-surface normal in 3D, a scratch BRDF with Fresnel effects can no longer be represented exactly with a 2D function.

To get a better idea of the influence of the Fresnel effect, we have simulated the full 3D scratch BRDFs with Fresnel attenuation using either the approximation of [START_REF] Schlick | An inexpensive brdf model for physicallybased rendering[END_REF] for dielectrics, or the approximation of Laznyi et al. [2005] for conductors. In practice, we compute a full scratch BRDF ρs(θo, φo, φi), where s simply stands for scratch, using our 2D ray-tracer for a range of θo values and storing the results in a 3D grid. To investigate variations in the θo dimension, we plot reflectance curves corresponding to fixed (φo, φi) pairs in Figure 6-left, where we have linearly remapped reflectance values to the [0, 1] range. Here we have used the approximation of Laznyi et al., but Schlick's approximation gives very similar results. As shown in Figure 6-right, the remapped reflectance curves are well approximated by a function of the form Rγ(θo) = 2|θo| π γ . We could have used a non-linear fitting procedure to recover the exponent γ for each reflectance curve; however, we found the following direct fitting solution to be sufficient and much faster. For each (φo, φi) pair, we evaluate ρs at θo = π 3 , apply the linear remapping with respect to ρs at θo = {0, π 2 }, and invert Rγ( π 3 ) to estimate the exponent γ. We chose not to use the exact Fresnel formula in our approach, since the gained accuracy would have been lost with our angular approximation anyway.

The full scratch BRDF ρs with Fresnel effects is then approximated using a simple linear blending of the form:

ρs(ω ω ωo, ω ω ω i ) ≈ Rγ (θo) ρs( π 2 , φo, φ i )+(1-Rγ (θo)) ρs(0, φo, φ i ), (2)
where γ stands for γ(φo, φi). This is illustrated in Figure 3b.

It can easily be shown that this equation ensures energy conservation as well as reciprocity (all involved terms remain the same when replacing θo by θi = -θo). We emphasize that our approximation does not require to simulate the full 3D scratch BRDF; only 3 slices at θo = {0, π 3 , π 2 } need to be computed. In addition, since θo only has to be considered for the evaluation of Fresnel formula, it does not affect ray propagation, which is thus performed in 2D as in Section 4.1. This avoids numerical issues that would otherwise occur when simulating a BRDF slice at θo = π 2 . Finally, since we have Fγ( π2 ) = 1, the scratch BRDF behaves as a mirror at grazing angles, hence ρs( π 2 , φo, φi) = ρm(φo, φi). As a result, our full scratch BRDF consists of a triplet of 2D functions: ρs(0, φo, φi), ρm(φo, φi) and γ(φo, φi), which we store in the 3 channels of a single 2D texture.

For a given view direction ω ω ωo, a slice of the scratch BRDF corresponds to a 1D function on the half-cone of reflection, as illustrated in Figure 3c-bottom (white peaks). Observe in particular the split off-specular peaks due to both profile shape and inter-reflections.

Comparisons Before adding spatial variations to our model, we compare it with two other anisotropic BRDF models in Figure 7 on a simple yet challenging scratch profile: the steep quartic function of Figure 10a. We compute a ground truth image by simulating the full 3D scratch BRDF with Fresnel effects thanks to our dedicated raytracer (Figure 7f) and provide difference images in the bottom right corner for each model. We use a refractive index of n = 2.91 and an extinction coefficient of k = 3.09 (this matches iron at a wavelength of 650nm). Figure 7a uses the Ashikhmin-Shirley model [2000], where we have set nu = 10 4 and nv = 2 in order to focus most of the energy into the half cone of reflection, which is appropriate for scratched materials. Any other single-lobe analytic BRDF model would produce similar results. Figure 7b uses a tabulated NDF as in previous works [START_REF] Yan | Rendering glints on high-resolution normal-mapped specular surfaces[END_REF][START_REF] Dong | Predicting appearance from measured microgeometry of metal surfaces[END_REF]]: we adapt our 2D raytracer to only compute the first bounce off the NDF of a scratch profile, and store the result in our 2D parametrization. In both cases, Fresnel and shadowing/masking effects are separately applied as multiplicative terms. Figure 7c demonstrates that our approach, even when limited to a single bounce, is superior to other models (see difference images), which is due to a more accurate treatment of shadowing and masking. When adding multiple bounces, as seen in Figure 7d, our model becomes nearly indistinguishable from the ground truth. The influence of inter-reflections depends on the choice of parameters for the Fresnel effect, which reduces intensity at each successive bounce. This is clearly seen in Figure 7e, where we use our mirror scratch BRDF, hence giving equal contributions to bounces of different orders. Inter-reflections also have a stronger visual impact when using steep instead of shallow scratch profiles, as shown in Figure 14. In any case, interreflections come for free in our approach.

Spatially-varying BRDF

Our full scratch BRDF may already be used to model uniform materials exhibiting densely packed parallel scratches (with no space in-between). As shown in Figure 1a, this reproduces well the directional smearing of reflections observed in Figure 2-right. However, variations in the spacing, width, height and orientation of scratches are necessary to reproduce glint lines and patterned highlights. Let us first consider a single scratch layer, which in our approach amounts to fix both the orientation and profile of scratches. An important property of our scratch BRDF is that it is invariant to isotropic transformations of the scratch profile. Hence a same layer may hold a distribution of parallel scratches of different sizes, provided they have the same profile and aspect ratio. Any such distribution of scratches may be used and we describe some of our implementation choices in Section 5.3. For now, it is sufficient to consider that scratches in a layer are localized by an indicator function 1S : R 2 → {0, 1}. Scratches of different sizes and spacings lead to different indicator functions 1S, which in turn modify the relative area α of scratches under a given pixel or ray footprint P:

α(x) = P(x)
1S(y) dy.

(3)

While ρs describes the angular properties of our model, α characterizes its spatial variations in a multi-scale manner (see Figure 3ctop). This explicit spatial and angular decoupling is central to the performance of our approach. An implicit assumption is that 1S agrees with the physical width of scratches. Likewise, the footprint P should enclose at least one scratch, which defines the minimum scale at which our SV-BRDF model operates.

Up until now we have considered a single scratch layer. To deal with multiple layers, we first compute a combination of scratch BRDFs weighted by relative area using ρ(x, ω ω ωo, ω ω ωi) = α k (x)ρ s,k (ω ω ωo, ω ω ωi), where we have added a k subscript to denote the k-th layer. This simple linear combination is not guaranteed to conserve energy though, since at locations where scratches overlap, the combined relative area ᾱ(x) = α k (x) might be greater than 1. In that case, we normalize the accumulated BRDF ρ by the sum of weights ᾱ. Otherwise, 1ᾱ represents the ratio of scratch-free surface for which we apply the base material ρ b . Putting it all together, our final SV-BRDF model boils down to:

ρ(x, ω ω ωo, ω ω ωi) = ρ/ ᾱ(x) if ᾱ(x) > 1 ρ + (1 -ᾱ(x))ρ b otherwise. ( 4
)
This equation avoids the computation of the complex BRDF occuring at scratch intersections, replacing it with a simpler average of scratch reflectances. As shown in Figure 3d, our model is sufficient to reproduce convincing glint lines and patterned highlights at multiple scales. Its simplicity is key to efficient evaluation and sampling in complex scenes as explained next.

Rendering with our SV-BRDF model

Our SV-BRDF model is naturally adapted to importance sampling (Section 5.1), but requires special considerations to work with light sampling (Section 5.2). We focus in these two subsections on the radiance reflected by a single scratch. Indeed, thanks to the simple formulation of our SV-BRDF model, the total reflected radiance is obtained by a linear combination of the radiances reflected by individual scratches, following Equation 4. Even if the choices of scratch profiles and distributions are arbitrary, we briefly describe our assets in Section 5.3 for completeness.

Importance sampling

Rendering with our scratch BRDF model amounts to a 1D integral along the half cone of reflections in the local upper hemisphere. We compute the radiance Lr reflected off a single scratch with importance sampling using numerical integration:

Lr(ω ω ωo) = 1 N N i=1 ρs(ω ω ωo, ω ω ωi)cos(φi) PDF φo (φi) Li(ω ω ωi), (5)
with N the number of samples and Li the incoming radiance.

The sample directions ω ω ωi are randomly generated on the reflection cone (i.e., where θi = -θo) following a 1D probability distribution function based on the mirror scratch BRDF: PDF φo (φi) = ρm(φo, φi)cos(φi). To this end, we pre-compute the 1D inverse cumulative distribution functions (CDF) for each φo-slice of the texture storing ρm, which is both fast and compact. We could have instead pre-computed and stored the inverse CDF for the full scratch BRDF ρs, but this would have been both time and memory consuming. Our approach still remains more accurate than relying on a NDF as done in previous work, since we are taking interreflections into account. Note that the cosine term in Equation 5cancels out with this choice of PDF. Moreover, when Fresnel effects are discarded, only the Li(ω ω ωi) term remains. When rendering with multiple scratch layers, we distribute N samples based on layer importance, which is simply given by α k in our model.

Light sampling

When light samples are provided by the renderer, as with (virtual) point light sources for instance, we are faced with a problem: since sample light directions are infinitely small, they will never exactly lie on the half cone of reflection. Following [START_REF] Yan | Rendering glints on high-resolution normal-mapped specular surfaces[END_REF], we assign a narrow reconstruction kernel function around a sample light direction ω ω ω ℓ to prevent this from happening. Any isotropic kernel would do; in practice we use a function fΨ : [-Ψ, Ψ] → [0, 1] with a compact support of Ψ equivalent to one degree as shown in figure 8. We then locate the direction ω ω ω ′ ℓ on the mirror cone that is closest to ω ω ω ℓ : it is given by (-θo, φ ℓ ) in our angular coordinates. We next compute the intersections of the half cone of reflection with the support of fΨ (see Figure 8):

φ± = φ ℓ ± Ψ 2 -(θ ℓ -(-θo)) 2 .
The radiance Lr reflected by a single scratch for a single light sample of intensity L ℓ is again computed by numerical integration:

Lr(ω ω ωo) = (φ+ -φ-)L ℓ πN N i=1
ρs(ω ω ωo, ω ω ωi) cos(φi)fΨ(ψ i,ℓ ), (6)

with N the number of samples, ω ω ωi a direction on the cone segment given by (-θo, φi) in our angular coordinates (with φi uniformly distributed in the [φ-, φ+] range), and ψ i,ℓ the angle between ω ω ωi and ω ω ω ℓ . We could alternatively use importance sampling of the BRDF on the [φ-, φ+] range, but we found it not to be necessary considering the small kernel support we use.

When dealing with multiple scratch layers, many combinations of light and scratch directions will make no contribution to reflected radiance. This happens when the half cone of reflections does not intersect the support of fΨ, or equivalently when Ψ < |θ ℓ + θo|.

It is thus straightforward to prune sample light directions and/or scratch layers prior to rendering by simply checking the latter inequality. This is a special case of the method of [START_REF] Yan | Rendering glints on high-resolution normal-mapped specular surfaces[END_REF].

Our solution is simpler only because it is tailored to the geometric properties of scratches, which shows the advantage of having a dedicated scratched material model.

Implementation details

Scratch profiles We provide users with two different controls over profile shapes, as shown in Figure 5: analytical functions and tabulated height fields from Merillou et al. [2001]. In the latter case, the profile might exhibit bulges that stick out of the base surface, which violates our hypothesis that scratches are only dug inside. We correct this situation by taking the highest profile point as the base surface height, and making the profile periodic during 2D simulation (ray tracing). This solution works well in practice, especially at grazing angles where profile bulges are naturally masked by nearby, repeated bulges.

Scratch distributions

Each layer of a distribution may have its own parametrization on the surface, which may be arbitrarily chosen by the user. In practice, we always rely on a single surface tangent field, which is by default aligned with the u coordinate in texture space. Curved scratches are thus easily achieved by modifying the parametrization as demonstrated in Figure 15. The scratch indicator function 1S k for a layer k is built by binarizing a noise function stretched along u (see figure 9):

1S k = N R(Θ k ) l 0 0 w -1 u v + ξ ξ ξ k d ( 7
)
where N is the Worley's cellular 2D noise function [START_REF] Worley | A cellular texture basis function[END_REF]], and ⌊x⌋ d is a thresholding function equal to 1 if x ≤ d and 0 otherwise. The density parameter d determines the proportion of scratchversus-base material. Scratch dimensions are controlled by the mean scratch length l along u and the mean scratch width w along v. Each scratch layer is rotated by an angle Θ k randomly picked in the user-specified [-Θ 2 , Θ 2 ] range. The offset vector ξ ξ ξ k ∈ [-1, 1] 2 randomly shifts layers with respect to each other. Table 1 collects distribution parameters for all our results, along with the number N of layers. Table 1: Distribution parameters used in all our results.

Image N Θ l w d 1 (b) 4 7 0.2 2 • 10 -4 0.6 1 (c) 16 360 0.1 1 • 10 -5 0.3 1 (d) 32
The computation of the multi-scale relative area α (Equation 3) varies depending on the type of implementation. In a GPU-based renderer, the simplest approach is to build a mipmap texture out of 1S and rely on graphics hardware to deal with aliasing issues for distant views. In an off-line renderer, the method will depend on the available texture filtering functionalities, which is most often transparent to the user. When working with multiple layers, we rotate and offset the same noise patterns or use different noise functions depending on the desired regularity of the distribution.

Results

Interactive material editing

We have implemented our SV-BRDF model in a custom OpenGL progressive rendering engine. Even though it only works with direct lighting, it provides a sufficiently rich and fast visual feedback for material design. Users may modify the scratch profile, Fresnel parameters, scratch distributions, camera or lighting interactively as shown in Figure 10 and the accompanying video. Performance range from real-time to a few seconds to reach convergence, depending on the number of scratch layers and the type of lighting.

We provide timings obtained with a Nvidia GTX680. Our 2D raytracer is implemented in CUDA to provide fast updates of the scratch BRDF (along with the PDF and inverse CDF) when the profile is modified. Preprocessing takes between 100ms and 300ms for filling a BRDF of resolution 256 × 256, sending N = 10K samples per φo slice. Our progressive renderer works in GLSL by accumulating the results of multiple passes, and works with both light sources and environment illumination. With a naive approach, performance scales linearly with the number of scratch layers. We have implemented a practical culling method that discards the k-th layer whenever α k < 10 -6 , with no perceivable difference in results. We obtain significantly better (sub-linear) timings with this optimization, as summarized in Table 2. Our brushed patterns usually require no more than 6 scratch layers in practice. To achieve fast feedback as demonstrated in the supplemental video, we distribute samples across multiple frames as is typically done in progressive rendering. Compared to the method of [START_REF] Yan | Rendering glints on high-resolution normal-mapped specular surfaces[END_REF] that requires to recompute and sample a NDF for each pixel, our solution is inherently much faster. The reason is that our BRDF amounts to a few texture lookups and simple analytic functions, and supports pre-filtering through the α-maps. We have also tested importance sampling on the scratch layers by redistributing the total number of samples based on α k . Figure 11 shows that with only 128 samples for all of the 32 scratch layers already produces clean results in practice compared to uniformly distributing the samples.

The memory footprint of our method is very small: we typically use a 256 × 256 three-channels texture for the scratch BRDF, a two-channels texture for the precomputed inverse CDF and a onechannel texture for the PDF (1.5 MB without compression). The relative area α(x) is currently stored in a texture, even though it could be computed procedurally if memory usage becomes critical. Note however that the relative area map is usually much smaller than would be a normal map with the same accuracy. More specifically, if the width of the biggest scratch in a high-resolution normal map is w (in texels), then a relative area map with a resolution w times smaller in each dimension will be sufficient. In practice, we observe that a minimum of w = 10 texels is necessary to capture a scratch profile, hence our relative area map is at least 100 times smaller than would be a normal map of similar quality.

Global illumination results

We have implemented our SV-BRDF model as a Blender plugin for the Cycle renderer. This not only permits to demonstrate scratched materials with global illumination, but also allows us to leverage the shading and texturing network functionalities of Cycles to quickly create interesting scratch distributions, as demonstrated in the supplemental video. Figures 1, 12 and 14 include insets that are actual 3D (camera) zooms, not image magnifications.

Figure 1a makes use of our scratch BRDF; it corresponds to a densely packed set of parallel scratches. Even in this simple setting, our model already reproduces the characteristic effects of a brushed metal, such as the directionally-smeared reflections observed in Figure 2 right. The main highlight also shows slight variadue to split off-specular peaks in the BRDF. Figure 1b instead uses our full SV-BRDF with a scratch distribution made of 4 layers oriented around a single main direction. This reproduces the appearance of brushed metal, as in Figure 2 middle: note in particular the patterned highlights as well as the mirror-like reflection due to the base BRDF showing up in-between scratches. Figures 1c-d instead use distributions where layers span all orientations, with two different densities of scratches (with 16 and 32 layers respectively) and deep versus thin profiles. This result is similar to Figure 2 left, where glint lines appear distinctly, especially around strong highlights. Notice also the overall hazy environment reflection obtained when increasing the density of scratches. Inset figures reveal the complex multi-scale behavior of our model: scratches not only get more resolved, but also let different reflection patterns emerge.

Our SV-BRDF is adapted to the modeling of a variety of scratched materials. Figure 12a reproduces the appearance of scratched gold, where scratches are so tiny that they are barely noticeable. Here layers are oriented around a main direction on the surface, which produces an elongated highlight blended with a mirror reflection. Such highlights are also visible in macro-scale inter-reflections, taking an orange tint, as is best seen in the supplemental video. Figure 12b uses two different BRDFs for the scratches and the base surface, imitating the appearance of a scratched painted metal. A sparse set of multi-directional layers reproduces the random deteriorations expected with such materials. Figure 12c makes use of a denser distribution of scratches in all directions, and adds a diffuse term common to the base and scratch BRDFs. Together with an appropriate choice of refractive index (n = 1.5), this reproduces the appearance of scratched plastic. Figure 12d uses a similar setting to imitate the appearance of finished wood, by simply using a textured diffuse term in place of a uniform color. Here again, inset figures reveal the behavior of our SV-BRDF model in close-up views, where different reflection patterns emerge compared to zoomed out views.

We have evaluated the performance of our plugin by comparing the anisotropic GGX model of Cycles to our scratch BRDF model (Equation 2). This amounts to consider a single layer of scratches that covers the entire object. Figure 13 shows renderings at 1000 samples per pixel with either BRDF model. Our approach produces a brighter result because of inter-reflections. Compared to the GGX model that takes 484 seconds, our method is slightly slower, taking 550 seconds. However, it should be considered that our implementation remains unoptimized.

We have found that the choice of scratch profiles has a significant impact on material appearance, as demonstrated in Figure 14. Nearly flat profiles (left column) tend to produce distinct offspecular peaks in the BRDF, which show up as multiple highlights on the surface both in densely brushed or sparsely scratched distributions. Steeper profiles (right column) lead to a different material appearance. The number of inter-reflections is increased, hence

Limitations and future work

Similar to previous work (e.g., [START_REF] Dupuy | Linear efficient antialiased displacement and reflectance mapping[END_REF]Jakob et al. 2014b;[START_REF] Dong | Predicting appearance from measured microgeometry of metal surfaces[END_REF]), our model makes some simplifications to enable efficient rendering. In particular, we consider idealized and repeated scratches, devoid of imperfections along the extrusion direction. Imperfections might be added by combining techniques dedicated to punctual glitter effects (e.g. [Jakob et al. 2014b]) with our solution. Pitted materials [START_REF] Pont | Bidirectional reflectance distribution function of specular surfaces with hemispherical pits[END_REF] could also be modeled and rendered at multiple scales with an approach similar to ours. Another way to increase visual richness would be to let all scratch parameters (profile, width, height, etc) vary per layer. This would require to devise an analytic scratch BRDF model that could be instanced on the fly, as opposed to running a simulation per scratch profile in preprocess.

In this work we have essentially focused on the SV-BRDF of the outer scratched interface of a material, which has two limitations. First, our spatial blending scheme becomes apparent for extreme close-up views where a scratch is significantly wider than a pixel (see Figure 14-left). This limitation is similar to that of [START_REF] Irawan | Specular reflection from woven cloth[END_REF], whose method also exhibits simplified fiber patterns in close-up views. We plan to extend our SV-BRDF model to a BTF model for extreme close-ups in future work, taking advantage of a scratch structure to avoid storing a full BTF. Second, we have not considered the complex SV-BTDF necessary to model scratched transparent materials, such as glass for instance. This is more challenging than reflection due to Snell's law that directly depends on the refractive index, hence making precomputations more difficult.

Ultimately, the decisions made in our work as well as in previous work on scratched and glittery materials should be compared to measured materials. This is a challenge in itself due to the high angular frequency of such BRDFs. Existing work either reduce the angular resolution to capture a full BRDF [START_REF] Filip | Effective acquisition of dense anisotropic brdf[END_REF] or measure a BRDF for a subset of directions only [START_REF] Dong | Predicting appearance from measured microgeometry of metal surfaces[END_REF]]. In the latter case, the authors also measured micro-geometries using an optical profilometer to directly fit their NDF model to brushed metal surfaces. Even though this is limited to single scattering effects (i.e., thin scratch profiles), we believe a similar approach could be used with our model. In our case, the challenge would be to convert the measured micro-geometry to a set of layered scratch primitives, a direction we plan to pursue in future work. We believe that 3D ray-tracing simulations could then be used to compare measured versus fitted micro-surfaces as a validation procedure.

Conclusion

We have shown that a dedicated SV-BRDF model permits to reproduce the high-frequency spatial and angular variations of scratched materials, without requiring heavy computational or memory resources. This has direct practical implications: scratched materials may be designed at interactive rates with fast feedback, and rendered in global illumination with a realistic final look. In particular, our model reproduces the visual appearance of scratched materials, such as the glint lines, patterned highlights and directionallysmeared reflections identified in Figure 2. This is made possible by exploiting the regular micro-structure of scratch distributions, finding a practical balance between precomputations (simulating scratch BRDFs with inter-reflections) and procedural distributions (linearly combining user-defined scratch layers). A benefit of our approach is that it is easily integrated in existing rendering packages and hence extensible to new scratch profiles and distributions. We believe that similar structured SV-BRDF models for other categories of anisotropic materials such as fabrics could also be devised.

Figure 2 :

 2 Figure 2: Photographs of scratched materials (top row) used as inspiration for some of our results (bottom row). Scratches appear as glint lines (left) at moderate distances, while at farther distances, they produce patterned highlight silhouettes (middle) as well as directionally-smeared environment reflections (right).

Figure 3 :

 3 Figure 3: (a) The BRDF for a single mirror scratch is simulated using 2D ray tracing (each ray is given a consistent color across reflections for illustration purpose). (b) The full 3D scratch BRDF with Fresnel effects is approximated using an adequate interpolation between a pair of 2D BRDF slices, controlled by an exponent map γ. (c) At render time, a SV-BRDF is constructed by linear combinations of scratch BRDFs (bottom) under a pixel or ray footprint based on the relative area of scratches (top). (d) Our SV-BRDF model is efficiently evaluated at multiple scales, exhibiting both high spatial and angular frequencies.

Figure 4 :

 4 Figure 4: (a) For a mirror scratch aligned with u, light reflected in the outgoing direction ω ω ωo lies on a half-cone of directions (in purple). (b) Expressing any direction by a pair of scratch-dependent angular coordinates (θ, φ), the half-cone becomes a θ-isocurve. (c) A mirror scratch BRDF is reduced to a 2D function of φ angles.

Figure 5 :Figure 6 :

 56 Figure5: Scratch profiles (shown on top) and their corresponding mirror scratch BRDFs (1st, 2nd and higher-order bounces in red, green and blue respectively). Top row: a pair of symmetric quartic profiles yield doubly symmetric BRDFs, with the steepest profile (right) yielding more higher-order bounces. Botton row: a pair of measured profiles yields more complex BRDFs with only one symmetry (due to reciprocity). Even though the two profiles on the right column have a similar global shape, their BRDFs differ significantly. In particular, the peak on the anti-diagonal (i.e., mirror reflection) is absent on the bottom right BRDF due to the lack of flat profile portions.

Figure 7 :

 7 Figure 7: Comparison of anisotropic BRDF models; the top-left insets show a BRDF slice at θo = 0 in our parametrization, the bottom right insets show the difference with ground truth (boosted 8 times). (a) Rendering with the Ashikmin-Shirley BRDF exhibits uniformly smeared reflections. (b) Using a data-driven NDF produces more realistic non-uniform smearing. (c) Our scratch BRDF model without inter-reflections is similar to (b), with more accurate masking and shadowing. (d) Adding inter-reflections produces a result nearly identical to the ground truth in (f). (e) Inter-reflections become strengthened when Fresnel effects are discarded.

Figure 8 :

 8 Figure 8: Geometric intersection between the support of a sample light direction ω ω ω ℓ (in orange) and the half cone of reflection (in purple). The intersected segment (solid purple curve, top right) is identified by its start/end angles φ±. For the falloff function (bottom right), we use fΨ(ψ) = (1 -(ψ/Ψ) 2 ) 2 .

Figure 9 :

 9 Figure 9: A 2D cellular noise is streched and rotated, then binarized to obtain the scratch indicator function. Here we show one layer with Θ = 5, l = 0.25, w = 2 • 10 -2 and d = 0.5.

Figure 10 :

 10 Figure10: Screenshots from an interactive material editing session, at a resolution of 615 × 780 using 64 samples per pixel and per scratch layer. Progressive rendering is done with direct lighting without shadows to obtain fast feedback. We start in (a) with the BRDF of a quartic scratch profile applied in (b) with a brush pattern (using 6 layers) on our test object lit by an environment map and a red light source. We then switch in (c) to a measured profile[Mérillou et al. 2001], recomputing the BRDF in 100ms. This changes the material appearance in (d): scratches are less pronounced and reflections less elongated. In (e) we rotate the red directional light source, which helps probe the response of the material to different light orientations interactively, taking 749ms. We swicth to a scratch pattern (using 32 layers) in (f), leaving other material parameters unchanged: despite the increase in the number of layers, rendering takes only 776ms thanks to layer culling (see text).

Figure 11 :

 11 Figure 11: Rendering with only 128 samples per pixel (for all of the 32 scratch layers) produces noisy scratches. If instead we use importance sampling on scratch layers with the same number of samples (right), scratches are much better resolved.

Figure 12 :

 12 Figure 12: Four materials rendered with global illumination, using different profile, Fresnel and distribution parameters. Top row: Separate contributions of base and scratch materials. Bottom row: final appearance combining both components. (a) Gold (4 layers) is obtained by a regular distribution of sparse tiny scratches that are combined with the base, mirror-like reflection. (b) Painted metal (16 layers) makes use of a sparse distribution of scratches in multiple orientations, and a paint BRDF for the base material. (c) Plastic (32 layers) uses Fresnel formula for dielectrics and adds a uniformly colored diffuse term to the dense multi-directional scratch layers. (d) Finished wood (16 layers) is similar to Plastic, except that it uses a sparser set of scratches and a textured diffuse term.

Figure 13 :

 13 Figure 13: Global illumination results used for performance comparison. Left: the GGX BRDF of Cycles. Right: our scratch BRDF.

Figure 14 :

 14 Figure14: The choice of scratch profile has a strong impact on material appearance, in either brushed (top row) or scratched (bottom row) patterns. In the first two columns, we use nearly flat profiles, which exhibit multiple highlights due to off-specular peaks in the scratch BRDF. The left-most column shows the limitation of our method for extreme close-up views focused on the area pointed by the red crosses: individual scratches span a few pixels. The third column uses a steeper scratch profile, which yields more constrasted reflections due to multiple bounces being attenuated by Fresnel effects. We compare our results to a NDF-based model in the last column, which shows the importance of interreflections. The material appears a lot dimmer when using a NDF since many scratches reflect little or none of the incoming radiance when higher-order bounces are discarded.

Figure 15 :

 15 Figure 15: Curved scratches are obtained by modifying the underlying parametrizations. Here we have used the same profile and distribution as in the middle two columns of Figure 14 (top row).

Table 2 :

 2 Performance

	Nb layers	1	6	16	32
	Brushed	166	749	1850	3540
	Scratched	93	185	419	776

(in ms)

for rendering a 615×780 image with 64 samples per pixel and per scratch layer, which provides good visual quality. We use culling based on α as explained in the text and render materials using both brushed and scratched patterns with a varying number of layers. Framed timings correspond to Figure

10a,d,e

(for brushed) and Figure

10f

(for scratched).

Acknowledgements

This work has been supported by the ANR RichShape project (ANR-14-CE24-0004), and by the EU Marie Curie Initial Training Network PRISM (FP7 -PEOPLE-2012-ITN, Grant Agreement: 316746). The environment maps are provided courtesy of the sIBL Archive for Figures 1 and7, and of the USC Institute for Creative Technologies for Figure 10. We thank the reviewers, Pierre Bénard and Romain Pacanowski, for insightful comments on the paper.