Open Access
Add to BibSonomy
Abstract
Imaging using terahertz time-domain spectroscopy is a valuable diagnostic tool for material inspection. However, in the case of samples with inhomogeneous shape and composition, the reliable extraction of spatially varying dielectric properties can be very challenging. Here, we demonstrate a new approach which combines THz-TDS with photogrammetric reconstruction. We show that this technique can be used to estimate the local refractive index of samples with a complex geometry. We employ this method to study samples of ancient pottery, and demonstrate that THz techniques can provide a valuable new tool for this branch of archaeological science.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
A major and ever-recurring task of archaeologists is to characterize discovered finds, especially ceramics, which in many parts of the world are the most ubiquitous category of finds [1]. Typological classification is a long-established and undoubtedly the most widely used method, which requires comparing vessel shape and decoration of newly found sherds with well-described specimens of known provenance and chronology. Because this method depends on diagnostic features such as rims, handles and bases, many sherds, invariably the majority of finds, cannot be identified. An alternative method that was first pioneered in the 1950s is based on the ceramic fabric or paste - which is present in every typologically undiagnostic sherd [2]. Descriptions of ceramic fabrics tend to be primarily qualitative, using microscopic petrographic identification of minerals and a range of chemical analyses, including trace element ‘finger printing’ [3,4]. Quantitative methods have also been proposed, relying mostly on low-power microscopic counts and visual estimates to characterize the overall make-up of the ceramic fabric [1,5]. One important variable is the density of mineral inclusions, which is usually roughly established through comparison to lab-prepared samples with known densities. Recent attempts to develop more quantitative metrics include CT and MicroCT [6,7], polarized light microscopy, and scanning electron microscopy [8]. Another technique that has great potential for characterization of these non-metallic ceramics is terahertz imaging. Terahertz time-domain spectroscopy (THz-TDS) is sensitive to density changes in ceramics [9,10], and can in principle provide sufficient spatial resolution to characterize samples with heterogeneous composition.
THz-TDS employs radiation between 0.1 and 5 THz, and can be used to investigate various dielectric materials [11]. The radiation is not ionizing and is therefore well suited for non-destructive inspection [12,13]. Moreover, in case of multilayered samples, the short coherence length and sub-cycle temporal resolution of the detection enables characterization of the properties of each distinct layer [14,15]. THz-TDS images can be acquired both in transmission and reflection geometries.
The most common method for acquiring THz images still relies on mechanical scanning of the sample through the focus of the THz beam, and acquiring image data one pixel at a time [12]. Recent advances in focal-plane arrays have led to dramatic new possibilities for video-rate imaging; however, in contrast to THz-TDS images, these camera systems typically do not provide information about sample thickness or frequency-dependent dielectric properties [16]. Yet, approaches which rely on TDS still face important challenges. In the case of reflection imaging, one well-known problem is that surfaces must be held nearly parallel to the beam propagation axis for each pixel measurement, or else the specularly reflected signal misses the detector. In the case of a sample with non-parallel front and back surfaces, it is difficult to orient the sample such that both the front-surface and back-surface reflections are simultaneously detected with equal detection efficiency. As a result of this sensitivity to sample orientation, observed changes in the amplitude and spectral content of measured THz waveforms may not be solely due to material properties.
As a first step to address this problem, Stuebling et al. [17] developed a THz tomographic system in reflection mode using a robotic arm. This system relies on an optical image-based reconstruction of the front surface of the object under study to reorient the THz transmitter such that the beam incident on the sample remains normal to the sample’s front surface for each measured pixel. This setup has been used to analyze ancient human remains [18] and irregularly shaped wooden samples [19]. The data analysis algorithms are based on identifying reflections due to internal interfaces. Sparse deconvolution [20] is used to output a sample response function in the time domain. The time delay between the peaks of the response function relate to the thickness of each layer. This procedure is based on two key assumptions; first, that the sample’s layers are homogeneous, and therefore that they have the same refractive index throughout the entire layer, and second, that the front and back surfaces of each layer are locally approximately parallel to each other. Also, the refractive index of each individual layer is assumed to be known from separate measurements.
Here, we propose an extension to this robotic imaging method which permits us to relax some of these assumptions. This enables us to estimate the refractive index of samples with arbitrary shape, assuming that they have only a single layer (i.e., only two reflecting surfaces, front and rear). Our approach relies on image fusion, combining the aforementioned robotic-arm-based THz-TDS tomography system with photogrammetric reconstruction. We illustrate the power of this technique by applying it to the study of Roman and Punic era ceramic potsherds excavated from an archaeological site in Sardinia, Italy.
2. Imaging procedure
The process of estimation of the refractive index can be divided into three major parts: acquisition of THz data, photogrammetric reconstruction, and point cloud manipulation of imaged 3D geometries to obtain local refractive index information. Here, we introduce each of these steps and describe how each procedure is applied.
2.1 THz tomography system
We use the setup introduced in [17] to acquire THz-TDS data. This setup consists of a robotic arm (ABB IRB 120) with emitter and detector photoconductive antennas mounted in a monolithic package at the end of it. The transceiver is a THz-TDS system in a non-collinear reflection geometry. The radiation is focused on a sample and the reflected signal is collected using a set of four off-axis parabolic mirrors, which are fixed relative to the antennas. The angle of incidence on the sample is around $9.5^{\circ }$. The diameter of the spot size is 0.8 mm. Prior to acquiring THz-TDS data, the sample is first scanned using structured light. A mesh of the surface is then available for simulation of the path of the robotic arm. To begin with, the mesh is divided into small regions according to the desired imaging resolution. Next, the normal vectors of all regions are calculated. Lastly, the robot’s control system constructs the optimal path to access all regions, avoiding collisions with the sample or mounts. After the path is known, the robotic arm approaches each position on the sample always normal to the local surface at this position and at the same distance from it. THz time-domain waveforms are acquired and saved alongside the XYZ-position of each measurement point.
2.2 Photogrammetry
Next, we use photogrammetry to recover the 3D shape of the object from its photographs. Photogrammetry is the process of reconstructing a 3D representation of an object from a group of photographs taken at multiple angles. It should be emphasized that the information about exact positioning of cameras and their directions is not needed. This information is determined during the reconstruction. Moreover, information about the camera settings, e.g., focal length, is usually saved as metadata by modern cameras, and is accessible to reconstruction software. In general, the accuracy of photogrammetric reconstruction depends on the resolution of the camera chip, working distance, and the distance between different viewpoints, known as the baseline [21]. The main challenge of the reconstruction is the arbitrary scale. The object needs to be rescaled by comparing the scale used in reconstruction to a known scale. Usually, one places a ruler alongside the imaged object. The ruler and the object of interest must form a solid object, and be reconstructed together. The uncertainty in overall scale represents the most significant source of error in our refractive index data.
In this work, the object was placed on a holder that could be rotated by hand. Therefore, we changed the camera positions virtually - we moved the object relative to the camera. We acquired 50 images (6000x4000 pixels) at approximately 15 cm from the object. The baseline was also approximately 15 cm. With these values, the relative accuracy of the reconstruction is more than an order of magnitude better than the resolution of the THz imaging system. Available software packages (e.g., 3DF Zephyr, [22]) can be used for 3D reconstruction, including the construction of point clouds and meshes to represent surfaces.
2.3 Point cloud manipulation and estimation of refractive index
The next step in our analysis is the calculation of the sample’s thickness in the direction of propagation of THz radiation inside the sample. One first registers the front surface using the structured light scan mentioned above. It is important to use this as the front-surface reference, because this structured-light mesh is registered to the location of the robotic arm/THz system for each pixel. Also, this mesh does not possess the overall scale ambiguity that is associated with the photogrammetric reconstruction. We then use another software package (CloudCompare [23]) for image fusion to align the structured-light mesh to the 3D reconstruction. We also use this software to sample the front and back surfaces with a high number of points. Then, at each point on the front surface the normal vector is calculated by considering neighboring points. These points form one end of the cylinder, which propagates towards the back surface. The points, which are inside the cylinder are used to calculate the distance between both ends [24,25].
The final step of our analysis is to convert the THz measurements into average refractive index values. Each THz waveform measured in a reflection geometry consists ideally of one pulse reflected from the front surface and a second time-delayed pulse reflected from the back surface. Given that we know the propagation angle of the incident THz pulse, and we also know the sample thickness along this direction, we can use the temporal delay between these two pulses to estimate the value of the refractive index of the sample, averaged along the line connecting the front and back surfaces at the measurement location. Of course, this procedure assumes that the material dispersion is negligible along the two-way propagation path of the THz pulse [10]. Inhomogeneities in the sample composition along this path are averaged, but inhomogeneities from pixel to pixel are preserved. As a result, this procedure produces a map of the column-averaged refractive index across the sample’s surface. This map can be used to estimate both the average refractive index of the sample and its standard deviation due to inhomogeneities in material composition.
3. Results
To illustrate the use of the technique in the evaluation of archaeological samples, we have first studied three potsherds recovered from an archaeological site in the Sardara district in west central Sardinia (Italy) as part of the ‘Riu Mannu’ survey project [26]. Two of these sherds had previously been identified as Roman fine ware of the so-called Terra Sigillata type, datable to the 3rd century CE, while the third is older, of a Punic cooking pot datable to the 4th-3rd century BCE. The difference in date corresponds to the two occupation phases of the site in question (see Table 1 for details). For reference, we also image a modern (20th century) sample of Sardinian pottery. Fragments 1-3 were collected at the site of Sant’Anna in transect 20 of the Riu Mannu survey; fragment 4 belongs to a reproduction of traditional 19-20th-century Sardinian pottery, made in Assemini (SW Sardinia). Fabric A (sherd #1) is described both qualitatively and quantitatively in [27]. Figure 1 shows the result of the photogrammetric reconstruction of sherd #2, from two different viewpoints.
Table 1. Summary of information about the four samples studied.
Figure 2 shows an exemplary THz-TDS trace. We apply a bandpass filter from 0.1 to 2 THz to eliminate low-frequency oscillations and high-frequency noise. The filtered waveform contains two peaks, corresponding to the front-surface and back-surface reflections. Obviously, the back-surface reflection is reduced in amplitude considerably, due to absorption and scattering losses during the two-way travel through the sample. For our analysis, we require that this second pulse’s amplitude exceeds a certain threshold, or else the signal from that pixel is rejected and not analyzed further. The standard deviation of the signal in the range of the expected arrival of the second pulse must exceed a certain level that is determined individually for each sample by investigating a few examples of traces with and without the second pulse. For sherd #2, approximately 43 percent of the measured pixels were rejected on this basis.
Figure 3 illustrates the thickness and refractive index values estimated using the procedure outlined above, superimposed on the 3D renderings of the sample. We clearly see a variation in the thickness profile from the left to the right. In Fig. 3(b), the grey pixels indicate locations where the second pulse in the THz waveform was too small, and thus failed to fulfill the requirements necessary to extract a value for the refractive index.
Fig. 3. (a) Estimated thicknesses in mm and (b) refractive indices. Grey points show positions where the refractive index could not be calculated.
From measurements on all four of the samples mentioned above, we can compare the extracted refractive index values, using all pixels for which a value could be obtained. These results are shown in Fig. 4. Table 2 provides the median and standard deviation for all four samples, extracted from these data. Interestingly, the THz measurements indicate that the two Roman samples (sherds #2 and #3) are quite similar, which is not surprising given that they are from the same fabric (fine wear); in contrast, the results from the Punic sample (sherd #1), which is several hundred years older than the Roman samples, has a larger median, as well as a larger standard deviation, consistent with its classification as a coarse fabric. Almost 50% of all refractive index values of sherd #1 lie above 75% of the values of both sherds #2 and #3. We note that the result for the modern sample (sherd #4) has about the same standard deviation as the two Roman samples, but a somewhat lower median. The lower median value obtained from this sample (n $\approx$ 1.979) is consistent with the refractive index obtained from a traditional transmission THz-TDS measurement [11], a technique which averages over a large area of the sample (n $\approx$ 1.93). We note that this 20th century fragment is the only one of the four samples which can be measured accurately using this traditional method of transmission THz-TDS, since it is the only one with any region that has nearly flat and parallel surfaces. This emphasizes the significance of developing a technique which is suitable for the study of samples with non-parallel facets.
Table 2. Median and standard deviation of the extracted refractive index values.
4. Discussion
In Fig. 3(b), we observe a slight increase of the refractive index towards the right side. Moving even farther to the right, we see an increasing number of points where the algorithm failed because the second reflected pulse is too small. One possible explanation for this could be a tilt of the bottom surface, so that the reflected radiation partially misses the detector. Even if the reflected radiation can be collected by the detector, its path inside the sample is not equal to twice the sample thickness, due to the tilted path of the returning ray. Therefore, we should consider the impact of a tilt of the back surface on the estimated refractive index.
To analyze this situation, we derive an expression for the propagation time inside the sample as a function of the tilt angle of the back surface. The schematic of our geometry is shown in Fig. 5(a). The radiation impinges on the sample’s front surface from air ($n=1$) at an angle $\alpha$. The transmitted radiation hits the opposite surface at an angle $\beta$. The surface is inclined at an angle $\phi$. Thus, $\theta =90-\left (90-\phi \right )-\beta =\phi -\beta$ and $\Psi =\beta +2\theta =2\phi -\beta$. The length of the path inside the sample is given by the sum of the incoming and outgoing paths. Assuming a uniform refractive index $n$, the resulting propagation distance inside the sample is given by:
Fig. 5. (a) Ray tracing in the case of a tilted surface. Red and blue lines are incoming and outgoing paths, respectively. (b) Percent variation in the propagation distance inside the sample, according to Eq. (1).
Figure 5(b) illustrates the percent variation in propagation distance, as a function of the inclination angle $\phi$ of the back surface, for an assumed (constant) value of the refractive index of $n$ = 2. It is important to note that the sign of the angle of inclination cannot be obtained by calculating the angle between the normal vectors of the front and rear surfaces. Therefore, we need to consider both possible inclination directions, $\phi$ > 0 and $\phi$ < 0.
The computed variation in propagation path length can be used to calculate the influence of a tilted back surface on the estimated refractive index. If $\Delta t$ is the time delay between the two measured reflected pulses, then we can use $n=c \Delta t / l_{total}$ to compute the refractive index values corresponding to the maximum positive and negative values of the tilt angle that can be found among the pixels in Fig. 3 for which a determination of the refractive index is possible. The difference between the lowest and highest values obtained in this way was 0.58%. This defines the maximum range of uncertainty due to the effect of a tilted back surface.
To compare this predicted variation to our measured results, we select nine regions on the sample, each containing at most 12 refractive index points. For each region, the mean and standard deviation were calculated. The refractive index increases from less than 1.98 to more than 2.02 as one moves from the left to the right side of the sample. This corresponds to variations of nearly 2% across the sample. This variation is more than three times larger than what could be accounted for by a tilt of the back surface of the sample. We therefore determine that the measured inhomogeneity in the refractive index reflects a true variation in the material properties, and not merely an artifact due to back-surface tilt. We thus conclude that the measured results reported in Table 2 accurately reflect the relative material properties of the four samples studied.
5. Conclusion
We have described an advanced image fusion method which combines terahertz imaging with photogrammetry to permit the examination of samples with complex shape and morphology and inhomogeneous material composition. This method combines highly accurate photogrammetric reconstruction and tomographic THz imaging methods to estimate the refractive index of arbitrary shaped samples in each pixel of an image. Our method can be used to compare samples by comparing distributions of refractive indices. We illustrate the concept through a study of several ancient potsherds, demonstrating that the coarse wear sherd has a larger average index in the terahertz range, and also a greater degree of inhomogeneity. The larger median refractive index could result from a greater density of mineral inclusions, which have a higher index than the background clay matrix. As such, the degree of homogeneity is likely to provide a quantitative measure of the texture and density of the ceramic fabric. In both cases, the results obtained here are consistent with previous classifications of these sherds as indicated in Table 1. Further investigation will correlate these measured properties with other features of the samples under study, and will build a larger database of results in order to further improve the statistical analysis.
Funding
Deutsche Forschungsgemeinschaft (KO 1520/12-1).
Acknowledgments
We acknowledge funding from the Deutsche Forschungsgemeinschaft through project DFG KO 1520/12-1. In addition, DMM gratefully acknowledges the support of the Alexander von Humboldt Foundation.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. C. Orton and M. Hughes, Pottery in Archaeology. Cambridge Manuals in Archaeology (Cambridge University, 2013).
2. A. Shepard, Ceramics for the Archaeologist (Carnegie Institution, 1956).
3. D. Albero Santacreu, Materiality, Techniques and Society in Pottery Production: The Technological Study of Archaeological Ceramics Through Paste Analysis (De Gruyter, 2014).
4. I. Whitbread, Fabric description of archaeological ceramics, in The Oxford Handbook of Archaeological Ceramic Analysis, A. Hunt, ed. (Oxford University, 2017).
5. P. Stienstra, “Systematic macroscopic description of the texture and composition of ancient pottery - some basic methods,” Newsletter of the Department of Pottery Technology 4, 29–48 (1986).
6. A. Barron, M. Turner, L. Beeching, P. Bellwood, P. Piper, E. Grono, R. Jones, M. Oxenham, N. K. T. Kien, T. Senden, and T. Denham, “MicroCT reveals domesticated rice (Oryza sativa) within pottery sherds from early Neolithic sites (4150-3265 cal BP) in Southeast Asia,” Sci. Rep. 7(1), 7410 (2017). [CrossRef]
7. A. Machado, D. Oliveira, H. Gama Filho, R. Latini, A. Bellido, J. Assis, M. Anjos, and R. Lopes, “Archeological ceramic artifacts characterization through computed microtomography and X-ray fluorescence,” X-Ray Spectrom. 46(5), 427–434 (2017). [CrossRef]
8. K. S. Park, R. Milke, E. Rybacki, and S. Reinhold, “Application of Image Analysis for the Identification of Prehistoric Ceramic Production Technologies in the North Caucasus (Russia, Bronze/Iron Age),” Heritage 2(3), 2327–2342 (2019). [CrossRef]
9. S. Niijima, M. Shoyama, K. Murakami, and K. Kawase, “Evaluation of the sintering properties of pottery bodies using terahertz time-domain spectroscopy,” J. Asian Ceram. Soc. 6(1), 37–42 (2018). [CrossRef]
10. S. Niijima, M. Shoyama, and K. Kawase, “Nondestructive inspection of sinterability of ceramic tiles by terahertz spectroscopy,” Electron. Comm. Jpn. 102(6), 19–24 (2019). [CrossRef]
11. D. Grischkowsky, S. Keiding, M. van Exter, and C. Fattinger, “Far-infrared time-domain spectroscopy with terahertz beams of dielectrics and semiconductors,” J. Opt. Soc. Am. B 7(10), 2006–2015 (1990). [CrossRef]
12. D. M. Mittleman, R. H. Jacobsen, and M. C. Nuss, “T-ray imaging,” IEEE J. Sel. Top. Quantum Electron. 2(3), 679–692 (1996). [CrossRef]
13. I. Duling and D. Zimdars, “Revealing hidden defects,” Nat. Photonics 3(11), 630–632 (2009). [CrossRef]
14. D. M. Mittleman, S. Hunsche, L. Boivin, and M. C. Nuss, “T-ray tomography,” Opt. Lett. 22(12), 904–906 (1997). [CrossRef]
15. S. Hunsche, D. M. Mittleman, M. Koch, and M. C. Nuss, “New Dimensions in T-Ray Imaging,” IEICE Transactions on Electronics E81C(2), 269–276 (1998).
16. N. Nemoto, N. Kanda, R. Imai, K. Konishi, M. Miyoshi, S. Kurashina, T. Sasaki, N. Oda, and M. Kuwata-Gonokami, “High-sensitivity and broadband, real-time terahertz camera incorporating a micro-bolometer array With resonant cavity structure,” IEEE Trans. Terahertz Sci. Technol. 6(2), 175–182 (2016). [CrossRef]
17. E.-M. Stübling, Y. Bauckhage, E. Jelli, B. Fischer, B. Globisch, M. Schell, A. Heinrich, J. Balzer, and M. Koch, “A THz Tomography System for Arbitrarily Shaped Samples,” J. Infrared, Millimeter, Terahertz Waves 38(10), 1179–1182 (2017). [CrossRef]
18. E.-M. Stübling, A. Rehn, T. Siebrecht, Y. Bauckhage, L. Öhrström, P. Eppenberger, J. C. Balzer, F. Rühli, and M. Koch, “Application of a robotic THz imaging system for sub-surface analysis of ancient human remains,” Sci. Rep. 9(1), 3390 (2019). [CrossRef]
19. K. Krügener, E.-M. Stübling, R. Jachim, B. Kietz, M. Koch, and W. Viöl, “THz tomography for detecting damages on wood caused by insects,” Appl. Opt. 58(22), 6063–6066 (2019). [CrossRef]
20. J. Dong, X. Wu, A. Locquet, and D. S. Citrin, “Terahertz Superresolution Stratigraphic Characterization of Multilayered Structures Using Sparse Deconvolution,” IEEE Trans. Terahertz Sci. Technol. 7(3), 260–267 (2017). [CrossRef]
21. T. Luhmann, S. Robson, S. Kyle, and J. Boehm, Close-Range Photogrammetry and 3D Imaging (De Gruyter, 2014).
22. 3Dflow srl, “3DF Zephyr Free Edition, .,”.
23. CloudCompare, “CloudCompare, GPL Software, .,”.
24. P. J. Besl and N. D. McKay, “A method for registration of 3-D shapes,” IEEE Trans. Pattern Anal. Mach. Intell. 14(2), 239–256 (1992). [CrossRef]
25. Y. Chen and G. Medioni, “Object modelling by registration of multiple range images,” Image Vis. Comput. 10(3), 145–155 (1992). Range Image Understanding. [CrossRef]
26. M. B. Annis, P. van Dommelen, and P. van de Velde, “Sardegna. progetto riu mannu. per unarcheologia del paesaggio,” Bolletino di Archeologia pp. 765–770 (2008).
27. P. Van Dommelen and M. Trapichler, “Fabrics of western central Sardinia,” in Fabrics of the Central Mediterranean: Provenance Studies on Pottery in the Southern Central Mediterranean from the 6th to the 2nd, V. Gassner, ed. (University of Vienna, 2011).
