Terahertz frequency-wavelet domain deconvolution for stratigraphic and subsurface investigation of art painting

Terahertz frequency-wavelet deconvolution is utilized specifically for the stratigraphic and subsurface investigation of art paintings with terahertz reflective imaging. In order to resolve the optically thin paint layers, a deconvolution technique is enhanced by the combination of frequency-domain filtering and stationary wavelet shrinkage, and applied to investigate a mid-20th century Italian oil painting on paperboard, After Fishing, by Ausonio Tanda. Based on the deconvolved terahertz data, the stratigraphy of the painting including the paint layers is reconstructed and subsurface features are clearly revealed, demonstrating that terahertz frequency-wavelet deconvolution can be an effective tool to characterize stratified systems with optically thin layers.


Introduction
Terahertz (THz) imaging, which can provide a nondestructive, noncontact, and nonionizing modality for characterizing materials, is emerging as a promising technique in the field of cultural-heritage conservation [1].The THz portion of the electromagnetic spectrum extends from ∼100 GHz to 10 THz, lying between the microwaves and infrared; and the wavelength range in this region is 3 mm down to 30 µm.THz waves can penetrate numerous nonmetallic materials which may be opaque in the range of visible and infrared light.Therefore, THz imaging can provide information in depth to investigate the stratigraphy and hidden features or damages in artworks, which is of great interest for conservation and restoration, as well as for art-historical studies.In recent years, THz imaging has shown great potential in the cultural-heritage arena, enabling the characterization of features in stone sculpture and in architecture [2], wooden objects [3] and tree rings [4], clay artifacts [5], written papyrus [6], and even ancient Egyptian mummies [7,8].Within the area of art paintings, THz imaging has been utilized to identify features hidden under the surface of a painting, such as the artist's signature [9], under-drawings and modifications due to the reuse of an earlier painting or canvas [10][11][12][13].In addition to revealing such subsurface features, the stratigraphic structure of art paintings is of general interest, since it may be used to understand the process of the painting's creation.
Our aim is to apply THz reflective imaging to reveal a painting's stratigraphy and hidden features.The approach relies on reconstruction of the buildup of layers of the painting based on monitoring the THz echoes reflected from the various interfaces [14].By layer, we mean apparent layers of paint, the apparent interfaces between which give rise to a reflected pulse that we register.Clearly, a painter does not apply uniform layers of paint, but uses multiple overlapping brushstrokes (or in our case swipes of a palette knife) using a continuous range of pigments that may be mixed one into the other.And if several pigments present similar refractive indices, we might never be able to distinguish consecutive layers containing such pigments.Moreover, if the application of paint is sufficiently thin, the reflected THz signal might not contain features that are sufficiently distinct to detect that "layer".Be that as it may, when the layers are optically thick (i.e., compared to the wavelengths within the THz spectrum) and the refractive-index contrast between layers is sufficient, the reconstructed profile in depth based on the raw THz waveforms is satisfactory.Usually, the interface between the paint and the substrate of the painting can be clearly identified.However, the assumptions above frequently do not hold when we only focus on the paint region and try to resolve the stratigraphic structure of paint layers, since the thickness of each paint layer may be tens of microns (or even smaller) and the difference in refractive index between different pigments or colors may be quite small [15].Therefore, the echoes will partially or totally overlap.Moreover, the signal-to-noise ratio may be quite low, leading to unsuccessful reconstruction in depth.THz frequency-wavelet deconvolution technique is utilized in this study to solve this problem, and is shown to be quite effective for stratigraphic and subsurface investigation of optically thin paint layers.In short, THz deconvolution retrieves the impulse response function from the inverse Fourier transform of the transfer function, which should ideally consist of a sequence of ideal impulses corresponding to reflections from the interfaces of a stratified system [16].The approach of frequency-wavelet deconvolution is to enhance the deconvolution process by first employing frequency-domain filtering and then further improving the signal-to-noise ratio by wavelet denoising techniques [17].The procedure is illustrated in detail in this study.In order to demonstrate this algorithm, a mid-20th century Italian painting, After Fishing by Ausonio Tanda [shown in Fig. 1(a)] is systematically investigated.This oil painting's specific pigments are unknown, and the substrate is unbleached paperboard.By applying THz frequency-wavelet deconvolution, the stratigraphy of the painting including the paint layers is successfully reconstructed and subsurface features are clearly revealed.

Theory
In the time domain, the THz reflected signal (electric field) r(t) is the convolution of the incident THz pulse i(t) with the impulse response function h(t), which corresponds to the structure and properties of the painting at a given two-dimensional position, Deconvolution retrieves the impulse response function h(t) by applying the inverse Fourier transform based on the convolution theorem, where FFT denotes the Fourier transform and FFT −1 the inverse Fourier transform.Frequently, successful deconvolution cannot be expected by directly applying Eq. ( 2), since division of small numbers will give rise to large abnormal values as high spikes, especially in the high frequency region, leading to severe ringing in the time domain.Therefore, the deconvolution process is usually further modified by applying frequency-domain filtering function to suppress the high frequency noise, which can be expressed as, with f (t) as the filter function in the time domain.By the convolution theorem, Eq. ( 3) can be transformed to h (t) can be considered as a reconstructed signal obtained by inputing a new function f (t) instead of i(t).In order to obtain a successful reconstruction, the time duration of the new input function f (t) should be short enough to resolve the time intervals between featured echoes, and the new input function should not contain extra signal cycles before or after the main peak, which will obscure the real featured echoes in the reconstructed signal; however, the selection of f (t) is also a compromise between time resolution and frequency-domain filtering [18].If the duration of f (t) is too short, its frequency spectrum will include large spikes at high frequencies, which will degrade the reconstructed signal in the time domain.The double Gaussian filter and Wiener filtering can be selected to serve as the frequency domain filtering [17], and a tapered cosine apodisation function has also been found to work well [19].Considering the complexity of implementation and effectiveness, the Hanning window function or von Hann filter is chosen as the new input f (t) in this study, and its frequency spectrum F(ω) can be expressed as where t o corresponds to the arrival time of the main peak in the time domain, f c is the cutoff frequency.The effectiveness of this frequency domain filtering is easy to manipulate just by changing the cutoff frequency f c .An example of f (t) and its frequency spectrum F(ω) for typical parameters (t 0 = 10 ps and f c = 4 THz) is shown in Fig. 2. Quite often, deconvolution only with frequency-domain filtering cannot guarantee a satisfactory signal-to-noise ratio when a relatively high value of f c is selected.Stationary wavelet shrinkage is applied to further attenuate the residual noise.This technique decomposes a 1D signal into the approximation coefficients vector and detail coefficients by convolving with a low-pass filter and a high-pass filter along the temporal axis at each level [20].Wavelet coefficients with small absolute values can be considered as noise, and wavelet coefficients with large absolute values are regarded as the main featured information of the signal [17].Removing the small absolute-value coefficients by thresholding and then reconstructing the signal is expected to produce a signal in which the contribution of noise has been reduced.Detailed information about the stationary wavelet shrinkage can be found in [20].A typical THz waveform reflected from the painting in our study is shown in Fig. 3(a).Echoes reflected from the top and bottom surfaces of the painting can be clearly identified, and overlapping echoes due to the internal structure just following the first peak can also be observed.THz frequency-wavelet deconvolution described above, is applied to this THz waveform, and a comparison is made with two different cutoff frequencies.For f c =2 THz, shown in Fig. 3(b), the signal after deconvolution with frequency-domain filtering exhibits a good signal-to-noise ratio; however, the peak indicating the underlying paint layer is not well separated; for f c =3 THz, shown in Fig. 3(c), the hidden peak can be well separated and clearly identified, although the signal-to-noise ratio is relatively low.The residual noise in both cases can be effectively attenuated with the wavelet de-noising.In our study, we choose the symlet (sym4) wavelets, which are a modified form of Daubechies wavelets.A maximum level of 7 is used for the wavelet decomposition as no significant improvement can be observed for higher levels to justify the extra computational expense.The positions of the noise intervals are set to the region of the fluctuations at the beginning and the end of the signal, which are boxed in Fig. 3. Sometimes, the signal after frequency-wavelet deconvolution contains slow fluctuations corresponding to the low frequency noise due to the THz source being deficient in the low THz frequency region.This kind of low-frequency noise can be canceled by subtracting the baseline of the deconvolved signal.The final deconvolution result of this typical THz waveform is shown in Fig. 3(d), which clearly exhibits the stratigraphy of the painting including one underlying paint layer at this single pixel.A THz time-domain spectroscopy (TDS) system (Teraview TPS Spectra 3000) [21] was employed to scan After Fishing.The GaAs photoconductive antenna is excited by an ultrafast laser to produce roughly single-cycle THz pulses with bandwidth extending from 60 GHz to 3 THz.The focus spot size of the THz beam is frequency-dependent, and is about 300 µm at 1 THz.THz imaging was performed in reflection at almost normal incidence from the top surface of the painting.Before scanning the painting, a THz reference signal was recorded by setting a metal plate at the sample position, which is shown in Fig. 4. The painting was raster-scanned by a set of motorized stages moving in the X and Y directions with a 1 mm spatial step over a 28 cm × 25 cm region of the painting, corresponding to 280 × 250 pixels.Each recorded reflected THz waveform contains 4096 data points, and the signal is averaged over 10 shots per pixel.Due to the curvature of the painting, two foldback clips, which can be seen in the following images, were carefully used at the edges to keep the surface relatively flat.The scanning of the painting was conducted in a temperature-controlled laboratory at 22 • C. The humidity in the laboratory was about 48 %.After performing the scanning, the positions of the maximum peak of all the waveforms, corresponding to the interface between the air and the paint surface, were aligned to the same position temporarily for getting rid of the surface curvature and the convenience of further signal processing.After this alignment, a 3D volume raw data set was acquired.

Results and discussion
First, we show the THz images based on the raw THz data directly obtained from the scan.The THz C-scan in the time domain is shown in Fig. 5, in which the contrast mechanism chosen is the maximum peak amplitude in the reflected signal.The THz C-scan shown here can provide information about the surface features of the painting, including the reflectivity of the surface paints and the surface roughness.THz C-scans in the frequency domain can also be obtained by taking the Fourier transform of the raw waveform at each pixel, shown in Fig. 6.The contrast mechanism chosen is based on the spectral power density at a given frequency (0.4, 0.6, 0.8, and 1.2 THz), then the contrast is normalized to one in order to make a clear comparison.Higher frequency, corresponding to shorter wavelength, is helpful for observing small and subtle features.For example, as we progress from Fig. 6(a) to Fig. 6(d), in the upper left region (between pixels 0 and 50 in the X direction, and 150 and 200 in the Y direction) a ripple pattern with periodicity approximately in the X direction emerges as the cutoff of the frequency window increases.Note that 1.2 THz corresponds to a wavelength of ∼ 0.25 mm, which is on the length scale of the periodicity of the ripples.These features are due to the clarity of the features in the THz C-scan at 1.2 THz, we conclude that the signal-to-noise ratio at 1.2 THz is still acceptable; consequently, we then can set the cutoff frequency f c higher than 1.2 THz to obtain higher time resolution in the following deconvolution procedure.By contrast with the X-ray transmission image, shown in Fig. 1(b), one striking feature we observe are several small and blurred spots randomly distributed in the frequency-domain THz C-scans; these spots are not observed in the X-ray transmission image.By carefully checking the waveforms at relevant pixels, small echoes can be found in the time interval 12 to 28 ps, corresponding to the paperboard support.The THz C-scan is plotted in Fig. 7(a) based on the peak signal value between 12 and 28 ps providing the contrast.THz B-scans can also be plotted in order to locate the spots in depth.Figure 7(b) shows the THz B-scan with the cross-section Y =176, which provides the locations in depth of four typical spots [highlighted in Fig. 7(a)].These spots are a result of defects in the paperboard, which we speculate may be due to oil and/or biological growth.Organic materials (which contain mainly light atoms like C, H, O, N) exhibit low radiological contrast [22]; therefore, these spots in the paperboard are difficult to detect by standard X-ray systems.THz B-scans based on the raw data cannot provide the stratigraphic detail, since separation of the paint layers above the paperboard cannot be resolved.Artifacts (weak horizontal lines) due to ambient water-vapor [20] and the system noise also lower the quality of the THz B-scans.
Next, we apply the THz frequency-wavelet deconvolution algorithm to each waveform in the raw THz data.The cutoff frequency f c chosen for all the waveforms is 3 THz.In order to reconstruct the stratigraphy and investigate subsurface features, peak detection is performed on the deconvolved signals, since the peak corresponds to the existence of an interface.A threshold value for peak detection is set for all pixels, above which we consider a feature as a valid peak.Based on the results of peak detection, the deconvolved signals can be mainly classified into four types.A typical waveform of each type is shown in Fig. 8.
For Type I, shown in Fig. 8(a), a small peak occurs before the arrival of the main peak corresponding to the air/paint interface.This feature is due to the diffuse reflection on a rough surface.Since the focus spot of the THz beam has limited size, if the surface within this spot  is rough and contains height variation, some portion of the THz beams will reflected earlier than the other portion of the beam, which leads to peaks in the reflected THz signal before the dominant peak.THz C-scan based on the amplitude of the peak before the main peak for Type I signals is plotted in Fig. 9(a).This C-scan mainly displays the outlines of the drawing objects in the foreground, the paint for which may contain angles and edges.THz C-scan based on the amplitude of the main peak corresponding to the air/paint interface is also plotted in Fig. 9(b).
By comparing these two C-scans, we see that they are complementary-positive and negative images.Areas with high values in Fig. 9(a) exhibit relatively low values in Fig. 9(b), which also supports our assertion concerning the role of diffuse reflection.For Type II and Type III, one featured peak is detected between the air/paint and paint/paperboard interfaces.This peak indicates the existence of an additional underlying paint layer; therefore, there are in total two paint layers on the paperboard at this pixel.The observation of this peak is important for a successful stratigraphic reconstruction, since it is hidden in the overlapping echoes in the THz B-scans based on the raw data.For peaks with positive sign, shown in Fig. 8(b), we conclude that the refractive index of the first paint layer is smaller than that of the second paint layer.For peaks with negative sign, shown in Fig. 8(c), the refractive index of the first paint layer is larger than that of the second paint layer.THz C-scans are plotted based on the amplitude of this featured peak, shown in Fig. 10. Figure .10(a) exhibits the regions with underlying layers, of which the refractive index is larger than that of the surface pigments; on the other hand, Fig. 10(b) displays the regions with underlying layers, of which the refractive index is smaller than that of the surface pigments.
THz B-scans can also be reconstructed based on the deconvolved signals to reveal the stratigraphy of the painting.Based on the peak-detection (both positive and negative peaks) procedure, binary THz B-scans (we consider the region with a valid peak as '1' and the other areas as '0') can be obtained.A more reasonable way to reconstruct the THz B-scan is to regard the bottom surface of the paperboard as flat; therefore, the deconvolved signals for binary THz For Type IV, shown in Fig. 8(d), no extra featured peak can be identified.Therefore, we conclude that there is only one paint layer on the paperboard or the thickness of the underlying layer is too thin to be detected in the THz regime.
The deconvolved data can also provide quantitative information on the painting, such as the layer thicknesses.The optical thickness of the paint layers can be acquired by calculating the optical delay between the air/paint and paint/paperboard interfaces.The distribution of the paint-layer thickness can reveal in which areas the artist applied more layers of paint.(In our case, the thick paint application with a palette knife resulted in considerable impasto that enables us to check the reasonableness of our THz results visually or tactilely.)If the refractive indice of the pigments of different colors, and of the oil binder, are known, the physical thickness of the  paint layers can be calculated.In practice, however, it is difficult to obtain the refractive index of the pigments in an actual art painting nondestructively.In this study, we assume the refractive index to be constant in different paint layers, which enables us to roughly estimate the physical thickness distribution of the applied paint layers on the basis of the optical thickness in the THz regime.The normalized THz image of the thickness distribution based on the optical thickness of the applied paint layers is shown in Fig. 12.A visible raking light image of After Fishing is shown in Fig. 13 to support our assumption.The raking light image provides information on the surface topography, which also highlights the physical thickness distribution of the applied paint layers on the paperboard.Areas with thicker applied paint layers, as determined in Fig. 13, also show high contrast in Fig. 12. Relatively small differences in the refractive indices of various paint layers at their mutual interfaces produce the Fresnel coefficients resulting in the echoes we observe.Still, as we argue, we are able to obtain a reasonable estimate of layer thicknesses and overall thickness of the paint application.Finally, we point out the high degree of correlation between the THz image of thickness distribution in Fig. 12 and the X-ray transmission image in Fig. 1(b).Whereas Fig. 12 provides a measure of the optical thickness of the painting and an estimate of the amount of paint applied on the paperboard, Fig. 1(b) provides a measure of the total X-ray scattering strength in a given pixel, which is also a measure of the amount of material present, but also must account for the atomic number of the constituents, core levels, etc.A fused image of Fig. 12 and Fig. 1(b) with normalized intensities is created using color to distinguish the similar intensity between the two images, which is shown in Fig. 14.By using red for Fig. 1(b) and green for Fig. 12, areas of similar intensity between the two images are shown yellow in the fused image.In Fig. 14, we can observe that the yellow areas are dominant and are mostly located in the regions with relatively thicker layers.Green areas shows the regions where Fig. 12 displays a higher response, corresponding to the regions with relatively thinner layers.Red areas, where the X-ray image shows higher response due to its shorter wavelength, highlight the regions containing angles and edges on the surface of painting, which are composed of more paint.On the contrary, regions in red are considered as rough in the THz regime, and due to the diffuse reflection, the received THz signals there belong to the Type I signals, which cannot provide accurate thickness information.Although the X-ray image reveals the features which are lost in the THz image of thickness distribution, it is also noticed that red areas in Fig. 14 present high similarity to the previous THz image in Fig. 9(a).Therefore, we can conclude that THz imaging can be utilized to obtain information that is similar to that contained in conventional X-ray imaging.

Conclusion
In this study, THz frequency-wavelet deconvolution is implemented specifically for the stratigraphic and subsurface investigation of art paintings.A mid-20th Italian oil painting on paperboard, After Fishing, by Ausonio Tanda is systematically examined with the method.Utilizing THz reflective imaging, the surface features and a distribution of spots in the paperboard are clearly revealed based on THz raw data.THz frequency-wavelet deconvolution is applied to the THz raw data, and the stratigraphy of the paint layers is successfully reconstructed across the painting based on the deconvolved data.THz C-scans and B-scans are analyzed based on different types of deconvolved signals to investigate the subsurface features of the painting, including the identification of regions with more than one paint layer, the refractive-index difference between paint layers, and the distribution of the paint-layer thickness.Based on our study, we conclude that THz frequency-wavelet deconvolution can be an effective tool for the stratigraphic and subsurface investigation of art paintings, and may also be applied for the characterization of other stratified systems.
In addition, one striking point comes out of this study is the high similarity between the THz and X-ray images.THz imaging of the thickness distribution of the paint exhibits a high degree of correlation with the X-ray transmission image.Moreover, THz imaging also reveals the spots in the paperboard which cannot be identified in the X-ray image.Compared with X-ray imaging, THz imaging can be carried out in reflection, so that there is no need access to both sides of a painting; THz imaging can provide information in depth and THz data contains spectral information which can potentially be used to identify different pigments.Therefore, our results open up the way for the use of non-ionizing THz imaging as a potential substitute for ionizing X-ray analysis in nondestructive evaluation of art paintings.

Fig. 1 .
Fig. 1.Images of After Fishing by Ausonio Tanda.(a) Visible photograph of After Fishing, and (b) X-ray transmission image of After Fishing.After Fishing is not copyrighted.Marcello Melis is the owner of After Fishing.

Fig. 2 .
Fig. 2. Hanning window function with typical values, t 0 =10 ps and f c =4 THz, in the time domain and its Fourier transform (power spectrum) in the inset.

Fig. 3 .
Fig. 3. THz frequency-wavelet deconvolution result for a typical THz waveform.(a) Original THz reflected waveform from one typical pixel; (b) deconvolution result with f c =2 THz; (c) deconvolution result with f c = 3 THz; (d) final deconvolution result to reveal the stratigraphy in this typical pixel.The insets in (b) and (c) show the zoom-in of the signals around the time scales of interest.

Fig. 4 .
Fig. 4. THz reference signal with its frequency spectrum in the inset.

Fig. 5 .
Fig. 5. THz C-scan image in time domain based on the raw data.The contrast mechanism chosen is the maximum peak amplitude in the reflected signal.

Fig. 6 .Fig. 7 .
Fig. 6.THz C-scan images based on the raw data in frequency domain at (a) 0.4, (b) 0.6, (c) 0.8, and (d) 1.2 THz.The contrast mechanism chosen is based on the spectral power density at the corresponding frequency, then the contrast is normalized to one in order to make a clear comparison.A few small and blurred spots are highlighted with red dotted circles.

Fig. 8 .
Fig. 8. Four different types of the deconvolved signals.(a) Type I: signal with a small peak occurring before the arrival of the main peak due to diffuse reflection; (b) Type II: signal with one featured peak indicating the existence of an addition layer with refractive index larger than that of the surface pigments; (c) Type III: signal with one featured peak indicating the existence of an additional layer with refractive index smaller than that of the surface pigments; (d) Type IV: signals with no additional peak identified.

Fig. 9 .Fig. 10 .
Fig. 9. THz C-scans based on the Type I signal and the main peak indicating the surface feature of the painting.(a) The THz C-scan based on the amplitude of the peak arriving before the main peak; (b) The THz C-scan based on the amplitude of the main peak corresponding to the air/paint interface.

Fig. 11 .Fig. 12 .
Fig. 11.Comparison between THz B-scans based on the raw data and binary THz B-scans based on the deconvolved data.(a) The THz B-scan based on the raw data with 'Cross Section 1' in Fig. 10(b); (b) The binary THz B-scan based on the deconvolved data with 'Cross Section 1' in Fig. 10(b); (c) The THz B-scan based on the raw data with 'Cross Section 2' in Fig. 10(a); (d) The binary THz B-scan based on the deconvolved data with 'Cross Section 2' in Fig. 10(a).

Fig. 14 .
Fig. 14.The fused image of the X-ray image in red and THz thickness distribution image in green to distinguish the areas of similar intensity in yellow.