Terahertz frequency modulated continuous wave imaging advanced data processing for art painting analysis

Reflection terahertz frequency modulated continuous waves scanner (300 GHz) has been proficiently optimized for imaging two easel paintings of different age. The information content of the obtained THz images has been fully inspected by selecting the appropriate THz image parameter. At the same time, a new data processing has been developed for improving the level of detail held by the axial parametric THz images by means of Gaussian fit of the reflected signals. By carefully weighting the reflected signals as a function of the optical path, the reflected amplitude has been corrected for the positioning of the object surface with respect to the beam focal point. The artifact affecting the THz images recorded from an uneven painting surface have been resolved and the obtained images fairly represent to the original painting. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (110.6795) Terahertz imaging; (100.2980) Image enhancement. References and links 1. K. Fukunaga, THz Technology Applied to Cultural Heritage in Practice (Springer, 2016). 2. W. Köhler, M. Panzer, U. Klotzach, and S. Winner, “Non-destructive investigation of paintings with THzradiation,” in 9th ECNDT (2006), pp. 1–7. 3. D. Mittelman, Sensing with Terahertz Radiation (Springer Berlin Heidelberg, 2003). 4. J. Bianca Jackson, R.M. Mourou, J. F. Whitaker, I. N. Duling, S. L. Williamson, M. Menu, and G. Mourou, “Terahertz Time-Domain Reflectometry Applied to the Investigation of Hidden Mural Paintings,” in Conf. Lasers Electro-Optics/Quantum Electron. Laser Sci. Conf. Photonic Appl. Syst. Technol. (2008), paper CThN3. 5. C. L. Koch-Dandolo, T. Filtenborg, K. Fukunaga, J. Skou-Hansen, and P. U. Jepsen, “Reflection terahertz timedomain imaging for analysis of an 18th century neoclassical easel painting,” Appl. Opt. 54(16), 5123–5129 (2015). 6. A. Younus, J. P. Caumes, S. Salort, B. Chassagne, C. Pradère, A. Dautant, A. Ziéglé, and E. Abraham, “A continuous millimeter-wave imaging scanner for art conservation science,” Adv. Opt. Technol. 2011, 1–9 (2011). 7. H. Zhang, S. Sfarra, K. Saluja, J. Peeters, J. Fleuret, Y. Duan, H. Fernandes, N. Avdelidis, C. Ibarra-Castanedo, and X. Maldague, “Non-destructive Investigation of Paintings on Canvas by Continuous Wave Terahertz Imaging and Flash Thermography,” J. Nondestruct. Eval. 36(2), 34 (2017). 8. L. Öhrström, B. M. Fischer, A. Bitzer, J. Wallauer, M. Walther, and F. Rühli, “Terahertz Imaging Modalities of Ancient Egyptian Mummified Objects and of a Naturally Mummified Rat,” Anat. Rec. (Hoboken) 298(6), 1135– 1143 (2015). 9. J. Guillet, K. Wang, M. Roux, F. Fauquet, F. Darracq, and P. Mounaix, “Frequency Modulated Continuous Wave Terahertz Imaging For Art Restoration,” in Proceedings of IEEE Conference on Infrared, Millimeter, and Terahertz Waves (IEEE, 2016), pp. 1–2. 10. J.-P. Guillet, M. Roux, K. Wang, X. Ma, F. Fauquet, H. Balacey, B. Recur, F. Darracq, and P. Mounaix, “Art Painting Diagnostic Before Restoration with Terahertz and Millimeter Waves,” J. Infrared Millim. Terahertz Waves 38(4), 369–379 (2017). 11. D. E. Barriok, “FM/CW radar signal and digital processing,” NOAA Technical Reports ERL283–WPL20, (1973) 12. C. A. Weg, W. V. Spiegel, R. Henneberger, R. Zimmermann, T. Loeffler, and H. G. Roskos, “Fast Active THz Vol. 26, No. 5 | 5 Mar 2018 | OPTICS EXPRESS 5358 #313724 https://doi.org/10.1364/OE.26.005358 Journal © 2018 Received 17 Nov 2017; revised 22 Dec 2017; accepted 23 Dec 2017; published 22 Feb 2018 Cameras with Ranging Capabilities,” J. Infrared Millim. Terahertz Waves 30(12), 1281–1296 (2009). 13. J.-J. Lin, Y.-P. Li, W.-C. Hsu, and T.-S. Lee, “Design of an FMCW radar baseband signal processing system for automotive application,” Springerplus 5(1), 42 (2016). 14. N. Palka and D. Miedzinska, “Detailed non-destructive evaluation of UHMWPE composites in the terahertz range,” Opt. Quantum Electron. 46(4), 515–525 (2014). 15. E. Cristofani, F. Friederich, S. Wohnsiedler, C. Matheis, J. Jonuscheit, M. Vandewal, and R. Beigang, “Nondestructive testing potential evaluation of a terahertz frequency-modulated continuous-wave imager for composite materials inspection,” Opt. Eng. 53(3), 31211 (2014). 16. J. Dong, A. Locquet, M. Melis, and D. S. Citrin, “Global mapping of stratigraphy of an old-master painting using sparsity-based terahertz reflectometry,” Sci. Rep. 7(1), 15098 (2017). 17. M. Stecher, et al., “14. Towards Industrial Inspection with THz Systems,” in Ultrashort Pulse Laser Technology, S. Nolte, F. Schrempel, F. Dausinger, eds (Springer, Cham, 2016) 18. H. Quast and T. Loffler, “3D-terahertz-tomography for material inspection and security,” in Proceedings of IEEE Conference on Infrared, Millimeter, and Terahertz Waves (IEEE, 2009), pp. 1–2. 19. K. Zuiderveld, “Contrast Limited Adaptive Histogram Equalization,” in Graphics Gems IV (Elsevier Inc., 1994),


Introduction
Terahertz radiation (0.1 -10 THz frequency range, 1 THz = 10 12 Hz) has been used for analyzing cultural heritage objects since 2006 [1,2]. The most common devices used for this purpose are pulsed terahertz Time-Domain systems (TD THz). In fact, in reflection configuration, they are able to provide non-invasive cross-sections of the object under investigation by imaging the time-of-flight of the reflected THz pulses [1,[3][4][5], thus giving information about the internal structure of the object. Without being capable of providing depth information about a sample, Continuous Wave terahertz (CW THz) imaging systems have been also tested for investigating cultural heritage objects such as panel and easel painting replicas, jars, mummies (0.11 THz [6]), and lately real easel paintings (0.1 THz [7]).
Just a few years ago, artworks started to be scanned by Frequency Modulated Continuous Waves terahertz (FMCW THz).These imaging systems have been used for scanning ancient mummified samples (0.1 THz and 0.3 THz [8]) and easel paintings (0.1 THz and 0.3 THz [9], [10]). Similarly to the time-domain systems, the reflection mode FMCW THz systems allow for the depth inspection of the object under investigation thanks to the frequency sweep of the emitted radiation [10][11][12][13]. TD and FMCW THz can collect significant information about the stratigraphy of multilayered samples and their internal structure [14,15]. TD THz can be applied to the examination of thin samples (~1 cm), but with high axial resolution (minimum layer separation distance ~30-50 μm) [14][15][16]. The FMWC method can be applied to the examination of thicker samples (~5-15 cm), but its axial resolution is only about 1-2 mm [14,15].
Having the capability of providing non-invasively surface and depth information about the scanned object and being faster (up to 10 kHz acquisition rate [12] versus the ~100 Hz-1kHz of the TD system [17]) and lower priced [15] than pulsed terahertz systems, they are likely to spread more into the conservation and restoration sectors than the time-domain systems.
In this paper, we will show how to optimize the information provided by the data acquired by the FMCW THz from highly heterogeneous objects, such as easel painting, by selecting the appropriate image parameters and by means of advanced signal processing.

The investigated paintings
We have investigated two paintings on canvas. Owing to the Gaussian nature of the beam, a scan performed by a XY scanning stage requires the surface of the object under investigation to be flat for positioning the surface plane within the Rayleigh range of the beam. For this reason, paintings can be considered a good subject for being scanned with this technology. The first painting represents some boats in a quiet lake [Figs. 1(a)-1(c), between the 20th and 21st centuries, acrylic on canvas]. The second one represents a bucolic scene where a woman holding a child is surrounded by two putti and a male figure [ Fig. 1(d), 18th century, oil on canvas].

3.Materials and methods
The paintings have been scanned with a FMCW imaging system (SynViewScan TRMF, by SynView GmbH, Bad Homburg -Germany) in reflection mode at normal incidence by using the 300 GHz system head. The achievable depth resolution of the system is in the first instance limited by the bandwidth of the frequency modulation (Δf), according to dZ = (cf b )/2nΔf, where c is the speed of light in air, f b the beat signal frequency, n is the refractive index of the medium. For obtaining an effective operating frequency of 300 GHz, the frequency of the voltage controlled oscillator (VCO) changes linearly between 0.23 THz to 0.32 THz, so that the Δf is equal to 0.09 THz [15,18].
The contemporary painting has been scanned from Z = −70 mm to Z = 70 mm (141 data points), with a dZ step of 1 mm, with Z being the distance between the transceiver and the object surface, Z = 0 the distance at the beam waist and Z > 0 for distances closer to the transceiver. For the 2-D raster scan the steps sizes (horizontal, dX and vertical, dY) have been set equal to 1 mm, obtaining a 594 x 499 pixels image.
The 18th century painting has been scanned with a dZ step of 0.2 mm from Z = −25 mm to Z = 25 mm (251 data points) and dX and dY steps equal to 1 mm, obtaining a 799 x 600 pixels image. While the minimum distance between two adjacent data points dZ is equal to 1 mm for the SynViewScan system configuration file, this last has been modified for reducing the dZ step along the Z axis.
The linear scale signal amplitude has been employed for calculating the parameters used for plotting the presented THz images. When specified, contrast limited adaptive histogram equalization (CLAHE) has been implemented for enhancing the contrast of the calculated THz image [19].

Source of contrast in the FMCW amplitude image
The profile and surface roughness of the object under investigation can be evaluated from the location of the maximum of the recorded signals. On the other hand, the shape of the recorded pulse can provide us important information about the stratigraphy of the object under investigation, especially considering that the resolution bandwidth is often unable to resolve two closely spaced peaks in the recorded signal. For evaluating the shape of the pulse, it is widely used the full width at half maximum (FWHM) parameter.   By comparing it with the surface topography of the painting obtained by illuminating the object with raking light [ Fig. 1(b)] it can be observed that the black lines of the THz image correspond to the sharp edges of the thickest brushstrokes used by the artist to model the subjects in relief. Figure 3(b) displays the signals recorded at (x, y) coordinates corresponding to low [in correspondence of the sharp edges, green and red curves in Fig. 3(b)], medium [yellow and orange curves in Fig. 3 Fig. 3(b)] reflectivity values, which have been collected from neighboring pixels in two different areas of the painting. The reflectivity value of the signals recorded in correspondence of the black lines [ Fig. 3(a)] drops notably off compared to that recorded in other region of the painting, possibly because of the edges diffraction effect.

(b)] and high [blue and violet curves in
The various color areas of the painting show differences in the reflectivity values such that a reasonable gray contrast in the THz image is provided. Looking at the painting, all the depicted clouds appear of the same white color whereas they seem differently and strongly contrasted in the THz intensity record, suggesting that the artist made them with different paint materials. The canvas texture is also identifiable in the THz image of the signal global maximum: it is particularly distinguishable in the area located under the foreground vessel on the left of the painting [ Fig. 3(a), red dashed square]. An anomalous rectangle-like area is found in the top-right corner of the THz record, located as indicated by the yellow arrow in Fig. 3(a). This anomaly corresponds to the white adhesive label attached on the back of the painting [ Fig. 1(c)]. Figure 4(a) shows the image of the signal amplitude integrated between Z = −13 mm and Z = −8 mm and it shows the wooden framework of the painting. In fact, the peak arising from the radiation back-reflected at the framework is mainly located within these Z positions, just before the main peak associated to the reflection at surface. Figures 4(b)-4(c) show the signals recorded at (100,300) and (100,200) (x, y) coordinates (black and blue lines), where no framework is present on the back of the painting, and at (300,300) and (300,200) (x, y) coordinates (red and magenta lines), where the framework is present. The peak associated with the reflection at the wooden structure (highlighted by black arrows in Figs. 4(b)-4(c) is clearly visible where the framework is present underneath the canvas. The second peak appears at minor Z values respect to those of the main one because the values of the Z axis increase for distances closer to the transceiver.
When present, the second peak is distinguishable in the recorded signals since the thick wooden frame is spaced several millimeters away from the painting surface (≤ 8 mm) and an air gap is present between the canvas support and the wooden frame. Differently, the thin white adhesive label imaged in Fig. 4(a) is attached directly on the back of the canvas. The width of the main peak at the baseline is estimated ~6 mm; for this reason, additional minor peaks associated with reflections at interfaces located at a distance from the surface minor than this value cannot be distinguished as separate peaks, as it is the case of the white label. In spite of this, as previously observed, the layer thickness and material variations associated to the presence of the white label on the back of the canvas have an influence on the amplitude of the main reflected peak. This is the reason why the white label has been imaged by using the parameter global maximum of the reflected signals [see Fig. 3(a)].  Fig. 4(a) as colored dots. c. The same signals of Fig. 4(b), rescaled for enhancing the visibility of the second peak registered for the signals recorded at (x, y) coordinates at which the frame is present on the back of the painting.

Significance of the calculated image parameters versus the numeric precision of the signal amplitude values: Gaussian curve fit
The step size dZ defines the numeric precision in the determination of the position of the maximum and FWHM from the recorded back-reflected signals. Even if selecting finer dZ during data acquisition could help in resolving two closely spaced peaks from the recorded signals, it will significantly and undesirably increase the acquisition and computation times.
Nevertheless, smaller step sizes dZ can be obtained after the data recording by fitting a mathematical model to the acquired signals. Compared to those obtained from the raw recorded data, more precise values can be determined from a fitted model for the location of the signal maximum and for the FWHM.
The recorded signals (main peaks, linear scale) are suitably described by the one term Gaussian model , where the coefficients a, b, c correspond to the amplitude, the centroid (location) and the standard deviation respectively, while x is the recorded data set [Figs. 5(a)-5(b) and Table 1].   Figure 6(a) show the THz image of the location of the maximum of the recorded signals along the Z axis, from −5 to 0 mm, with a dZ step of 1 mm (i.e. the data acquisition step). It can be observed how the key information provided by this image is that the main plane of the painting is tilted from the focal plane, being the top right corner the closest and the bottom left part the farthest from the focal plane. We can also observe some differences in the image parameter owed to the surface roughness of the object, particularly evident for the outlines of the two boats in the foreground. They are located closer to the transceiver compared to the neighbor pixels, possibly meaning that the artist deposed a thicker paint layer for let them stand out. The farthest area of the painting from the transceiver is located in the top-left corner of the painting.  The key information obtained by the image of the b coefficient is no longer the main tilting of the panting from the focal plane, as previously observed [ Fig. 6(a)], but the fine surface roughness of the painting. Apart from the foreground boats, further details of the painting seem to be in relief, such as the clouds and the other boats. The presence of the label on the back of the painting had an influence on the shape of the recorded signals and it can be observed how the Gaussian b values in that region differ from that of the adjacent pixels, especially at the label edges. In addition, a periodic oscillation of the b value parallel to the Y axis of the image is found, which is possibly related to a linear motor stage misalignment. Figure 7(a) shows the image of the FWHM (mm) of the acquired signals. Three different values are found for the FWHM of the recorded signals, either 1 mm, 2 mm or 3mm. The majority of the recorded signals have a FWHM equal to 2 mm [gray zones in Fig. 7(a)], while the signals recorded where a major change in the tilting of the painting surface is observed or at the borders of the greater reliefs of the paint material [compare with Fig. 6(a)] are the only ones which shows a ± 1 deviation from this value [black or white areas in Fig. 7(a)]. It follows that the deviation from the 2 mm FWHM value found is mainly due to signal losses and distortions owed to edge effects while the variations of the FWHM related to the painting thickness and/or optical properties do not appear because minor of the adopted dZ step (1 mm).     Fig. 8(b), image of the position at which the maximum of the signals is found along the Z axis].
In fact, when the painting surface is displaced from the beam focal point, the collection of the reflected signal suffers severe losses. This is because the light cone of the reflected radiation exceeds the acceptance angle (numerical aperture) of the collecting lens, so that the exceeding part is lost. This is shown in Fig. 9(a), which shows the reflected signals recorded by the SynViewScan system at normal incidence from a metal plate for 21 different Z positions above and below the beam focal point. For the data acquisition, the metal plate has been positioned perpendicularly to the incident beam and it has been fixed on a manual graduated linear stage for changing its position along the Z-axis. The amplitude signal drop from the beam waist is in this case well described (SSE: 0.1792, R-square: 0.9928, RMSE: 0.1093) by a two term Gaussian function ( ) , with a n , b n and c n being the Gaussian coefficients. This double Gaussian fit improves by a factor of 5 the SSE and RMSE with respect to a single Gaussian fit (not shown in this work). The deviation of the experimental data from the Gaussian model are primarily related to the geometrical optics, as well as the shifting of the focal point from the Z = 0 position (beam waist shifted from Z = 0 to Z = 1.58 mm). The position of the maximum is approximately located at the focal point for most of the signals recorded at each (x, y) coordinate but, because of the unevenness of the painting surface, this parameter varies following the distribution shown in Fig. 9(b). It follows that the amplitude of the signals recorded from the painting varies with the Z positions according to the Gaussian function found for the reference metal plate. Thus, the same calculated function has been used for weighting the signals recorded from the painting, thus correcting the THz image for the amplitude variations owed to the unevenness of the painting [ Fig. 10(b)]. This helped in retrieving the details of the painting depiction that were lost before the THz image correction, e.g. the flesh and dresses of the depicted man [compare

Conclusions
Terahertz frequency modulated continue waves system (300 GHz) has been proficiently applied to the investigation of two easel paintings of different ages. Both the thickness variations of the paint layer and the nature of the artistic materials used are capable of providing a contrast in the recorded THz images, thus giving important information about the painting structure and composition. The information content of the obtained THz images has been fully inspected by selecting the appropriate THz image parameters. In addition, we have demonstrated how the Gaussian fit of the recorded signals provides more precise values for the determination of the location of the signal maximum and for the signal FWHM, thus significantly improving the level of the detail held by the THz images. The dependency of the amplitude of the reflected signal on the beam optical path has been characterized and modeled for the scanner used during the experiments. Thus, we have successfully corrected the amplitude of the recorded signals for the displacement of the object surface from the beam focal point by means of the calculated Gaussian model. The access to subtle hidden technical information about painting gives a new interest for FMCW imaging in heritage science.