Data intrinsic correction for working distance variations in MA-XRF of historical paintings based on the Ar signal

The nonplanar shape of a painting as well as practical constraints often result in the painting's surface not being parallel to the plane in that the measurement head of a MA-XRF scanner is being moved. Changing the working distance affects the measurement geometry, so that the sensitivity for the same element may vary throughout the investigated area and induce visible artifacts. These artifacts are especially visible when different scans of the same painting are stitched together. In this article, we present an approach to correct for the variation of the measurement distance. We explored using an intrinsic part of the XRF data set, the Ar signal from the air, to estimate the distance between surface and instrument. The model is developed based on fundamental parameter calculations and a measurement of a NIST 610 standard and is verified on a set of scans of Rembrandt's ‘ Portrait of Oopjen Coppit (1611 – 1689) ’ .

The nonplanar shape of a painting as well as practical constraints often result in the painting's surface not being parallel to the plane in that the measurement head of a MA-XRF scanner is being moved. Changing the working distance affects the measurement geometry, so that the sensitivity for the same element may vary throughout the investigated area and induce visible artifacts.
These artifacts are especially visible when different scans of the same painting are stitched together. In this article, we present an approach to correct for the variation of the measurement distance. We explored using an intrinsic part of the XRF data set, the Ar signal from the air, to estimate the distance between surface and instrument. The model is developed based on fundamental parameter calculations and a measurement of a NIST 610 standard and is verified on a set of scans of Rembrandt's 'Portrait of Oopjen Coppit (1611-1689)'.

| INTRODUCTION
Macroscopic X-Ray Fluorescence imaging (MA-XRF) has been established in the last years for the investigation of historical paintings, providing representative information about the pigment use and to reveal later changes to a paintings composition. In some cases it was possible to revisualize abandoned and overpainted compositions that are not accessible via other techniques. [1,2] The information obtained from investigating a sample by any spectroscopic technique can be separated into the signals resulting from the sample and artifacts that are contributions from the instrument and the surrounding of the experiment. The distribution of pigments in a painting is very heterogeneous and the contribution of many artifacts is commonly negligible compared to the variance induced by the painting's heterogeneity. This allows seeing paintings during the analysis of XRF data often as perfectly flat objects and the measurement geometry for every point as identical.
However, the plane in that the measurement head of the scanner is moved and the plane of the painting's surface are not always perfectly parallel. While small variations are negligible, large ones of a few millimeters can contribute to visible variations in sensitivity for certain elements and result in visible artifacts, especially when different subscans of a painting are stitched together, as shown below.
The fact that instrument and painting are not well aligned is not the result of carelessness of the operators. Most historical paintings are slightly bend, so that a measurement head fixed to a two dimensional plane cannot follow the surface. Further, moving large historical paintings with millimeter precision is not always possible outside of dedicated workshops, so that aligning them to the scanner is difficult. Finally, it is a general practice to have a 'Safety First' approach to measurements in that scan patterns are chosen that produce elemental distribution images with stronger artifacts while minimizing the risk of collision with the painting. Only few XRF scanners allow following the surface of a painting, such as the LANDIS scanner making use of a laser distance sensor. [3] Other instruments described are limited to moving the measurement head in a two dimensional plane. [4][5][6][7][8] One of the reasons for this is that automatic distance correction requires perfect implementation beyond doubt, in order to not contradict the 'Safety First' approach.
The effect of an uneven surface in MA-XRF is well known, albeit until now mostly the effect of sample curvature and thus changing excitation and detection angles is discussed. Consequently, the main effect considered is self-absorption in a paint layer. This problem is of special interest when studying painted marble from antiquity. Approaches described so far used the difference of the signals of two detectors, [9] the variance of the Ca signal in a nearly pure marble sample induced by topography [10] or the 3D information obtained by photogrammetry. [11] The last study included a correction for the absorption for X-rays in air.
The approaches so far either required external information, such as photogrammetry or laser based distance measurements, or made simplifying assumptions about the sample, that is, assuming a thin layer of paint on homogeneous marble.
To estimate and correct for the distance between measurement head and painting surface we will use an intrinsic part of the data acquired by MA-XRF: the Ar signal resulting from the air.
In order to make our approach easily applicable by other groups, we combined a limited analytical model, including solid angle and quantum efficiency of the detector as well as absorption in the air and sample and fit it to measurements of a reference material. It would be feasible to build an entire fundamental parameter (or Monte Carlo) simulation of the measurement geometry, but this would require high effort in modeling and characterization of the instrument, especially its polycapillary optic. Acquiring correction factors for each element separately would require comparable effort to the approach described here but be highly specific for the instrument used.
In what follows we will first describe the data acquisition of reference material and verification data, followed by the theoretical model used before demonstrating the validity of the calculated correction factors in a real application.

| DATA ACQUISITION AND EVALUATION
For the development of the model a number of spot measurements at different distances were done on a standard reference material (SRM) of the National Institute of Standards and Technology (NIST) 610 'Trace Elements in Glass'. The NIST 610 is a 1 mm thick blue disk of Ca rich glass, which contains the majority of the elements of the periodic table at a nominal concentration level of 500 ppm. [12] In order to test the model and verify its effectiveness data of a historical painting acquired independent from the measurements of the NIST 610 was needed. A previously acquired data set from Rembrandt's 'Portrait of Oopjen Coppit' from 1634 (inv. No. SK-C-1768) was selected as in here artifacts resulting from the variation of the working distance were clearly discernible. The 210.5 × 134.7 cm sized oil on canvas painting was investigated before a conservation treatment after the purchase of it and its pendant 'Portrait of Marten Soolmans' by the Rijksmuseum (Amsterdam, the Netherlands) and the Louvre (Paris, France). [13] Both objects were investigated with the prototype of the Bruker M6 Jetstream, which was characterized previously in detail. [4] The measurement head consists of a 30 W Rh anode X-ray tube with a polycapillary optic, a silicon drift detector (SDD) and two cameras for sample observation. The focus of the polycapillary is approximately 5 mm in front or the measurement head and has an (energy dependent) focus diameter of less than 50 μm. However, the common measurement distance is 10 mm where the beam has widened to around 300 μm. The measurement head is mounted on a 80 × 60 cm 2 motorized stage, which is used to scan the surface of a painting. A third motorized stage allows changing the distance of measurement head and surface between scans.
The NIST 610 reference material was brought in physical contact with the protective casing of the measurement head. In this position a spectrum was acquired and the measurement head moved backwards in steps of a millimeter and further spectra were acquired up to a distance of 50 mm. Each spectrum was acquired for 60 s with X-ray tube settings of 50 kV and 600 μA. The surface of the painting was scanned with the same tube settings, a step size of 700 μm and a dwell time of 70 ms. Nine scans were needed to scan the entirety of the painting of which four are used to verify our approach.
The spectra of the NIST 610 were processed in PyMCA [14] and normalized to live time by in-house written Python scripts. The data of the painting were processed using PyMCA and datamuncher [15] and also normalized to the live time per pixel, a feature added to the latter for this study.

| MODELING
A correction for the working distance requires a measure for it, for which we will use the Ar signal. The primary radiation does not pass through the air without interaction. It is scattered and excites the air to emit fluorescence radiation on its own. The scattered Bremsstrahlung continuum of the X-ray tube contributes to the background in the acquired spectra but of the elements of the air only Ar has sufficiently high characteristic radiation and abundance in order to be detected. Both signals are considered artifacts, that is, not belonging to the sample, and thus the detector is commonly collimated to maximize its solid angle for the designed distance of analysis. Further, it covers the edge of the detector's active area where a partial loss of the charge generated by the photoelectric effect might occur and result in stronger tailing of peaks. Hence, the amount of Ar in the beam is proportional to, but not linearly dependent of, the working distance.
Like any other spectrometer the M6 consists of a radiation source, sample and detector, as shown in Figure 1. The primary radiation is focused (and filtered) by the polycapillary optic, passes the air, interacts with the sample before being recorded by the detector after passing through the air again. As we are searching for relative correction factors the absolute intensity of the primary radiation is not relevant for these calculations. The recorded intensity of element i, dependent on the working distance is given in Equation (1).
Ω(d) is the solid angle from that the detector is recording the fluorescence radiation emitted by the beam matter interaction. A i (d) is the absorption of the fluorescence radiation in the air and is described in Equation (2). ϕ i (d) is the quantum efficiency of the detector. It is explained in Equation (3). k i is an element specific constant used to normalize the curves.
To calculate the solid angle a common approximation is dividing the active area of the detector by the surface of a sphere with the radius of the detector-sample distance. This, however, does not take the collimator that reduces the air scatter, and the tilt of the detector into account, which is not straightforward. It requires calculating the angle between vectors from the edge of the active area of the detector that is exposed to the fluorescence radiation to the origin of the fluorescence radiation. This calculation includes next to the working distance also a number of geometric model parameters that are shown in Figure 1d and summarized in Table 1. These are the radius and height of the collimator (r c , h c ), the offset of the detector from the point where the primary radiation leaves the protective casing of the measurement head in horizontal and vertical direction (o v , o h ) and the tilt of the detector relative to the primary beam (α m ). In these calculations, we assume the active area to have the shape of an ellipse, which is a simplification toward its actual shape, described in Wrobel et al. [16] The geometric operations performed to calculate the solid angle are given as DSI.
In Figure   shown. For Zn(−K α : 8.63 keV) dashed and dotted line are nearly identical and describe the recorded intensity well. This is not true for Si(−K α : 1.74 keV) and Ca(−K α : 3.69 keV), as the less energetic fluorescence radiation is much stronger absorbed in the air than that of Zn. The absorption can be modeled by Equation (2).
d is the distance from polycapillary exit window to beam/sample interaction, d eff is the distance from beam/ sample interaction to detector, ρ is the density of air, E pi is the energy of the primary beam, E fi is the energy of the fluorescence radiation, μ is the mass absorption coefficient, dependent on these energies. Assuming dry air near sea level, as defined by NIST, these variables can be easily obtained from xraylib. [17] Instead of calculating the entire primary spectrum emitted from the X-ray tube and weighting it by the photo electric cross section of element i we used a single energy per element, which is the edge energy of the fluorescence line in question times the Edge Factor given in Table 1. A similar approach was described earlier. [18] This correction allows to also properly model the distant dependent intensities of Si and Ca, as shown in Figure 2.
While Si, Ca and Zn are well modeled by this approach, Zr(−K α : 15.75 keV) is not, albeit its fluorescence line has a considerably higher energy than that of the other elements and it should be less affected by absorption. This is true, but we are observing a different effect: the distance dependent quantum efficiency of the detector. In order to detect an X-ray photon, it needs to be absorbed within the active volume of the detector, which is comparably thin (initial estimate 450 μm) and absorbs only 65% of the incoming Zr-K lines. If the detector is not perpendicular to the sample's surface, as it is closer than the design distance, the path of the incoming X-rays is longer in the active volume and thus the quantum efficiency raises and even overcompensates for the loss in solid angle. This is, however, more a curiosity than an exploitable feature as the 'gain' from tilting the detector is less than 5%. The absorption of radiation in the detector's Be window is also included in the curve. d Si is the thickness of the detectors active volume, ρ Si is the density of Si, E fi is the energy of the fluorescence radiation, α is the angle between the surface of the detector and the incoming fluorescence radiation from the sample. d Si is a parameter of the model, given in Table 1. ρ Si and E fi are fundamental parameters, while α can be calculated from the geometric parameters in Table 1 and the working distance. The more energetic radiation of Zr can escape from deeper inside the NIST 610 and reach the detector than that of lighter elements and thus the solid angle of the detector should not be calculated from the surface of the sample/beam interaction but deeper.
If one looks at the information depth for Zr in the NIST 610, the depth that contributes together 99% of the recorded signal, one obtains 1,210 μm, actually behind the 1 mm thick sample. However, 50% of the recorded radiation is emitted from the first 180 μm, so that it is not responsible for the shift of the Zr line, which is 1.7 mm compared to the maximum in the solid angle curve. However, this distance is included in the results shown below by shifting the curves accordingly.
As the model allows satisfyingly describing the relationship between distance and recorded intensity for most elements in the NIST 610, it is also possible to integrate the air path to the sample and describe the curve of Ar (see Figure 3). Here, the calculated curve was normalized to the recorded curve to have the correct absolute intensity of Ar. From this Ar curve and the predicted curves for the other elements one can calculate correction factors. These work well for the range between 0 and 20 mm, so well around the design distance of 10 mm. The model is neglecting the fact that fluorescence and scattered primary radiation also can excite the Ar, which means that the Ar intensity is, especially for long distances slightly underestimated. But as the measured Ar intensity varies only slightly with the intensity of the emitted fluorescence radiation, as can be seen in the DSI, this limitation was considered acceptable. After developing the model and selecting suitable initial values for all parameters we used a numerical optimizer to find the predictions made by our model to the measured curves. Not all recorded intensity curves follow the model as well as those shown in Figure 2, most contain considerable statistical noise. Fitting the model to all curves simultaneously allows using the information from other elements to improve the statistics of individual correction curves.

| VERIFICATION
In Figure 4, a part of Rembrandt's 'Portrait of Oopjen Coppit (1611-1689)' is shown together with its corresponding Ca and Fe distribution images. As the painting is larger than the range of the motorized stage of the Bruker M6 it was scanned in several parts that were later stitched together. The upper right scan 'B' employed due to geometric constrains in its lower left corner a shorter working distance than the other scans, resulting in a reduced recorded Ca intensity. To the right the same area after application of the correction factor is shown in Figure 3, which results in the section blending in much better. The same can be observed for the Fe image.
For applying above calculated correction factors, the Ar signal for each pixel needs to be calculated. Give that each pixel was acquired with a dwell time of 70 ms the statistical noise is significant. Adding this noise to the image is not desirable. As the surface of a painting is F I G U R E 3 Ar intensity as measurement (triangles) and simulated (dashed line) and values and correction factors derived from the simulation comparably smooth one can assume neighboring pixels to have similar signals, but one also needs to take the edge of a painting properly into account where sudden changes of Ar intensity can be observed. The best compromise found was a two dimensional Savitzky-Golay filter that reduced the noise and preserved the edges. In two homogeneous areas of the scans A and B this allowed reducing the standard variation of the Ar signal considerably (Scan A: 242 ± 87 ! 242 ± 21 cps, Scan B: 155 ± 79 ! 155 ± 30 cps). The Ar intensity maps of the scans A and B are given as DSI.
To go beyond a purely visual comparison of images, the overlapping regions of the upper left scan (A) and the upper right scan (B) were extracted and histograms of their Ca signals were calculated. As it can be seen in Figure 5, the Ca signal of A is indeed corrected to the F I G U R E 5 Area of joint with histograms of Ca distribution, before and after correction by Ar signal but including in both cases dwell time correction same average intensity as B. However, both peaks were slightly broadened as the statistical noise of the Ar signal, albeit reduced, is added to that of the Ca.
While the intensities are well adjusted, upon close visual inspection the upper right square in the elemental distribution images in Figure 4 is a bit off. To explain this, the area of the sitter's proper left shoulder was magnified in the same figure. The upper part of the image is much sharper than the lower part. This is due to the fact that in the upper scan the instrument/painting distance was around 5 mm, while in the lower part it was around 10 mm, consequently the beam in the upper part was of the size of 100 μm, a factor of 3-5 smaller than in the lower scan. Thus, the upper scan is under-sampled and the resulting sharper image contains some contribution from the canvas structure without really resolving it, while in the lower scan these are averaged out. A similar effect of under-sampling and canvas was observed earlier in a painting by Vincent van Gogh. [19]

| CONCLUSION
We have shown that it is possible to correct for variations in the working distance by making use of the Ar signal intrinsic to the acquired MA-XRF data. The approach uses a simple model that can be easily calibrated, making use of a data set of the NIST 610, which can be acquired in less than 2 hr. Consequently, it is possible to also correct already acquired data without revisiting the object if the suitable calibration can be obtained.
While suitable for the presented case, the model used has a few limitations. It does not take the effect of a strongly tilted surface, especially for self-absorption, into account, that previous studies observed to be significant in strongly curved objects. [9][10][11] It needs to be explored how well the Ar signal can be used to follow a complex contour and what the limits in resolution of the method are. Also, it is assumed that all elements are in a same homogeneous layer in the sample that is perpendicular to the incoming beam and that no absorption effects are enhanced by changing the detection angle.
Clearly beyond the focus of this study is to smooth out the differences between over-and under-sampled images and to use the 3D height information obtained to project the images into the same plane for stitching. The routines used for correcting these images are available upon reasonable request to the first author.