Vincent Mouchi
email: vmouchi@gmail.com

Quentin G Crowley

Teresa Ubide

AERYN: a simple standalone application for visualizing and enhancing elemental maps

Keywords: LA-ICP-MS, EDS, WDS, element distribution maps, processing, Matlab, biominerals

Interpretation of high spatial resolution elemental mineral maps can be hindered by high frequency fluctuations, as well as by strong naturally-occurring or analytically-induced variations. We have developed a new standalone program named AERYN (Aspect Enhancement by Removing Yielded Noise) to produce more reliable element distribution maps from previously reduced geochemical data. The program is Matlab-based, designed with a graphic user interface and is capable of rapidly generating elemental maps from data acquired by a range of analytical techniques. A visual interface aids selection of appropriate outlier rejection and drift-correction parameters, thereby facilitating recognition of subtle elemental fluctuations which may otherwise be obscured. Examples of use are provided for quantitative trace element maps acquired using both laser ablation (LA-) ICP-MS and electron probe microanalysis (EPMA) of the cold-water coral Lophelia pertusa. We demonstrate how AERYN allows recognition of high frequency elemental fluctuations, including those which occur perpendicular to the maximum concentration gradient. Such data treatment compliments commonly used processing methods to provide greater flexibility and control in producing elemental maps from micro-analytical techniques.

Introduction

Element distribution maps or images acquired using in-situ analytical techniques are widely used in Geology and Earth Sciences (e.g. Cusack et al., 2008; Rittner & Müller, 2012; McGowan et al., 2014).

Several software packages exist to process element distribution maps [START_REF] Paul | CellSpace: A module for creating spatially registered laser ablation images within the Iolite freeware environment[END_REF][START_REF] Rittner | 2D mapping of LA-ICPMS trace element distributions using R[END_REF]; [START_REF] Ortolano | X-Ray Map Analyzer: A new ArcGIS based tool for the quantitative statistical data handling of X-ray maps (Geo-and material-science applications)[END_REF], with some recent studies having investigated multi-trace element maps processed from data obtained by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) (e.g. [START_REF] Petrus | A novel 2D LA-ICP-MS data analysis and visualization solution[END_REF][START_REF] Paul | Towards a method for quantitative LA-ICP-MS imaging of multi-phase assemblages: mineral identification and analysis correction procedures[END_REF][START_REF] Ubide | High-resolution LA-ICP-MS trace element mapping of igneous minerals: in search of magma histories[END_REF].

Common issues in micro-analytical techniques linked to data acquisition, such as analytical drift and outliers, can greatly influence the quality of resultant elemental maps. We propose a user-friendly method created using the programming software Matlab (MathWorks; http://www.mathworks.com/products/matlab) to deconvolute complex data and produce more readable elemental maps. LA-ICP-MS data must have been previously reduced to account for baseline correction and reference material calibration, e.g. by using common data reduction software. Our program, which we term AERYN (Aspect Enhancement by Removing Yielded Noise), allows the filtering of outliers and smoothing of values included in a data matrix, resulting in an increased readability and importantly permitting more representative and accurate interpretation of heterogeneity in elemental distribution. In addition, we propose a de-trending option to highlight subtle geochemical features which may be obscured by strong signal fluctuations, either natural in origin or caused by instrumental drift. AERYN has been compiled and exported as a standalone application for Windows, allowing its use on a PC without the requirement of a Matlab licence (http://www.tcd.ie/Geology/staff/crowleyq/AERYN/). The only requirement for this program to perform correctly is to download and install the free Matlab Runtime R2015a (8.5) from the MathWorks website (http://www.mathworks.com/products/compiler/mcr).

To illustrate how AERYN operates, we present elemental maps acquired by LA-ICP-MS and electron probe microanalysis (EPMA) from the wall of the skeleton of the cold-water coral Lophelia pertusa.

Virtually any type of numeric compositional map can use part or all of the protocol described here, as long as the data fit to a square or rectangular shape (in laser-based systems, this can be achieved by analysing either parallel equivalent line scans or a grid of regularly-spaced spots). We provide the corresponding Matlab code and the standalone application as well as a "Read Me" file with some specifications as supplementary material, allowing the software to be adapted by other users.

AERYN description

User interface

AERYN is presented with a control panel on the left-hand side and three graphic areas on the righthand side (Fig. 1). The upper graphic displays an elemental map (using Matlab's "imagesc" function) whereas the lower-left graphic displays a 3-dimensional view of the data (using Matlab's "surf" function). The lower-right graph represents a histogram of the data (number of bins calculated using Yule's rule) to graphically illustrate the distribution of elemental concentration represented in the other graphic displays.

Importing data

AERYN uses text files (.txt, .tsv) exported from a variety of elemental map processing software or from Microsoft Excel. We tested exports from LA-ICP-MS, Scanning Electron Microscopy with electron-dispersive spectroscopy (SEM-EDS) and EPMA maps, but text files from any software should be compatible. If not, we advise the user to follow instructions provided in the "Read Me" file.

AERYN does not process raw data, and a preceding data reduction step is still required. We used Iolite [START_REF] Paton | Iolite: Freeware for the visualisation and processing of mass spectrometric data[END_REF]v. 2.5 and 3.0) for data reduction prior to importing in AERYN. To import the numeric matrix corresponding to the map file, AERYN passes any number of lines composed of characters and reads only the part of the file containing the concentration values. The size of the imported matrix is only restricted by computer performance for very large datasets.

Once a numeric matrix has been generated, it is automatically displayed as a map on the upper graph and as a 3D surface plot of the 2D data on the lower-left graph (Fig. 1). It is therefore possible to rotate the lower-left graph or to automatically orient the view to a side section, in order to have better visibility of outliers and fluctuations. Moreover, if required for better display, it is possible to transpose the matrix by clicking on "Transpose" in order to invert lines with columns.

From this point and at any stage, the user can cancel any number of operations made on the data by clicking on the "Undo" button.

Outlier removal

Outliers can be either extremely high or low values in comparison to the rest of the observed data and introduce bias in interpretations. AERYN allows the user to locate and remove all outliers at once. All removed outliers are replaced by empty spaces in the matrix (NaN). Processing of these spaces is explained in the following section. A statistical test is proposed to the user to define automatic thresholds of outliers based on the following equations [START_REF] Tukey | Exploratory Data Analysis[END_REF][START_REF] Carling | Resistant outlier rules and the non-Gaussian case[END_REF]:

(1) Tup = P90 + 1.5 * IQR

(2) Tdown = P10 -1.5 * IQR With P90 the percentile at 0.90 for the dataset, P10 the percentile at 0.10 and IQR = P90 -P10. All values not occurring within the interval between Tdown and Tup are considered as dataset outliers.

Such a statistical treatment will likely be sufficient for most datasets, though for the data presented here the cut is too drastic and the upper end of some high value areas is lost, preventing an accurate reading of the corresponding concentrations. For this reason, a second technique, based on iterative visual definition of outliers, is made available for the user.

Firstly, for the high value outliers, it is necessary to estimate the upper limit of data to be retained. For this, the "surf" function of Matlab is extremely powerful. By displaying the elemental map in 3D on the lower-left graph (with the concentrations in the z axis) and with the possibility of rotating the image, one can clearly (i) identify variations in concentration, (ii) assess the potential presence of trends that would require a correction at a later stage, and (iii) identify obvious outliers and their approximate values. At this stage, successive outlier removal steps can be achieved by searching for values in the dataset lower or higher than a chosen specific value (Fig. 2). Moreover, the histogram on the third graph (lower-right area) allows a better estimate of the minimum and maximum values to be chosen for the colour palette in order to highlight most of the compositional information.

Spatially-restricted high concentration areas are maintained as a whole with this method, though more subjectively compared to the automatic definition of outliers as presented above. It should be noted that a possibility of inadvertently "over-trimming" the dataset exists, especially when controlled by the user, and therefore outlier removal should be done with caution (e.g., Fig. 2c).

For the second phase of outlier removal, a lower minimum limit value expected for quantitative or semi-quantitative data is zero. This minimum value can be higher if the concentrations of the studied elements are high, but in the present example the concentrations are low and can easily be cut to positive values without further consideration. If the mapped area is not entirely filled with sample material (i.e., if a sample was mounted in epoxy resin and this resin was intentionally ablated in order to obtain a rectangular map), the user must ensure that the resin value, which is generally constant, is still included in the data to prevent incorrect resampling of these areas by the outlier removal step.

This step can have multiple iterations as successive outlier removals allow an improved visualisation of geochemical fluctuations. The automatic Side View is extremely valuable at this stage as the outlier values as well as the signal to be retained are indicated by the vertical axis. During the whole process however it is necessary to keep in mind that some valuable geochemical information can reside in high-frequency variability. It is thus unwise to crop the highest or lowest values too close to the main visible fluctuations (see Fig. 2a andcaption). Moreover, one must be aware that the larger the holes caused by outlier removal, the less accurate the following steps will become. We suggest that the number of outliers compared to the size of the data must be reasonably low. As an indication, the number of points removed in Figure 2 represent approximately 0.5 % of the dataset (see Fig. 2c).

Each time the user filters the dataset, a temporary message is displayed in the history box to indicate the percentage of data cropped compared from the entire matrix. The cut values are displayed on the map as white areas. Once outlier removal has been validated, the side view and the histogram can provide preliminary boundaries for the colour palette (Fig. 2b). The filters can be cancelled to retain the entire initial dataset while maintaining the palette range values.

Filling gaps following outlier removal

In the outlier removal step, the outliers have been automatically replaced by the specific passive value NaN. The "inpaint_nans" function [START_REF] D'errico | Inpaint_nans. MATLAB Central File Exchange[END_REF] allows replacement of NaNs by analyzing surrounding values for interpolation and extrapolation. With this step completed, the maps are efficiently "cleaned" of outliers.

This gap-filling step should not be performed after further steps such as smoothing, as the method used here substantially increases the number of NaN values in direct connection to pre-existing NaN values. For this reason, we chose to automatically fill the NaNs at the validation of the outlier removal.

Smoothing

An option to perform median filtering is available at this stage of data correction. Numerous smoothing functions and methods are available in the Matlab environment. In order to take into account the possible variability present on both x and y axes, the "medfilt2" function was chosen. This function calculates the mean of all values included in a moving square centred on a value to smooth, in order to replace that value with the mean, before moving the mask to the adjacent value. A larger mask results in a stronger applied smoothing function. If the dataset presents a strong background noise, it may be valid to repeat this step a second time, but it should be noted that if this option is performed, absolute values of concentration will not be accurate. The mean will be identical, but the absolute values of the extremes will no longer correspond to those given by initial data reduction and the map will therefore be semiquantitative. Relative fluctuations, however, will remain the same.

Detrending

The detrending option is useful in two different scenarios. Firstly, if a strong naturally-occurring trend is measured in the sample, the colour palette (commonly used to present fluctuation patterns) will only show this particular trend. Any higher frequency fluctuations along this trend or variations orthogonal to this direction (commonly observed in biominerals due to non-linear growth) and/or fluctuations of lower amplitude are thus invisible. Removing this main trend will allow qualitative study of such fluctuations. Additionally, if the background noise or sensitivity gradient is skewed as a result of experimental drift, a correction must be applied to stabilize the response across the elemental image. The technique presented here makes it possible to use data, even when a dramatic instrumental drift occurs during an analytical session.

As is the case for smoothing, the process of detrending may change absolute values, so any interpretations can only be proposed for elemental distribution and relative fluctuations in concentration, in a semiquantitative fashion. Even in the case of naturally occurring trends, this option allows observation of lateral variations in concentrations that would otherwise be concealed.

The particular data set presented here shows a general trend of a supposedly decrease in Sr all along acquisition time, with data acquisition progressing in successive lines from left to right then top to bottom.

A first option to deal with this trend is to execute an automatic detrending. By clicking the corresponding button, AERYN calculates the best fitting linear plane to the data and applies a detrending corresponding to this plane.

Alternatively, to improve visibility in this case we can select a transect representative of the main fluctuation. By calculating the equation of the linear trend of this transect, it is possible to apply this equation as a correction to all lines (or columns) and reveal geochemical fluctuations across the entire map (Fig. 3). With this simple operation and even in the case of instrumental drift, it is possible to work on the map without necessarily having to return to data acquisition. When selecting a line or a column and clicking on "Display", two new windows present the selected transect and its calculated linear trend on one panel (Fig. 3a), and the same transect along with its detrended transect on the other (Fig. 3b). No operations to the map have been done at this stage. If satisfied with these detrending parameters, the user can close both computer windows and Validate. This will detrend the matrix and display the detrended map on all graphs (Fig. 3c). If the calculated transects are not satisfactory, they can be cancelled by resetting the detrending operation.

The corrected map now presents fluctuations in Sr which are visible as coherent bands both at the top and bottom of the map, showing continuity of these features. The map can be smoothed to clean the high-frequency fluctuations using a 5x5 mask (Fig. 4).

Export map

At any stage, it is possible to export the matrix. In order to maintain optimum picture quality, the maps are exported as TIFF files without compression with a uniform size of both dimensions of the map. The user is asked to enter the path and file names. The matrix itself and the historic record of operations made are saved as tab delimited text files, and the data are subsequently importable in AERYN for further modification if required.

Figure 5 presents the Sr map used as an example here as well as a Li map from the same acquisition session after data reduction (Iolite v. 3.0) and after post-processing using AERYN. Whereas maps processed only with data reduction software are noisy and hamper visibility of concentration bands, these bands are clearly visible after processing with AERYN. If required, in order to counteract the standardization of the dimension of the map after saving with AERYN, it is possible to import the TIFF file in another image software (e.g., Adobe Photoshop) to resize the image to the proper aspect ratio.

Example materials

Specimens of the cold water coral L. pertusa were prepared at Trinity College Dublin using the protocols described in [START_REF] Mouchi | Potential seasonal calibration for palaeoenvironmental reconstruction using skeletal microstructures and strontium measurements from the cold-water coral Lophelia pertusa[END_REF]. Cleaned aragonitic skeletons were embedded in epoxy resin and cut longitudinally with a Buehler Isomet Low Speed Saw. Sections obtained were polished and subsequently etched with 2 % formic acid for 50 s. The etching protocol allowed direct acquisition of measurements on LA-ICP-MS without pre-ablation.

LA-ICP-MS

Geochemical data were collected using an Analyte Excite 193 nm excimer ArF* laser ablation system (Photon-Machines Inc.) coupled to an iCAP-Q ICP-MS unit (Thermo Scientific) at the Trinity College Dublin geochemistry facility. Details on the instrument setup for elemental mapping can be found in [START_REF] Ubide | High-resolution LA-ICP-MS trace element mapping of igneous minerals: in search of magma histories[END_REF]. Parallel and slightly overlapping (1 µm) line scans were acquired as suggested by [START_REF] Ubide | High-resolution LA-ICP-MS trace element mapping of igneous minerals: in search of magma histories[END_REF] to avoid gaps during data acquisition without notable contamination from adjacent lines. A square beam of 12 µm was selected with 10 Hz repetition rate at 8 µm.s -1 scan speed and a laser fluence of 3.9 J.cm -2 . Ion intensities of 7 Li, 43 Ca, 86 Sr and other isotopes not discussed here were measured for a total dwell time of 146 milliseconds. Individual dwell times were 70 ms for 7 Li, 5 ms for 86 Sr and 10 ms for 43 Ca. Data acquisition procedures ran for a total of 2.5 h.

Either NIST SRM 610 or NIST SRM 612 was ablated for approximately 120 s using the same ablation parameters as the unknown line scans. Two calibration lines were used at the beginning and end of each experiment and one additional line every 10 unknown sample lines which equates to 120 s of data acquisition on the reference material every approximately 25 minutes. Collected data were then processed for quantitative values with the Iolite software [START_REF] Paton | Iolite: Freeware for the visualisation and processing of mass spectrometric data[END_REF]v. 2.5 and 3.0) using the Trace_Elements data reduction scheme (included in the Iolite program package) and the recently accepted values of NIST SRM 610 and NIST SRM 612 [START_REF] Jochum | Determination of reference values for NIST SRM 610-617 glasses following ISO guidelines[END_REF]. 43 Ca was used as the internal standard, assuming a stoichiometric Ca content. Reduced data were then exported to be readable and further treated by AERYN; instructions to export text files readable by AERYN from Iolite are in the "Read Me" file provided as supplementary material.

Depending on the size of the acquisition area for each specimen, typical maps have between 950 and 1950 columns for 30 to 60 lines, corresponding to approximately 1-2 mm by 330-660 µm, respectively.

The description of AERYN presented above (Fig. 12345) used the map of distribution of Sr in L. pertusa as example. The map of distribution of Li (Fig. 6a) does not require the outlier removal step as the Side View indicates that the majority of measured fluctuations are below the resin value (see the arrows on Fig. 6a) and visible. Indeed, though some positive outliers are present on the right-hand side, they do not prevent reading of the minimum and maximum values on the histogram for palette definition. The map is however difficult to interpret in terms of the location of fluctuations in Li concentrations. The top part of the map is depleted in Li, as indicated by blue areas in opposition to red and yellow areas on the lower part of the map. An automatic detrending applied to reduce the amplitude of the fluctuations from top to bottom (Fig. 6b) aids in the visualization of several areas of higher and lower Li concentration oriented N-S in this view.

From this step, a light smoothing using a 3 by 3 points mask further highlights the fluctuations. A wider mask was not chosen due to the reduced number of lines in this map (only 30). Except for the high value areas (in red) on the rightmost part of the map, both the map and the Side View show higher amplitudes of fluctuations on the left-hand part compared to the right-hand part of the map (Fig. 7). As this side of the coral skeleton displaying elevated elemental concentrations corresponds to the initial stage of skeletal growth, this map can be interpreted in terms of environmental forcing or kinetic incorporation, as suggested by [START_REF] Mouchi | Potential seasonal calibration for palaeoenvironmental reconstruction using skeletal microstructures and strontium measurements from the cold-water coral Lophelia pertusa[END_REF] from a simple transect (i.e. the elemental variation of Li shown here is real and not an analytical or data processing artefact).

Electron microprobe

Another L. pertusa specimen was mapped with a Cameca SX-Five electron microprobe at the Camparis Service, Paris, France. Acquisition parameters were 15 keV, 200 nA, 500 ms dwell time per point with a 4 µm beam diameter and step size on a surface of 580 by 1352 µm. The map presented here corresponds to Mg wt% (Fig. 8) after calibration using internal reference materials for each measurement. The analytical uncertainty for this element using this protocol is approximately 100 µg.g -1 (2 SD).

Although a dwell time of 500 ms was used for this analytical protocol, which is high compared to common trace element mapping protocols (usually set at around 300 ms), fluctuations in Mg concentrations are still difficult to observe (Fig. 8a) and seem to be lower than the 100 µg.g -1 uncertainty. The visibility of these subtle fluctuations in Mg concentration can be improved using AERYN (Fig. 8b).

The corresponding export files were written with transposed matrices compared to the ones generated by Iolite. The "Transpose" button in AERYN was used to get the same orientation as for the LA-ICP-MS maps. After filtering the data to cut the outliers, the map remains very noisy due to low concentration values. We applied a 7x7 mask to strongly reduce the high-frequency fluctuations and applied a weak detrending as the first lines presented lower values. The resulting map is shown in Figure 8b.

The limits of the bands of high and low concentrations remain difficult to establish, but AERYN permits areas of lower concentrations within high concentration bands to be more clearly visualized, in particular around columns number 40, 80, 170 and 225. The map is therefore ready for interpretation. High concentration bands are wider on the left-hand side (corresponding to the initial stages of skeletal growth) than on the right-hand side of the map. Whereas a detailed interpretation of the geochemical variations presented in these maps is beyond the scope of this paper, if Mg incorporation is, as Sr incorporation, governed by kinetics in the skeleton of L. pertusa, this observation would indicate in the present dataset that early stages of skeletal growth may feature enhanced growth fluctuations. The large bands of high Mg concentration also tend to indicate that most of these early stages present high growth rates and that the low growth rate periods are temporally extremely restricted. These observations are in accordance of those of [START_REF] Mouchi | Potential seasonal calibration for palaeoenvironmental reconstruction using skeletal microstructures and strontium measurements from the cold-water coral Lophelia pertusa[END_REF] and [START_REF] Rollion-Bard | Possible controls on Li, Na, and Mg incorporation into aragonite coral skeletons[END_REF].

Conclusion

AERYN uses data generated from micro-analytical techniques and processed by data reduction software (such as Iolite) to improve the interpretation of elemental distribution maps. This paper has presented an alternative way to the trial-and-error routine of modification of palette thresholds by allowing the user to display the geochemical dataset in 3D and observe data extremes systematically. Moreover, AERYN allows the user to counteract any apparent changes due to instrumental drift, which might not be possible to correct with traditional data reduction schemes. AERYN makes it possible to observe discrete areas of high or low concentration which may be hidden by natural strong fluctuations or analytical artefacts. This program can highlight subtle variations in elemental maps in just a couple of minutes, thereby enhancing the interpretation of elemental maps obtained with a variety of analytical techniques such as LA-ICP-MS, EPMA and SEM-EDS. the use of his Matlab function inpaint_nans. We acknowledge detailed feedback from one reviewer and the editor which helped us improve the original version of the manuscript.

Figure captions

Figure 1 :

 1 Figure 1: Main window of AERYN showing a Sr distribution map measured by LA-ICP-MS on a

Figure 2 :

 2 Figure 2: Coral skeleton Sr content measured by LA-ICP-MS as line scans at the outlier removal step (colour bar represents concentration in µg.g -1 ) using Matlab's "surf" function. a: Map after filtering data out below 0 and above 100,000 µg.g -1 . b: Side view rotated to expose the surface, showing the

Figure 3 :

 3 Figure 3: Detrending applied to the Sr distribution map. a: Selected transect from the map (in black)

Figure 4 :

 4 Figure 4: Sr distribution map after automatic detrending and smoothing with a 5 by 5 mask.

Figure 5 :

 5 Figure 5: Element distribution maps resulting from data reduction using Iolite v. 3.0 (left) and after

Figure caption

  Figure caption

Figure 6 :

 6 Figure 6: First steps for the Li map. a: Initial data after import of the Iolite export file. The resin is

Figure 7 :

 7 Figure 7: Smoothed map of Li concentration, using a 3x3 mask (each value replaced by the mean of

Figure 8 :

 8 Figure 8: Mg (wt.%) map of a L. pertusa specimen acquired by electron microprobe, as obtained by

  Figure 1 Click here to download high resolution image

Acknowledgements

The work presented in this paper was made possible by the ENS Ph.D. programme. The Earth and Natural Sciences Doctoral Studies Programme is funded under the Programme for Research in Third-Level Institutions Cycle-5 and co-funded under the European Regional Development Fund. The authors would like to thank Cora McKenna for her help operating the LA-ICP-MS, and Eric P.

Verrecchia for his advice on some algorithms for the code of AERYN. John D'Errico kindly offered


K Blaum

D Beck

G Bollen

P Delahaye

C Guénaut

F Herfurth

A Kellerbauer

H.-J Kluge

U Köster

D Lunney

S Schwarz

L Schweikhard

C Yazidjian

Study of ground-and low-lying nuclear states by combining

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Study of ground-and low-lying nuclear states by combining laser ionization and Penning trap mass spectrometry 1 Introduction

For the first time we have demonstrated that element-selective laser ionization in combination with high resolution mass spectrometry can be used to prepare isomerically pure beams of 68 Cu and 70 Cu. To this end, the resonance ionization laser ion source (RILIS) [1] and the triple-trap mass spectrometer ISOLTRAP [2,3] at ISOLDE/CERN [4] were used. Together with β-γ coincidence studies this method allowed us to determine the low-energy structure of 70 Cu and to unambiguously identify three β-decaying isomers. By selective laser ionization and high-resolution Penning trap mass measurements it was possible to identify each state and to determine its mass with a relative uncertainty of δm/m ≈ 5 • 10 -8 . Thus, a clear identification was possible, which allowed us to resolve the assignment puzzle in 70 Cu [5].

Production and Ionization of Radionuclides at ISOLDE

The short-lived radionuclides are produced at the on-line isotope separator ISOLDE [4] by bombarding a thick target of a few g/cm 2 with 1.4-GeV proton pulses containing in average about 3 • 10 13 protons. The exotic nuclides diffuse out of the target, which is heated up to 2000 • C and are ionized either by surface, plasma or resonant laser ionization. The 60-keV ion beam is then mass-separated in a magnetic separator and delivered to different experiments.

For the work presented here the laser ion source was applied. The atomic excitation path for the copper isotopes is shown in Fig. 1. It is a resonant excitation scheme followed by ionization through a strong autoionizing (AI) resonance. Two copper vapor lasers with a repetition rate of 11 kHz are used to pump the frequency-doubled dye lasers. The bandwidth of the laser for the first transition at λ = 327.4 nm was 1.2 GHz, which made it possible to resolve the hyperfine structures of the 68,70 Cu isomers [6,7]. For the ionization step a broadband (≈ 24 GHz) dye laser was used.

High-Precision Mass Measurements with ISOLTRAP

The mass measurements on the copper isotopes were performed with the Penning trap mass spectrometer ISOLTRAP [2,3]. It measures the mass m of ions with charge q stored in a homogeneous and stable magnetic field B via the determination of the cyclotron frequency ν c = qB/(2πm). As shown in Fig. 2 ISOLTRAP consists of three traps, a linear radiofrequency Paul trap [8] and two Penning traps [2] located in superconducting magnets with field strengths of 4.7 T and 5.9 T, respectively. The first two traps are operated with helium buffer gas and serve for deceleration, cooling, bunching and isobar purification of the continuous 60-keV ion beam delivered by ISOLDE. The third trap is a precision Penning trap for the mass measurements. The ions' cyclotron frequency ν c is probed by excitation of the ions' motion with a radiofrequency signal [9] and measurement of the time of flight to the micro-channel-plate detector MCP5. The cyclotron resonance (see inset of Fig. 2) is determined by repetition of this sequence and measurement of the time of flight as a function of the frequency of the applied field. The magnetic field strength B is determined via the cyclotron frequency of the well-known 85 Rb + mass for the presented copper measurements. In total, the masses of close to 300 radionuclides have been measured by use of ISOLTRAP with a relative uncertainty of typically δm/m = 10 -8 -10 -7 [10]. The inset shows a cyclotron resonance for 85 Rb + which was used here for the magnetic-field calibration. The line is a fit of the theoretical curve to the data points [11].

Preparation of Isomerically Pure Ensembles

Due to the broad bandwidth and thus the overlapping hyperfine structure the selectivity of the RILIS was not sufficient to separate the 70 Cu isomers effectively. A cyclotron excitation time of T RF = 0.9 s was used in the precision Penning trap. This results in a linewidth of ∆ν c (FWHM) ≈ 0.9/T RF ≈ 1 Hz and a resolving power R = ν c /∆ν c (FWHM) of more than 1 • 10 6 , sufficient to clearly resolve the 68,70 Cu isomers. In order to avoid systematic errors it is not only important to resolve ions with different masses but also to selectively remove the unwanted species, as for instance remaining isobaric contaminations or isomers. This can be achieved with appropriate excitation of their motion. In the experiment presented here, radiofrequency cleaning excitation periods of 3 s were used to prepare isomerically pure ensembles.

Results and Discussion

By scanning the laser frequency of the first transition and by use of a βγcoincidence setup the radionuclides were detected and identified [6,12]. Investigation of the intensities of the individual γ-rays in the β-decay of 70 Cu revealed distinct groups of γ-rays belonging to three different hyperfine-structure patterns as shown in the top part of Fig. 3. This is a clear evidence for the existence of three β-decaying isomers in 70 Cu. Tentative spin values were deduced from the magnetic moments with spin (6 -) for the lowest, (3 -) for the intermediate and 1 + for the highest lying isomer [7].

For the separation of the three isomers and the subsequent mass measurements the laser frequency was tuned to the positions (a, b, c, d) as indicated by arrows in the upper part of Fig. 3. The obtained cyclotron resonances are shown in the lower part of Fig. 3. For position a and d the selectivity of the RILIS did not suffice to prepare isomerically pure samples. Since the resolving power of ISOLTRAP is sufficient to clearly separate the different cyclotron resonances (see a in the lower part of Fig. 3), we applied the above mentioned mass selective cleaning procedure to one or two of the isomers to obtain an isomerically clean cyclotron resonance of the isomer of interest (see resonance d). For the positions b and c the selectivity of the laser ionization was high enough to obtain almost pure ensembles of the (6 -) and (3 -) isomer.

Since the cyclotron frequency is inversely proportional to the mass of the ion, our results unambiguously determined the order of the levels and assigned spins to states, i.e. the (6 -) is the ground state, the (3 -) the first and the 1 + the second isomeric state, as proposed in [7]. The mass value differences are in excellent agreement with the excitation energies obtained in the decay studies The assignment of the different spin states to ground and isomeric states results from the ISOLTRAP measurement (see text). The tentative spin values were taken from [7]. The solid lines are fits of the theoretically expected line shape [11] for two (a) or one (b, c and d) resonance curves to the data points.

and are given in [5]. The same procedure was applied to 68 Cu to prepare pure ensembles of the ground and isomeric state and to perform background-free mass measurements. These data exemplify the strength of the unique combination of resonance ionization and Penning trap mass spectrometry to produce isomerically pure ensembles of short-lived radionuclides. A further exploration of these techniques, including post-acceleration with e.g. REX-ISOLDE [13], will open a whole field of radioactive beam research with isomer beams.

Figure 1 :

 1 Figure 1: Excitation scheme used for double-resonance ionization of copper.

Figure 2 :

 2 Figure 2: Schematic drawing of the mass spectrometer ISOLTRAP including the RFQ trap and the preparation and precision Penning traps. Micro-channel-plate (MCP) detectors are used to monitor the ion transfer (MCP1 and MCP3) as well as to record the time-of-flight resonance (MCP5) for the determination of the cyclotron frequency.The inset shows a cyclotron resonance for 85 Rb + which was used here for the magnetic-field calibration. The line is a fit of the theoretical curve to the data points[11].

  d

Figure 3 :

 3 Figure 3: Top: Measured (βγ) radioactivity of the ground (squares), first (bullets) and second (triangles) isomeric state of 70 Cu + versus laser frequency of the first transition. Bottom: Time-of-flight resonance curves as a function of cyclotron frequency for the laser settings marked with arrows (a, b, c, d) in the top figure. (a) shows the isomers resolved but not separated while (bd) shows the sequential separation of the isomeric states.The assignment of the different spin states to ground and isomeric states results from the ISOLTRAP measurement (see text). The tentative spin values were taken from[7]. The solid lines are fits of the theoretically expected line shape[11] for two (a) or one (b, c and d) resonance curves to the data points.

Acknowledgments

The authors would like to thank the ISOLDE staff for their continuous support. This work was supported by the BMBF (06GF151) and by the European Commission (IHRP HPRI-CT-1999-00018, TARGISOL HPRI-CT-2001-50033, and NIPNET HPRI-CT-2001-50034).


Nicolas Holzschuch

Romain Pacanowski

A Physically-Based Reflectance Model Combining Reflection and Diffraction

Keywords: Appearance, material, Reflectance, BRDF, Measured materials, Diffraction

Reflectance properties express how objects in a virtual scene interact with light; they control the appearance of the object: whether it looks shiny or not, whether it has a metallic or plastic appearance. Having a good reflectance model is essential for the production of photorealistic pictures. Measured reflectance functions provide high realism at the expense of memory cost. Parametric models are compact, but finding the right parameters to approximate measured reflectance can be difficult. Most parametric models use a model of the surface micro-geometry to predict the reflectance at the macroscopic level. In this paper, we show that this micro-geometry causes two different physical phenomena: reflection and diffraction. Their relative importance is connected to the surface roughness. Taking both phenomena into account, we develop a new reflectance model that is compact, based on physical properties and provides a good approximation of measured reflectance.

Un modèle de réflectance combinant diffraction et réflexion Résumé :

Les propriétés de réflexion expriment la façon dont les objets interagissent avec la lumière dans les scènes virtuelles. Elles contrôlent l'apparence de l'objet: s'il apparaît brillant ou non; son aspect métallique ou plastique ? Le modèle de réflexion (BRDF) est primordial afin d'obtenir des images photo-réalistes. Les modèles de réflexion mesurés fournissent un haut degré de réalisme au détriment du coût mémoire. Les modèles paramétrés sont compacts mais il reste difficile de trouver les bonnes valeurs des paramètres à partir des réflectances mesurées. Dans ce rapport de recherche, nous montrons que deux phénomènes physiques sont présents dans les mesures de réflexion : la réflexion à proprement parler auquel s'ajoute le phénomène de diffraction. En prenant les deux en compte nous présentons un modèle de BRDF fournit une très bonne approximation des mesures tout en restant compact en mémoire. Les artistes peuvent modifier les paramètres du modèle qui sont reliés aux propriétés de la surface afin de créer de nouveau matériaux.

Mots-clés : Modèles de matériaux, simulation de l'éclairage, réflectance, diffraction

Introduction

Reflectance properties express how materials interact with light. They are responsible for the overall aspect of objects in virtual scenes: whether they look shiny, metallic, plastic. . . Measured reflectance functions [START_REF] Matusik | A data-driven reflectance model[END_REF] reproduce the appearance of the material well, but have storage issues (33 MB for a single isotropic material). They also make it more difficult to edit materials or to use importance sampling.

Parametric reflectance models, such as Phong, Lafortune or Cook-Torrance, are more compact and easier to edit. They are also easier to combine with importance sampling in Monte-Carlo integration. On the other hand, it is difficult for users to find the parameters that correspond to the visual characteristic of a specific material. To combine the advantages of both approaches, we want to compute the parameters for one BRDF model so it produces the same appearance as a given measured material. This would provide a compact representation and allow importance sampling. It would also make it easier to edit material properties, to make them smoother or rougher.

Many parametric models are based on a physical representation of the surface micro-geometry and how it interacts with incoming light. The Cook-Torrance model [START_REF] Cook | A reflectance model for computer graphics[END_REF] is very popular in Computer Graphics. Based on the assumption that surface micro-geometry acts as a set of specular mirrors, it connects the normal distribution of the surface at the microscopic level with its visual appearance at the macroscopic level. Another model, widely used in the Optics community, computes diffraction effects caused by differences in height in the surface micro-geometry, and predicts visual appearance from the frequency content of the height distribution. In this model, the width of material response intrinsically depends on the light wavelength.

It is interesting to compare measured reflectance functions with what the models predict. Butler et al. observed that measured material response is well approximated with a micro-facet model [START_REF] Butler | Experimental analysis of bidirectional reflectance distribution function cross section conversion term in direction cosine space[END_REF] but that the width of material response appears to vary with the wavelength [START_REF] Butler | Experimental measurement and analysis of wavelengthdependent properties of the BRDF[END_REF], in contradiction with the model. Löw et al. [START_REF] Löw | Brdf models for accurate and efficient rendering of glossy surfaces[END_REF] observed that the behaviour of many measured reflectance corresponds to the predictions of the diffraction model for low intensities. Several researchers, e.g. [NDM05, [START_REF] Burley | Physically-based shading at Disney[END_REF][START_REF] Löw | Brdf models for accurate and efficient rendering of glossy surfaces[END_REF], have observed that some measured materials appear to be the sum of two different behaviors.

In this paper, we provide a physical explanation to all these observations, even where they appear to contradict the physical models: when modelling the influence of surface micro-geometry on material reflectance, diffraction and reflection are not mutually exclusive. The Harvey-Shack diffraction model [START_REF] Harvey | Light-Scattering Characteristics of Optical Surface[END_REF] explicitly predicts the energy repartition between specular reflection and diffraction effects. Integrating this energy repartition, we propose a new reflectance model that combines microfacet reflection and height-based diffraction.

With our model, reflection effects, governed by the micro-geometry normals, dominate near the specular direction. Diffraction effects, governed by the micro-geometry height distribution, are dominant for wide-angle scattering, away from the specular direction. The relative importance of these two phenomena is governed by the variance of the height distribution. All these effects come directly from the full application of Generalized Harvey-Shack theory. For smooth surfaces, the two effects are well separated and appear as two different lobes. For rougher surfaces, the two effects overlap and result in a single lobe, whose width changes with the wavelength.

For conductors, our model express reflectance properties without a diffuse term. We introduce two variations of our model for materials that explicitly cause a diffuse effects: sub-surface scattering and multi-layered plastic materials.

All parameters in our model are connected to physical quantities. The color of the material comes solely from its index of refraction. The shape of the reflectance lobe comes from the normal distribution function, the variance of the surface height and the frequency spectrum of the height distribution. We validate our model by using it to approximate measured material reflectance functions. It is a good predictor of material behavior for smooth and moderately rough surfaces.

In the next section, we review some of the previous work on material reflectance, both in the Computer Graphics and Optical Engineering communities. We then describe in Section 3 the scientific background: the micro-facet reflectance models and the Generalized Harvey-Shack diffraction theory. In Section 4, we describe our model, combining reflection and diffraction, and fully describe each component. In Section 5, we present our method to extract model parameters from measured reflectance functions, and discuss the results. In Section 6, we conclude and present avenues for future work.

Previous Work

Micro-Facet Model

The micro-facet model [START_REF] Torrance | Theory for off-specular reflection from roughened surfaces[END_REF][START_REF] Cook | A reflectance model for computer graphics[END_REF] assumes that surface micro-geometry is made of specular micro-facets. It predicts the overall surface appearance from the probability distribution of these micro-facet normals. Existing models include Gaussian [START_REF] Cook | A reflectance model for computer graphics[END_REF], rational fraction [TR75, WMLT07], fraction to the power p [LKYU12, Bur12], Shifted-Gamma Distribution [START_REF] Mahdi | Accurate fitting of measured reflectances using a Shifted Gamma micro-facet distribution[END_REF], exponential of a power function [START_REF] Brady | genBRDF: Discovering new analytic BRDFs with genetic programming[END_REF].

The shadowing and masking term is essential for energy conservation in the micro-facet model. Assuming that microfacets position and orientation are independent, Smith [START_REF] Smith | Geometrical shadowing of a random rough surface[END_REF] computes the shadowing/masking term from two successive integrations of the slope distribution. Heitz [START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF] shows that this is the most physically consistent method to compute the shadowing and masking term, and provides an improved shadowing term taking into account correlation between input and output directions.

Diffraction Models

Diffraction models focus on diffraction effects caused by variations in surface height. They are widely used in the Optical Engineering Community, for example to measure imperfections on smooth surfaces, such as mirrors. There are several models expressing the relationship between surface height distribution and reflectance properties: the Rayleigh-Rice vector perturbation theory is designed for smooth surfaces and wide-angle scattering; the Beckmann-Kirchhoff [START_REF] Beckmann | The scattering of electromagnetic waves from rough surfaces[END_REF] and Modified Beckmann-Kirchhoff [START_REF] Vernold | A modified beckmann-kirchhoff scattering theory[END_REF] theories express the BRDF as an infinite sum; they are valid for a large range of surfaces, but for small angles. For very smooth or very rough surfaces, the infinite sum is approximated by a single term. The Harvey-Shack model [START_REF] Harvey | Light-Scattering Characteristics of Optical Surface[END_REF] express the diffraction effects using a surface transfer function. For smooth surfaces, this transfer function is the Power Spectral Distribution of surface height. Krywonos [START_REF] Krywonos | Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction[END_REF] extended this with the Generalized Harvey-Shack model, where the surface transfer function is recomputed for each incoming direction. This model is valid for all surfaces, rough and smooth, and for all angles. Krywonos [START_REF] Krywonos | Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction[END_REF] compared models with measured reflectance and showed that the Generalized Harvey-Shack model provides a better approximation than other diffraction models.

He et al. [START_REF] He | A comprehensive physical model for light reflection[END_REF] derived the most comprehensive reflectance model for isotropic surfaces, based on Beckmann-Krichhoff theory; it takes into account all scales in the micro-geometry in a single model. Stam [START_REF] Stam | Diffraction shaders[END_REF] introduced an anisotropic reflectance model, also based on the Beckmann-Krichhoff model. His model reverts to He's for isotropic surfaces. These models remain impractical for rendering and fitting since their number of terms grows quickly with the surface roughness.

In this paper, we build on the Generalized Harvey-Shack model, as it models diffraction effects using a single term instead of an infinite sum and has been shown to agree better with actual reflectance measurements.

Comparison with Measured Materials

Butler et al. studied the connection between these two models. They found that measured reflectance functions are well approximated by the micro-facet model in the incident plane [START_REF] Butler | Experimental analysis of bidirectional reflectance distribution function cross section conversion term in direction cosine space[END_REF], but that the normal distribution appears to change with the wavelength [START_REF] Butler | Experimental measurement and analysis of wavelengthdependent properties of the BRDF[END_REF]. They studied the connections between the two models and proposed an approximation for the ratio of the Fresnel term for reflection and the color term for diffraction [START_REF] Butler | Comparison of microfacet BRDF model elements to diffraction BRDF model elements[END_REF].

The two models have different parameters and thus predict different shapes for reflectance lobes outside of the incident plane. Löw et al. [START_REF] Löw | Brdf models for accurate and efficient rendering of glossy surfaces[END_REF] showed that the shape of the reflectance lobe in many measured materials is consistent with the predictions of diffraction theory. They also showed that the diffraction model provides a good fit with measured reflectance. To conduct their fitting, they removed the wavelength dependency, in contradiction with the diffraction model. They also replaced the color term for diffraction Q with the Fresnel term for reflection F.

Holzschuch and Pacanowski [START_REF] Holzschuch | Identifying diffraction effects in measured reflectances[END_REF], through an empirical study, showed that measured reflectance appear to be a combination of the two phenomena: diffraction for wide-angle scattering, micro-facet reflection close to the specular direction. In this paper, we provide a physical explanation for their empirical observation. In later work [START_REF] Holzschuch | A physically accurate reflectance model combining reflection and diffraction[END_REF], they showed that summing the two components provides a good approximation to all measured reflectance. Compared to their work, we remove the need for a diffuse lobe, explicitly model the energy repartition between the two lobes and predict the color of each lobe using only the material index of refraction, reducing the number of parameters strongly.

Dong et al [START_REF] Dong | Predicting appearance from measured microgeometry of metal surfaces[END_REF] compared measured material reflectance and measured micro-geometry; they found that measured normals are a good predictor of the specular peak of the reflectance lobe, if the surface features are filtered to remove subwavelength details. They also compared measured reflectance with Kirchhoff diffraction theory, using a single term.

Matusik et al. [START_REF] Matusik | A data-driven reflectance model[END_REF] measured and released reflectance properties for a large range of materials. We use their database in our tests. Ngan et al. [START_REF] Durand | Experimental analysis of BRDF models[END_REF] tried fitting parametric BRDF models to this measured data. They found the best fits for the Cook-Torrance [START_REF] Cook | A reflectance model for computer graphics[END_REF] and Lafortune [START_REF] Lafortune | Non-linear approximation of reflectance functions[END_REF] models. They also reported that many specular materials are poorly approximated with a single lobe. The quality of the fit improves with multiple lobes, but the fitting process becomes unstable.

Ashikhmin and Premože [AP07] approximated measured BRDFs using back-scattering: if input and output directions are equal, the entire BRDF can be expressed as a function of the half-vector. By storing this function, they get a compact BRDF model, that fits measured data very well. We use their intuition in our fitting process to provide a first approximation. Brady et al. [START_REF] Brady | genBRDF: Discovering new analytic BRDFs with genetic programming[END_REF] applied genetic programming to discover new analytical BRDF models from measured BRDFs. We use their distribution as the basis of our reflectance lobe. Our contributions include the normalization constant, the Smith shadowing term and the importance sampling method for this distribution.

3 Scientific Background

Parameters

We consider a given material; the surface of this material is described by a height field h, as well as a normal distribution function D(m). We are interested in the reflectance properties and whether they can be predicted from the surface properties, h and D.

The BRDF ρ(i, o) encodes the way the material reflects light. It is defined as the ratio of outgoing radiance in direction o by incoming irradiance along i.

The two competing models for reflectance properties, micro-facet model and diffraction models, use different parameters, making direct comparison difficult. Micro-facet models depend mostly on the half-vector h, defined as:

h = (i + o)/ i + o , (1)
while diffraction models depend on a vector f , defined as the projection on the tangent plane of the difference between the reflected incoming direction and the outgoing direction, divided by the wavelength (see Figure 2):

f = 1 λ proj(refl(i) -o) (2)
For isotropic materials, diffraction effects only depend on f , the norm of f , which can be expressed using the half-angle parametrization [START_REF] Rusinkiewicz | A new change of variables for efficient brdf representation[END_REF]:

f = 2 λ sin θ h cos θ d (3) n i λf R(i) o
Figure 2: Diffraction effects depend on f , projected vector between reflected incoming direction and outgoing direction, divided by the wavelength.

Micro-Facet Model

The micro-facet model [START_REF] Cook | A reflectance model for computer graphics[END_REF] is widely used in Computer Graphics. It is based on the hypothesis that surface details are much larger than light wavelength, so we can apply the rules of optical geometry. The model has been described in great details in the literature [CT82, WMLT07, Hei14] so we focus here on the key points.

The model assumes that the surface is made of micro-facets, with a normal distribution function D. Each micro-facet acts as a specular mirror; its reflectance is a Dirac delta function, multiplied by the Fresnel term for reflection:

ρ µF (i, o) = F(i, o)δ(refl(i), o) (4)
The Fresnel term F is responsible for the color of the facet. It depends on the index of refraction η of the material; the index of refraction depends on the wavelength, and this dependency is responsible for the color of the material. For dielectrics (transparent materials), η is a real number, related to the speed of light inside the material. For conductors, such as copper and nickel, η is a complex number, with real and imaginary parts:

η = n + ik.
Summing the contributions of all micro-facets results in the Cook-Torrance reflectance model:

ρ CT (i, o) = F(θ d )D(θ h )G(i, o) 4 cos θ i cos θ o (5)
The shadowing-masking term G expresses the probability for light to be blocked before or after the specular reflection on the micro-facet. It is usually computed using Smith's method, through a double integration of D [Smi67, [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF][START_REF] Heitz | Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs[END_REF]. First, we approximate G as the product of separate, one-dimensional shadowing and masking functions, then compute these using Λ function:

G(i, o) ≈ G 1 (θ i )G 1 (θ o ) (6) G 1 (θ) = 1 1 + Λ(θ) (7)
Λ is computed through a double integration of D, through auxiliary functions P 22 and P 2 :

P 22 tan 2 θ = D(θ)χ [0, π 2 ] (θ) cos 4 θ (8) P 2 (r) = ∞ -∞ P 22 (r 2 + q 2 ) dq (9) Λ(θ) = ∞ 1/ tan θ (r tan θ -1)P 2 (r) dr (10)
Note that D and G only depend on the geometry of the micro-facets, not on the wavelength of the incoming light. The only wavelength dependency comes from the Fresnel term F. For this model, the reflection lobes for the different wavelengths should be scaled version of each other, with the same relative width.

θ i + cos θ o )h(x, y) [Kry06]
. Heights h are relative to the average surface level.

Color and Fresnel term for reflection

The material color comes from the Fresnel term. For unpolarized light, the Fresnel term is computed using these equations:

F s (θ) = cos θ -η 2 -sin 2 θ 1/2 cos θ + η 2 -sin 2 θ 1/2 2 (11) F p (θ) = η 2 cos θ -η 2 -sin 2 θ 1/2 η 2 cos θ + η 2 -sin 2 θ 1/2 2 (12) F = F s + F p 2 (13)
These equations remain valid even if η is a complex number, as happens for conductors. F s , F p and F are always real numbers.

In the presence of polarized light, F s gives the material reflectance for perpendicular polarized light, and F p the reflectance for parallel-polarized light.

Diffraction: Generalized Harvey-Shack Theory

In this section, we review the basis of the diffraction theory for predicting material reflectance. We focus on the Harvey-Shack theory, because it is applicable to a wide range of configurations and has been found to provide a good approximation to measured reflectance. In the diffraction theory, we consider the distance travelled by light as it is reflected on the surface. Because the surface is not flat, light being reflected at a point of height h(x, y) has travelled a shorter distance than if it was reflected at the reference height (see Figure 3). This difference in optical path length (OPD) is:

OPD = (cos θ i + cos θ o )h(x, y) (14)
This results in phase variations at the reflected wavefront of (2π/λ)(cos θ i + cos θ o )h(x, y). Averaging these phase variations gives the reflectance function.

According to the Harvey-Shack theory, reflectance is separated in two different lobes: a peak in the direction of specular reflection, corresponding to a Dirac delta function, surrounded by a halo of scattered light [START_REF] Harvey | Light-Scattering Characteristics of Optical Surface[END_REF][START_REF] Krywonos | Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction[END_REF]. The relative intensity of these two components depend on the surface roughness σ 2 s , defined as the variance of height distribution h:

ρ diff (i, o) = Aδ(refl(i), o) + (1 -A)S GHS ( f ) (15)
with

A = e -(2π σs λ (cos θ i +cos θ o )) 2 (16)
where S GHS is an angle spread function. By construction, this equation guarantees conservation of energy: if the surface is rough, σ s /λ is large; the energy of the diffraction lobe increases, and the energy of the specular peak decreases. The relative importance of the specular peak and diffraction lobe depends on the directions: at grazing angles, cos θ i ≈ 0 and A is almost equal to 1; most of the energy goes into the specular reflection. The influence of diffraction is stronger at normal incidence, where cos θ i ≈ 1.

Krywonos [START_REF] Krywonos | Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction[END_REF] showed that the Generalized Harvey-Shack theory provides a good fit with measured reflectance, outside of the specular peak. It provides a better match than competing diffraction models, such as Rayleigh-Rice and the original Harvey-Shack theory. It provides equivalent results to the Modified Beckmann-Kirchoff theory.

Smooth Surface Linearization

If the reflecting surface is smooth, we can assume that σ s /λ 1 and linearize the terms in Equation 16:

A ≈ 1 -2π σ s λ (cos θ i + cos θ o ) 2 (17)
The diffraction lobe becomes:

ρ diff. (i, o) = 4π 2 σ 2 s λ 2 (cos θ i + cos θ o ) 2 S GHS ( f ) (18)
Under the smooth surface linearization, we compute the diffraction function from the Fourier transform of the autocovariance function of the surface height, which is the Power Spectral Density (PSD) of the surface:

ρ diff. (i, o) = 4π 2 λ 4 (cos θ i + cos θ o ) 2 PSD( f ) (19)
We recognize here the usual formula for Harvey-Shack diffraction for smooth surfaces. This expression depends on the wavelength in two places: the intensity of the diffraction lobe decreases with 1/λ 4 ; the width of the lobe increases with λ, since f contains a 1/λ term. The diffraction lobe for blue colors is sharper, but with higher intensity, while the diffraction lobe for red colors spreads more, but with lower maximum intensity. The Generalized Harvey-Shack diffraction theory is not limited to this smooth surface approximation. In the generic case, the angle spread function S GHS changes with incoming direction θ i and is computed by a Fourier transform of surface height.

K-correlation Model

If the surface height distribution follows a Gaussian, then the PSD should also be Gaussian, making computations simpler. Most optical surfaces however do not have a Gaussian distribution of heights. They tend to have an inverse power law falloff at large spatial frequencies. For these surfaces, we use the K-correlation model, or abc model, to express the power spectral density:

PS D( f ) = 1 2 √ π Γ ((c + 1)/2) Γ (c/2) ab 1 + b 2 f 2 c+1 2 (20)
Integrating this distribution gives the surface roughness σ 2 s . It is finite only for c > 1 and is equal to:

σ 2 s = √ π c -1 Γ ((c + 1)/2) Γ (c/2) a b (21)
We use this surface roughness to express the Power Spectral Density in a simpler expression:

PS D( f ) = c -1 2π σ 2 s b 2 1 + b 2 f 2 c+1 2 (22)

Renormalization

The angle spread function S GHS used in Equation 16 is defined over the entire plane for its parameter, f . But f is limited in its variation, as the projections of incoming and outgoing directions must fall over the unit disc in the plane. For energy conservation, S GHS must be renormalized by dividing by the integral over the set of possible values for vector f :

Ŝ GHS = σ 2 0 σ 2 rel S GHS (23)
σ 2 0 = ∞ -∞ ∞ -∞ S GHS ( f ) d f x d f y (24) σ 2 rel (θ i ) = o ≤1 S GHS ( f ) d f x d f y (25)
This renormalization ensures that the energy of the diffraction lobe is constant. Usually, angle spread functions are normalized so that σ 2 0 = 1. In that case, Equation 23 is equivalent to:

Ŝ GHS = 1 σ 2 rel S GHS (26)

Color and Polarization

As for specular reflection and micro-facet models, the color of the diffraction lobe comes from an index of refraction dependent term. This term is noted Q; it is different from the Fresnel term for reflection, F, defined in Section 3.2.1 [START_REF] Krywonos | Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction[END_REF][START_REF] Butler | Comparison of microfacet BRDF model elements to diffraction BRDF model elements[END_REF].

For unpolarized light, it is defined as:

φ = φ o -φ i -π Q ss = (η 2 -1) cos φ (cos θ i + η 2 -sin 2 θ i )(cos θ o + η 2 -sin 2 θ o ) Q sp = (η 2 -1) η 2 -sin 2 θ o sin φ (cos θ i + η 2 -sin 2 θ i )(η 2 cos θ o + η 2 -sin 2 θ o ) Q ps = (η 2 -1) η 2 -sin 2 θ i sin φ (η 2 cos θ i + η 2 -sin 2 θ i )(cos θ o + η 2 -sin 2 θ o ) Q pp = (η 2 -1)× η 2 -sin 2 θ i η 2 -sin 2 θ o cos φ -η 2 sin θ i sin θ o (η 2 cos θ i + η 2 -sin 2 θ i )(η 2 cos θ o + η 2 -sin 2 θ o ) Q = |Q ss | 2 + |Q sp | 2 + |Q ps | 2 + |Q pp | 2 (27)
Q depends on the azimuthal angle between reflected incoming direction and outgoing direction. With polarized light, part of the parallel-polarized light is reflected at perpendicular-polarized, and reciprocally. The Q sp and Q ps terms express this. Along the specular reflected direction, where θ i = θ o and φ o = φ i + π, we have Q sp = 0 and Q ps = 0, and Q = 2F. Outside of this specific case, the two color terms are substantially different.

Our Model

Our reflectance model is based on the following observation: reflectance and diffraction effects are not mutually exclusive. The Generalized Harvey-Shack theory explicitly allocates energy to the reflection and diffraction lobes, combining a specular reflection with a diffraction lobe (see Equation 16).

Most surfaces in Computer Graphics are not smooth, at least not in the sense used by the Optics community. They have variations in surface normal. These variations cause a larger distribution of specular reflections, instead of just a Dirac. The Cook-Torrance model is designed precisely to connect the variations of surface normals with the general reflectance properties of the surface.

Derivation from micro-surface properties

More precisely, we assume that the diffraction model is valid at the micro-geometry level. Each micro-facet reflects light according to Equation 15. We assume that we have a micro-facet normal distribution D(m), describing the statistical distribution of micro-facet normals m. This micro-facet normal distribution has a larger scale than the micro-geometry responsible for diffraction.

We compute the surface reflectance at the macroscopic level by integrating the contributions over all visible microfacets:

ρ(i, o, n) = i • m i • n ρ diff (i, o, m) o • m o • n G(i, o, m)D(m) dω m (28)
The integral gives a reflectance model with two contributions:

ρ(i, o, n) = ρ refl. (i, o) + ρ diff. (i, o) (29)
The reflection part is simpler; the dirac inside the integral transforms the integral into a function evaluation, after taking into account the proper Jacobian for a change of variable (see, e.g. [START_REF] Walter | Microfacet models for refraction through rough surfaces[END_REF]). Angles between micro-facet normal and incoming or outgoing direction become angles between half-vector and incoming or outgoing directions, that is θ d :

ρ refl. (i, o) = e -(4πσ s cos θ d ) 2 F(θ d )D(θ h )G(i, o) 4 cos θ i cos θ o (30)
The diffraction part is expressed as:

ρ diff. (i, o) = i • m i • n (1 -A)S GHS ( f ) o • m o • n G(i, o, m)D(m) dω m (31)
We cannot compute this integral; we approximate it using a different angle spread function, S GHS ( f ):

ρ diff. (i, o) ≈ (1 -A)S GHS ( f )G(i, o) (32)
(33)

This gives the complete model:

ρ(i, o) = ρ refl. (i, o) + ρ diff. (i, o) (34)
ρ refl. (i, o) = e -(4πσ s cos θ d ) 2 F(θ d )D(θ h )G(i, o) 4 cos θ i cos θ o (35) ρ diff. (i, o) = (1 -A)S GHS ( f )G(i, o) (36)
where

A = e -(2π σs λ (cos θ i +cos θ o )) 2 (37)
Compared to the Generalized-Harvey-Shack diffraction model, we make the following assumptions: we do not use the smooth surface approximation, A is expressed as an exponential instead of the linearized version; we represent the combined effect of surface normal and diffraction using an angle spread function S GHS and assume it remains the same for all θ i . However, it is not necessarily the Power Spectral Distribution of the surface height.

The color for both micro-facet and diffraction terms comes from the index of refraction of the material, through the Fresnel term F for the micro-facet component, and through the reflectance Q for diffraction. Figure 1 shows examples of the contributions of diffraction and reflection for several materials. You can see slight differences in color between the components.

In the following sections, we review the components of the model: the micro-facet distribution for the Cook-Torrance model and the angle spread function for the diffraction component. We also combine our surface reflectance model with two multi-layer models: sub-surface scattering and plastics with specular varnish.

Micro-Facet distribution

For our micro-facet distribution, we use a Generalized Beckmann distribution; it can model the sharp peak we observe on measured reflectance for many surfaces (see Figures 11 and12):

D(θ) = χ [0, π 2 ] (θ) cos 4 θ P 22 tan 2 θ (38) P 22 (x) = p πβ 2 Γ(1/p) e -x β 2 p (39) n θ i θ' i θ o θ' o

Multiple scattering events

Participating media

Single scattering

Figure 4: We also consider multiple layer materials. In the subsurface model, light is reflected on the surface as well as refracted inside the underlying participating media, where it is scattered before leaving again through the refractive interface.

It has two parameters: β controls the width of the peak, while p controls its kurtosis, that is whether it drops sharply at the origin or not. This distribution is similar to the "BRDF model A" obtained by Brady et al. [START_REF] Brady | genBRDF: Discovering new analytic BRDFs with genetic programming[END_REF] by analysing measured BRDFs with genetic programming; we added the 1/ cos 4 θ term, which is required for analytical integration. This allows us to provide normalization, Smith shadowing and importance sampling.

We compute the shadowing term G for this distribution using Smith's method [Smi67, WMLT07, Hei14], as described in Section 3.2. For this distribution:

P 2 (r) = p πΓ(1/p) ∞ 0 e -(r 2 +q 2 ) p dq (40) Λ(β tan θ) = ∞ 1/β tan θ
(rβ tan θ -1)P 2 (r) dr (41)

G 1 (θ) = 1 1 + Λ(β tan θ) (42)
We observe that Λ, and therefore G 1 , can be expressed as a function of β tan θ and p, reducing the number of parameters to two. We precompute G 1 and store it in a 2D array, to be accessed at render time (a single 2D array for all values of the three parameters, with a memory cost of 10 kB) Finally, we compute G as:

G(i, o) = G 1 (θ i )G 1 (θ o ).

Diffraction component

We did not use the smooth surface approximation of the Generalized Harvey-Shack model; as a consequence, S GHS is not necessarily the Power Spectral Distribution of the surface height. We nevertheless chose to express it using the K-correlation model described in Section 3.3.2, as it has good qualities. We used the simplified form defined in Equation 22, and renormalize it using Equation 26:

S GHS ( f ) = σ 2 s σ 2 rel c -1 2πλ 2 b 2 1 + b 2 f 2 c+1 2 (43)
This component has three parameters: b, corresponding to the width of the diffraction lobe, c, corresponding to the falloff and σ 2 s , the surface roughness term.

Multi-layer Models

Some materials have multiple layers: for example diffuse plastic with a layer of transparent varnish, or a semi-precious stone with participating media. These layers contribute to the overall response of the material. We combine our surface model with two different multi-layer materials: sub-surface scattering BRDF and diffuse plastics. The first is the BRDF approximation of subsurface scattering, introduced by Jensen et al. [START_REF] Jensen | A practical model for subsurface light transport[END_REF] (see Figure 4). The second is the multi-layered model introduced by Weidlich and Wilkie [START_REF] Weidlich | Arbitrarily layered micro-facet surfaces[END_REF] with a transparent layer of varnish over a diffuse layer (see Figure 5).

In both models, we use our reflectance model only for direct reflection by the top interface. For transmission through this interface, we approximate it as a flat specular surface. The transmission term for light passing through this interface at i and With Importance Sampling Without Importance Sampling leaving in direction o is:

n θ i θ' i θ o θ' o

Transparent varnish layer

Diffuse opaque layer

T = (1 -F(θ i ))(1 -F(θ o )) ( 44
)
where F is the Fresnel term for reflection, defined in Equation 13. In both models, we also use the Fresnel diffuse reflectance, F dr , the average of Fresnel reflectance times the cosine of the angle with the normal over all directions. We use the rational approximation given by d'Eon and Irving [START_REF] D'eon | A quantized-diffusion model for rendering translucent materials[END_REF].

The sub-surface BRDF model combines single and multiple scattering. Assuming a material with albedo α and an isotropic phase function, the resulting lobe is:

single = αT 4π 1 cos θ i + cos θ o (45) A F = 1 -F dr 1 + F dr (46) multi = αT 2π e - √ 3(1-α) 1 + e -4 3 A F √ 3(1-α) (47)
ρ SSS = single + multi (48)
where θ i and θ o are the angles of the refracted rays with the surface normal (see Figure 4).

In the plastic BRDF model, we assume that we have a diffuse layer of reflectance ρ d under a transparent layer. We account for light entering through the interface, being reflected on the diffuse layer, and leaving through the interface. We also account for possible multiple bounces between the diffuse layer and the refractive interface. The resulting component is:

ρ plastic = T ρ d π 1 - ρ d π F dr 1 η -1 (49)

Importance sampling

The A coefficient express the energy repartition between the reflection and diffraction lobes. It depends on the incoming and outgoing directions and on the σ s /λ ratio (see Figure 8). For σ s /λ > 0.3, it is close to 0 on most of the parameters range, and close to 1 for grazing angles. We need to take it into account for importance sampling.

We perform importance sampling in three steps: first, we select whether to importance sample the diffuse lobe (if there is one) or the glossy lobe. For the glossy lobe, we select between the micro-facet and the diffraction lobe depending on the incoming direction. Once we have picked a lobe, we sample it using the appropriate strategy. Figure 7 displays the impact of importance sampling.

Approximated Energy for each Lobe

To compute the approximate energy for each lobe, we begin by computing the values of the Fresnel reflection coefficient at normal incidence:

F 0 = n 2 + k 2 -2n + 1 n 2 + k 2 + 2n + 1 (50) (51)
We use this as the approximate energy of the combined micro-facet and diffraction lobes:

E C = F 0 .
For multi-layers material, we also compute the energy of the diffuse lobe: For the plastic model:

E D Plastic = (1 -F 0 ) 2 ρ d 1 - ρ d π F dr 1 η -1 (52)
For the sub-surface scattering model:

E D SSS = α(1 -F 0 ) 2 1 + R d 2 ( 53
)
where

R d = e - √ 3(1-α) 1 + e -4 3 A F √ 3(1-α) (54)

Lobe Selection

To importance sample the BRDF, we begin by selecting which lobe to sample. We take a random variable u 1 in [0, 1) to select between the diffuse lobe and the combined micro-facet and diffraction lobes:

p C = E C /(E CT + E D ) ( 55
)
p D = E D /(E C + E D ) (56)
We sample the combined micro-facet and diffraction lobes if u 1 < p C and the diffuse lobe otherwise. We use cosine-weighted sampling for the plastic or sub-surface scattering diffuse lobes.

n i d p d p R(i) o i p o p i p -i p r φ i φ
Figure 9: To importance sample the diffraction lobe, we start from the projection of the reflection of incoming direction i, -i p . We compute the difference between projected incoming direction and outgoing directions, d p and use it to compute the outgoing direction o.

We then compute A imp , using the known incoming direction:

A imp = e -(2π σs λ cos θ i ) 2 (57)
We take a second random variable u 2 in [0, 1), and sample the micro-facet lobe if u 2 < A, and the diffraction lobe otherwise.

To avoid using too many random variables, we renormalise each random variable after use.

Importance sampling for the micro-facet lobe

To importance sample the micro-facet lobe, we take two uniform random variables (u 1 , u 2 ) in [0, 1) and use them to build the half-vector h through (θ h , φ h ):

φ h = 2πu 1 ( 58
)
θ h = arctan β γ -1 u (1/p, u 2 ) p 2 (59)
where γ -1 u is the inverse of the normalized upper incomplete gamma function. We then build the outgoing direction o from h.

Diffraction Lobe Importance Sampling

To importance sample the diffraction lobe, we work in the unit disk corresponding to the projections of all directions (see Figure 9). We begin with the projection of the incoming direction, i p = (sin θ i cos φ i , sin θ i sin φ i ). The projection of the reflected incoming direction is -i p . We build a 2D vector d p , the difference between -i p and the projection of the outgoing direction o p . From this, we extract the projection o p = -i p + d p and then the outgoing direction o = o p + 1 -o p 2 n.

Our importance sampling method does not give equal probability for all points in the disc. We compensate for this point by dividing the sampling weight for each point by its probability. We take two uniform random variables (u 1 , u 2 ) in [0, 1) and use them to build d p in polar coordinates (r, φ):

r = λ b (1 -Mu 1 ) -2 c-1 -1 ( 60
)
where

M = 1 -1 + b 2 λ 2 (1 + sin θ i ) 2 -c-1 2 (61)
This gives r ∈ [0, 1 + sin θ i ]. We then select φ so that o p remains inside the unit disc:

φ max = arccos max -1, r 2 + sin 2 θ i -1 2r sin θ i (62) φ = φ i + π + (2u 2 -1)φ max (63)
Our importance sampling of the diffraction lobe is simpler than the method provided in [START_REF] Löw | Brdf models for accurate and efficient rendering of glossy surfaces[END_REF], and works for all valid values of c. It does not provide an equal coverage of the unit disk, so each sample must be multiplied by φ max /π.

Validation with measured materials

To test the validity of our model, we have used it to fit all 100 materials in the MERL database [START_REF] Matusik | A data-driven reflectance model[END_REF]. We searched for optimal parameters using Levenberg-Marquadt optimisation, as implemented by Lourakis [START_REF] Lourakis | levmar: Levenberg-marquardt nonlinear least squares algorithms in C/C++[END_REF]. We conduct our fitting over the entire parameter domain (θ i , θ o , φ o -φ i ), using an L 2 norm:

E = θ i Ω measured(i, o) -predicted(i, o) 2 dω o cos θ i dθ i (64)
We restrict the fitting to θ i > 1.4, θ o > 1.4, to avoid grazing angles, where measured data is unreliable. We fit over all three channels simultaneously. We enforce physical consistency by restricting these parameters: c > 1, |η| > 1.

Parameters

All our models have 5 parameters related to the surface micro-geometry: the surface roughness σ s , 2 parameters for the normal distribution β and p and 2 parameters for the angle spread function b and c. Two of these parameters correspond to distances in the diffraction micro-geometry: b is its characteristic length, σ s is the variance of its height. As these parameters only appear in fractions with the wavelength λ: / ¯λ, σ s /λ, they are expressed in the same unit as λ (micro-meters in our implementation, to avoid numerical issues).

Our single layer model has 6 additional parameters, corresponding to the index of refraction η = n + ik of the material for each color channel, with real and imaginary parts. For our multi-layer models, there are only 4 extra parameters: the albedo or diffuse color of the underlying material, and the index of refraction η. For the multi-layer models, we made the assumption that the index of refraction would be a constant, independent of the wavelength. Index of refraction for transparent materials in the visible domain is well approximated using Cauchy's formula: η = B + C/λ 2 . In practice, the variations with wavelength appear to be negligible compared with other sources of imprecision.

Our single layer model thus has 11 total parameters, the multi-layer models have 9. For a more compact representation, it is possible to remove the imaginary part k of the index of refraction, as its influence over material behaviour is relatively small. This would bring the total number of parameters to 8 for the single layer model, at the expense of physical consistency.

Our model is extremely constrained: the micro-facet and diffraction lobes are connected by the surface roughness and the index of refraction. The number of degrees of freedom is close to a minimum. By enforcing the connection between the different components, we sometimes make it harder for the fitting algorithm to converge.

The supplemental material includes a detailed comparison between measured data, our model and two state-of-the-art material fits: Shifted Gamma Distribution [START_REF] Mahdi | Accurate fitting of measured reflectances using a Shifted Gamma micro-facet distribution[END_REF] and Smooth reflectance model [START_REF] Löw | Brdf models for accurate and efficient rendering of glossy surfaces[END_REF]. We provide lobe shape visualization, both in the incident plane and outside, rendering comparisons and error curves. Both models we are comparing against are not physically consistent: the normal distribution in [START_REF] Mahdi | Accurate fitting of measured reflectances using a Shifted Gamma micro-facet distribution[END_REF] changes with the color channel, while [START_REF] Löw | Brdf models for accurate and efficient rendering of glossy surfaces[END_REF] models diffraction effects without wavelength dependency. We get similar results while being physically consistent.

Please note that Bagher et al. [START_REF] Mahdi | Accurate fitting of measured reflectances using a Shifted Gamma micro-facet distribution[END_REF] use a total of 18 parameters, significantly more than our model; Low et al. [START_REF] Löw | Brdf models for accurate and efficient rendering of glossy surfaces[END_REF] use 9 parameters, but do not use the imaginary part of the index of refraction. They also have a diffuse lobe, with different color. Differences in fitting quality are not related to our model having more parameters. Holzschuch and Pacanowski [START_REF] Holzschuch | A physically accurate reflectance model combining reflection and diffraction[END_REF] provide a much better fit with measured materials, but use 15 parameters and separate the materials for the specular lobe and diffraction lobe.

Implementation Details

Two-Step Fitting Process

The time it takes to run an optimisation process is correlated with the number of function evaluations: fitting over the entire domain (θ i , θ o , φ o -φ i ) requires O(n 3 ) function evaluations. To speed-up the fitting, we designed a two-steps fitting process: in a first step, we fit only over the backscattering data, with i = o. As pointed out by Ashikmin and Premože [START_REF] Ashikhmin | Distribution-based BRDFs[END_REF], this slice of the data contains most of the information about BRDF behaviour. This first step is relatively quick as we only fit over a 1-dimension slice of data.

In a second step, we fit over the entire data, using as starting points the values we found in the first step. This second step converges in a small number of iterations (less than 20), as we are already close to a minimum.

We also tried a single step fitting process. It converges for some materials but not all, and tends to leave some parameters unaffected.

Color Channels Wavelengths

The most common wavelengths for the red, green and blue channels in color displays are 645 nm, 526 nm and 444 nm. We used these values in all our computations, with positive results.

The BRDFs in the MERL database were captured using a QImaging Retiga 1300 camera [START_REF] Matusik | A data-driven reflectance model[END_REF]. This camera has continuous color sensitivity for each channel, with peaks at roughly 610 nm, 540 nm and 470 nm. The raw numerical data stored in the MERL database is scaled before use, with coefficients of 1 for the red channel, 1.15 for the green channel and 1.66 for the blue channel. We found that we also obtain good fits with measured data if we use the peak sensitivity wavelengths of the Retiga 1300 camera and the unscaled data. We used the standard wavelength and the scaled data for compatibility with other systems.

Renormalization

The renormalization step is costly, but important for accuracy outside of normal incidence. We have precomputed the value of σ 2 rel for all values of b, c and sin θ i and stored them in a large look-up table, which we access at runtime.

Analysis of the Fitting

The fitting results bring the following observations:

• Materials named after semi-precious stones (aventurine, alumina-oxide, white-marble. . . ) are well approximated with the subsurface model.

• Materials named after specular plastics (specular-blue-phenolic, blue-acrylic. . . ) are well approximated with the plastic model.

• Diffraction effects are visually more important for smooth surfaces, such as chrome, nickel, specular-blue-phenolic. . .

• Several materials previously regarded as multi-lobes are well approximated by our model: red-metallic-paint, alumbronze, two-layer-silver and two-layer-gold.

• Our model does not perform so well for very rough and diffuse surfaces, including diffuse paints and fabrics. This may be due to two different cause: other effects are present including backscattering and multiple scattering inside the micro-geometry, and the fitting process has difficulties separating between the different phenomena.

• An interesting effect appears for rough surfaces (where σ > 0.5λ). For these materials, A is negligible for almost all directions (see Figure 8). The impact of the micro-facet lobe is only felt at grazing angles. As we restricted fitting to angles smaller than 1.4 radians, we do not get any input on the micro-facet distribution from the fitting. This results in surprising values for the micro-facet distribution, including highly specular lobes on some fabrics. The only solution here is to fit over the entire interval. The fitting performs better on materials where the effects of both lobes are being felt over the entire domain, such as alum-bronze.

Editing Material Parameters

Our reflectance model introduces a separation between two shape components: diffraction and reflection. Figure 10 shows the effects of editing each parameter independently:

• Reducing σ s reduces the importance of the diffraction lobe. With a low value, we get only the specular reflection. The color shift between the diffraction and micro-facet lobes appears clearly on Figure 10(b). This color shift is also present in Figure 10(a): the material is bluer at normal incidence, more orange at grazing angles, but as the transition is smooth, the effect is less perceptible.

• The presence of a diffraction lobe, combined with the reflection, explains this variation of lobe width. The diffraction lobe intensity increases as the wavelength decreases (from red to blue).

Importance of the Diffraction Component

Material Features

We visualise material response by plotting BRDF values over backscattering data (where i = o), as a function of θ h . This provides a good visual summary of each BRDF. We visualize this data using a logarithmic scale for the y axis, as BRDF values typically span several orders of magnitude (from 10 3 to 10 -4 for tungsten-carbide).

For smooth materials, such as tungsten-carbide (see Figure 11), we see two distinct behaviours: a sharp peak near the origin, a smoother tail after the peak. Our model explains these two behaviours as two different physical phenomena: reflection for the peak, diffraction for the tail, and computes the relevant parameters.

The visual impact of the diffraction lobe depends on the surface properties. For smooth surfaces, including metals, it explains most of the behaviour. The reflection lobe only becomes active at grazing angles. This is the observed behaviour for nickel (see Figure 1, top row), but also chrome and most smooth plastics. For rougher surfaces, the two lobes contribute together, sometimes with very different effects. That is the behaviour for alum-bronze (see Figure 1, middle row). Finally, for other materials, the impact is less visible, and most of the behaviour is well approximated with a Cook-Torrance lobe. That is the case for green-metallic-paint2 (see Figure 1, bottom row).

Lobe Width Variation with Wavelength

Rougher materials present a more interesting case. Measured nickel, for example, appears to have a single behaviour, with monotonic decrease (see Figure 12(a)). Renormalizing each lobe with the value at origin, we see that the lobes for the red, green and blue channels have different widths, in contradiction with the micro-facet model . Butler et al. [BNM15c] observed a similar variation in lobe width in other measured reflectance data.

Our model explains this variation of the lobe width with wavelength by the presence of a large diffraction lobe, combined with the micro-facet reflection lobe (see Figure 12(b), (c) and(d)). Summing these two lobes gives the appearance of a single lobe with varying width.

Efficient Fitting

The presence of the diffraction component also helps when fitting measured materials to the model. When we fit material parameters without a diffraction component, we get two possible behaviours: either the fitting process ignores wide angle scattering effects and focuses on the specular peak, or it tries to fit the wide angle scattering effects with the micro-facet model, with unsatisfying results (see Figure 13).

We used here our fitting algorithm, working over the entire data (θ i , θ o , φ o -φ i ). Other fitting approaches [AP07, BSH12, DHI * 15] avoided the issue by explicitly fitting only over θ h with the back-scattering data. By construction, this gives a micro-facet model that approximates the BRDF well, but it does not take into account data outside of the incident plane.

Behavior at Grazing Angles

Several papers [BSH12, [START_REF] Durand | Experimental analysis of BRDF models[END_REF][START_REF] Burley | Physically-based shading at Disney[END_REF] have noted the discrepancy between the behaviour of measured materials and the behaviour predicted by the Cook-Torrance model at grazing angles. Two components of the Cook-Torrance model rise at grazing angles: the Fesnel term F(θ d ) and the 1/(cos θ i cos θ o ) term. For rough surfaces, the latter is partially counterbalanced by the shadowing and masking term G. For smooth surfaces, we have G ≈ 1 for almost all directions. Figure 14 shows the measured data for chrome, as a function of θ d , along with a fit using only Cook-Torrance model and a fit using our model, a combination of diffraction and Cook-Torrance. For this material, the fitting algorithm predicted that most of the energy was going into the diffraction lobe. The Cook-Torrance approximation begins rising quickly for θ d > π/4. For comparison, we plot 1/cos(θ d ) 2 . The two curves are parallel, showing that most of the rise comes from the 1/(cos θ i cos θ o ) term. The Fresnel term explains the difference between these two curves. By comparison, our model rises later, around θ d = 1.4. For most of the interval, A ≈ 0, and the diffraction lobe explains most of the behaviour. It is only for grazing angles that A ≈ 1 and we revert to the Cook-Torrance model.

Multiple Scattering

Our model is limited to a single light interaction with the surface: both reflection and diffraction are modelled by a single bounce. This approximation is valid for smooth materials: metals, shiny plastics, polished stones, but breaks for rougher materials, such as glossy paints. Reflectance values for a glossy paint such as green-metallic-paint show a large lobe, followed by a plateau for θ h > π/4 (see Figure 15). Our model is unable to capture this behaviour, although our approximation is visually acceptable.

Conclusion and Future Work

We have presented a physical explanation of the different phenomena involved in material reflection. Based on the Generalized-Harvey-Shack theory, our model combines together diffraction and reflection caused by surface micro-geometry: the reflection lobe is responsible for the specular peak of the reflectance, while diffraction is responsible for wide-angle scattering. The energy repartition between the different lobes depends on the surface roughness, which is defined as the variance of surface height variations. We validated our model by comparing it with measured materials.

The relative importance of the diffraction lobe depends on surface roughness: it can be an important part of the surface aspect for semi-rough materials. Our analysis explains many observed behaviours of measured reflectance, such as the difference between the peak and the tail, and the size of the lobe changing with the wavelength. This new explanation of the influence of surface micro-geometry on material reflectance opens the way for new research on material acquisition and representation. In future work, we want to exploit the knowledge provided by this new model to speed-up material acquisition, possibly using polarized light.

Our model is only based on physical parameters: the index of refraction of the materials, the albedo of participating media. These parameters are useful for physical explanation of surface properties, but not practical for material editing by artists. In the future, we would like to connect these parameters to more convenient parameters for artists, such as color or roughness.

Finally, We have focused on a single light interaction with the surface. This is consistent for relatively smooth surfaces, but rougher materials break this assumption. For these materials, we would like to model multiple interactions with the surface, as in [START_REF] D'eon | A comprehensive framework for rendering layered materials[END_REF][START_REF] Heitz | Multiple-scattering microfacet BSDFs with the Smith model[END_REF] but taking into account the color change caused by diffraction.

Figure 3 :

 3 Figure3: Light reflected off the surface at a height h has travelled a shorter distance. This Optical Path Difference (OPD) causes interference between parallel rays and is responsible for diffraction effects. Here, OPD = (cos θ i + cos θ o )h(x, y)[START_REF] Krywonos | Predicting Surface Scatter using a Linear Systems Formulation of Non-Paraxial Scalar Diffraction[END_REF]. Heights h are relative to the average surface level.

Figure 5 :Figure 6 :

 56 Figure5: In the plastic two-layer model, we assume we have a transparent layer of refractive varnish over a diffuse material. We model light bouncing several times between the two interfaces before leaving the object.

Figure 7 :

 7 Figure 7: Importance sampling reduces the noise in global illumination computations (1000 samples per pixel in both pictures).

Figure 8 :

 8 Figure 8: The A coefficient express the ratio of energy going into the reflection and diffraction lobes. It depends on the incoming angle and the σ s /λ ratio.

Figure 10 :

 10 Figure 10: Effects of parameter edition for micro-facet geometry for our reflectance model. (a) nickel, with parameters from fitting. (b), (c) and (d): effect of modifying each of the parameters, in a grid pattern. We used different parameters values for the squares and the borders of the grid. For (d), σ s = 0.05µm to increase the impact of the micro-facet lobe.

Figure 12 :

 12 Figure 12: (a) Measured reflectance data for nickel as a function of θ h , renormalized by the value at the origin. The width of the BRDF lobe changes with the wavelength, in contradiction with the micro-facet theory. (b) (c) and (d): The presence of a diffraction lobe, combined with the reflection, explains this variation of lobe width. The diffraction lobe intensity increases as the wavelength decreases (from red to blue).

Figure 13 :

 13 Figure13: We tried fitting measured materials without a diffraction lobe. For steel, left, the fitting process ignores wideangle scattering, but provides satisfying results. For brass, right, the fitting process tries approximating wide-angle scattering with a larger micro-facet distribution, resulting in unsatisfying results.

Figure 15 :

 15 Figure 15: Our model is limited to a single interaction between the light and the material. It has problems approximating more glossy materials, such as green-metallic-paint, probably because it does not model multiple light bounces inside the micro-geometry.

  Figure11: Measured reflectance data for tungsten-carbide as a function of θ h and our fitting. This reflectance has two separate behaviours: a sharp peak and a softer tail. Our model provides a simple explanation: these correspond to different physical phenomena, micro-facet reflection and diffraction.

	���� ��			������������ �������������� �������������		0.1 1			Red (measured) Together Cook-Torrance Diffraction		0.1 1			Green (measured) Together Cook-Torrance Diffraction		0.1 1			Blue (measured) Together Cook-Torrance Diffraction
	�����						0.01						0.01						0.01
	������						0.001						0.001						0.001
	�������	�	���	���	����	���	0.0001	0	π/8	π/4	3π/8	π/2	0.0001	0	π/8	π/4	3π/8	π/2	0.0001	0	π/8	π/4	3π/8	π/2
				� �						θ h						θ h						θ h
				1 10 100 1000								Red (measured) Together Cook-Torrance Diffraction
					0.1
				0.001
				0.0001
			0.00001 0					π/8			π/4			3π/8				π/2
													θ h