Phase-contrast THz-CT for non-destructive testing
Peter Fosodeder, Simon Hubmer, Alexander Ploier, Ronny Ramlau, Sandrine van Frank, Christian Rankl
PPhase-contrast THz-CT for non-destructivetesting P ETER F OSODEDER , S IMON H UBMER , A LEXANDER P LOIER , R ONNY R AMLAU , S ANDRINE VAN F RANK AND C HRISTIAN R ANKL Research Center for Non Destructive Testing GmbH (RECENDT), Altenbergerstraße 69, A-4040 Linz,Austria Johann Radon Institute Linz, Altenbergerstraße 69, A-4040 Linz, Austria Doctoral Program Computational Mathematics, Johannes Kepler University Linz, Altenbergerstraße 69,A-4040 Linz, Austria Institute of Industrial Mathematics, Johannes Kepler University Linz, Altenbergerstraße 69, A-4040 Linz * [email protected] Abstract:
Computed tomography (CT) is a well-established method for non-destructive testing(NDT), typically performed with X-radiation. Here, an alternative approach for CT based on THzradiation is presented. Optimized for the potential use-case of NDT on extruded plastic profiles,an imaging model and different resultant reconstruction approaches are derived. For experimentalvalidation, 3D printed plastic profiles were tested with a THz time-domain spectrometer intransmission geometry. On the basis of differently complex samples, it is shown that the derivedreconstruction approaches are highly robust. Quantitative analysis has further shown that theresults are highly accurate for samples without curved structures. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Over the past years, THz technology has evolved into a more and more recognized technology forindustrial non-destructive testing (NDT). Some of the advantages over other complementary NDTmethods involve the fact that THz systems are typically operated contact-less and are completelyharmless in any human environment without the need for costly shielding. One major reasonfor the increasing value of THz NDT methods is the recent development of fast and compactmeasurement systems [1]. In general, THz radiation sources are subdivided into pulsed andcontinuous wave (cw) types. Applications such as layer thickness control [2] based on pulsedsources and security control [3] based on cw sources are often mentioned as prime examples inthis context. With this work, we present and evaluate the performance of a further NDT method,namely THz computed tomography (THz-CT). In the past, THz-CT was performed using bothcw [4] and pulsed [5] sources. However, due to their shorter wavelength and therefore betterresolution, we focus on the use of a pulsed THz time-domain spectrometer (THz-TDS) [6].One potential applications for this technology is for instance in-line NDT of extruded plasticprofiles, effectively reducing the amount of produced plastic waste in case of production failure.In contrast to the well established imaging technique X-ray CT (X-CT), a THz-TDS setupenables one to not only detect the amplitude but also the phase information. In addition toreconstructing the sample geometry from its absorption coefficient, this opens up the new pathof extracting data based on the refractive index of the sample. For this purpose, we present aTHz-CT imaging model specifically suited to the application case of quality control of plasticprofiles.Despite the fact that THz imaging is strongly influenced by electromagnetic wave propagationeffects (scattering, refraction), our imaging model is based on a variant of the Radon model [7],which is typically applied in conventional X-CT. This enables us to adapt existing concepts of a r X i v : . [ phy s i c s . a pp - ph ] F e b -CT to the THz-CT problem. Note that effort has been taken previously in order to overcomethe limitations of Radon model based image reconstruction by using Helmholtz’s equation as theimaging model [8]. Other approaches deal with refraction by means of ray-tracing [9]. However,the increased mathematical complexity of the problem makes it necessary to introduce preliminaryknowledge about the sample and/or to perform computationally expensive simulations. Thesecan largely be avoided with our approach, while still achieving quantitatively accurate results.Based on our derived imaging model, we present two methods for calculating a sinogram fromthe measured THz signals. Analogously to X-CT, the first method utilizes amplitude informationonly while the second method takes advantage of the fact that in THz-TDS the signal is sampledwith sub-ps resolution, which allows us to extract a pulse phase sinogram. In the next step, weapply different image reconstruction approaches to the sinogram, namely filtered backprojection,Tikhonov regularization and Landweber iteration [7, 10, 11]. For performance testing, we workedwith 3D-printed plastic profiles and compared the outcome of the above procedure qualitativelyfor the same sample. In addition to that, two more 3D printed samples were measured andevaluated quantitatively by means of measuring their wall thickness from the reconstructed imageand comparing it to a mechanical reference measurement.
2. Materials and Methods
For our measurements, the commercially available THz-TDS system TeraFlash smart [1](TOPTICA, Munich) was used. This system uses two photoconductive antennas (PCA) and twosynchronized femtosecond lasers, operating at a center wavelength of 1560 nm for generatingand detecting pulsed THz radiation.As illustrated in Fig. 1, the THz pulses are guided and focused by 4 off-axis parabolic mirrors(OPM), and the sample is mounted in the focal plane. Sample positioning for the CT measurementis performed by one translation and one rotation stage. For imaging in transmission it is crucialto find a reasonable trade-off between depth of field and focal spot-size, therefore the effectivefocal length of the focusing OPM was dimensioned with 101.6 mm. A Gaussian focal spot witha FWHM of approximately 1 mm (measured by the knife-edge technique) was achieved whilestill preserving enough depth of field for sample diameters of several centimeters.Fig. 1 also illustrates two THz signals, one reference measurement in air and one measurementon a 3D printed plate of Polypropylene (PP) with a thickness of 2 mm. Typically, the bandwidthof a single shot THz reference signal exceeds 3.5 THz at a measurement rate of 1600 signals persecond. The total length of one THz signal is 152.3 ps, sampled with a resolution of ≈
49 fs.
The samples for this study were 3D printed from natural PP at a temperature of 240 ° C. Due tominor visible structural inhomogenities from the 3D printers layer resolution of 150 µ m, a THzspectroscopy experiment on a printed plate of PP was performed. The results in Fig. 2 werecalculated, using a frequency-domain model-based parameter estimation [12, 13]. Despite thevisible surface structure of the printed samples, no significant absorption features compared toother THz spectroscopy experiments performed on PP were found [14]. In conventional X-CT measurements, detector arrays are used in order to minimize acquisitiontime. Since our THz-TDS setup employs two PCAs, effectively acting as a static point-likeemitter/detector pair, we perform continuous line scans by moving the sample through the THzfocus, at typical velocities of 150 mm/s. After each linescan, the sample is rotated by 1 ° until thefull 360 ° scan is completed. Typical acquisition times for a full 360 ° scan are below 5 minutes x Rx O b j ec t T r an s l a t i on s xyz R o t a t i on θ T x S i gna l R x S i gna l OPM OPMOPMOPM
Fig. 1. Schematic drawing of the measurement setup. The THz radiation is generatedby a PCA (Tx). Two Off-Axis Parabolic Mirrors (OPM) are used to form a focused THzbeam. After interacting with an object, the THz beam is guided to the detector PCA(Rx) by two OPMs again. Sample movement is performed by a motorized translationand rotation stage. Exemplarily, a measured reference signal through air and a measuredsignal through a 2 mm plate of PP is shown.
Absorption Coefficient
Refractive Index
Fig. 2. Material parameters of a 3D printed plate of PP (thickness = 2 mm). Theabsorption coefficient 𝛼 and the refractive index 𝑛 were calculated from a THz-TDSmeasurement in transmission. No significant absorption features caused by potentialstructural inhomogenities were found. for the samples presented in this work, while achieving an angular resolution of 1 ° and a spatialresolution smaller than 0.1 mm. In order to improve the signal-to-noise ratio (SNR) of the acquired THz signals, a digital FIRlow-pass filter [15] with a passband frequency of 2.5 THz and a stopband frequency of 3.5 THz(attenuation: -80 dB) was applied to the raw data.Subsequently, consecutively acquired THz signals were filtered by a moving average filter witha window size of five signals. This effectively corresponds to a spatial low-pass filtering of thedataset and therefore further reduces the SNR.Due to the continuous motion of the sample, the acquired measurement positions are availableon scattered positions only. After calculating the sinogram values on the corresponding scatteredositions (see section 3.1), the values were interpolated on a Cartesian grid with a spatialresolution of 0.1 mm and an angular resolution of 1 ° .After image reconstruction, the retrieved images were thresholded in order to remove valuessmaller than zero. These values typically appear due to numerical reasons and are nonphysical.
3. The THz-Imaging Model
In its most general form, the propagation of electromagnetic radiation through an object isdescribed by the macroscopic Maxwell equations [16]. In this framework, the common phenomenaof wave propagation, such as refraction and diffraction, are considered. Recently, efforts havebeen taken in order to perform THz CT image reconstruction within this framework [8].Since our work is specifically dedicated to 3D imaging of plastic profiles, some assumptionswill be introduced in the following, leading to a more compact imaging model for our specificproblem. The basis of our imaging model is the well-known regime of geometrical optics, wherea light beam is treated as a mathematical ray, undergoing reflection, refraction and absorption bythe sample. In this regime, the material-light interaction is commonly assumed to be linear andisotropic. Furthermore, the appearance of multiple reflections within a sample is neglected, sincethey result only in minor contributions to the measured signal.Since the samples in this work represent a special case of sparsely populated objects, weintroduce a significant simplification to the imaging model by neglecting diffraction and refractionin order to preserve a reasonably complicated framework for computationally efficient imagereconstruction.Using the assumptions introduced above, the remaining radiation-matter interaction can befully described by the complex and frequency dependent refractive index of the material (cid:101) 𝑛 ( 𝜔 ) = 𝑛 + 𝑖𝜅 , (1)where 𝑛 denotes the real refractive index and 𝜅 denotes the absorption index. With the notation inEq. (1), the measured THz signal, after travelling a distance 𝑑 through the object, can be writtenas 𝐸 ( 𝑑, 𝑡 ) = +∞ ∫ −∞ ˆ 𝐸 ref ( 𝜔 ) exp (cid:104) 𝑖 (cid:16) 𝜔𝑡 − 𝜔𝑐 ( (cid:101) 𝑛 ( 𝜔 ) − 𝑛 ) 𝑑 (cid:17)(cid:105) d 𝜔 . (2)Here 𝑐 refers to the speed of light in a vacuum, 𝑛 denotes the real refractive index of air and ˆ 𝐸 ref is the Fourier transform of the reference signal 𝐸 ref measured in air.The complex refractive index (cid:101) 𝑛 is commonly approximated by its Taylor series (cid:101) 𝑛 ( 𝜔 ) = (cid:101) 𝑛 ( 𝜔 ) + dd 𝜔 (cid:101) 𝑛 ( 𝜔 ) | 𝜔 ( 𝜔 − 𝜔 ) + O (cid:16) 𝜔 (cid:17) . (3)A comparison with Fig. 2 shows that 𝑛 is only marginally varying over the frequency bandwidthof our measurement system and is therefore sufficiently approximated by a constant 0-th orderterm 𝑛 = 𝑛 ( 𝜔 ) .We will further rewrite the absorption index 𝜅 ( 𝜔 ) using the absorption coefficient 𝛼 ( 𝜔 ) = 𝜅 ( 𝜔 ) 𝜔𝑐 . Fig. 2 shows that the absorption coefficients frequency dependence isnegligibly low and we therefore introduce the constant absorption coefficient 𝛼 = 𝛼 ( 𝜔 ) as afurther simplification to our model.Inserting the average quantities 𝑛 and 𝛼 into Eq. (2) leads to 𝐸 ( 𝑑, 𝑡 ) = +∞ ∫ −∞ ˆ 𝐸 ref ( 𝜔 ) exp (cid:104) − 𝑖 𝜔𝑐 ( 𝑛 − 𝑛 ) 𝑑 (cid:105) exp (cid:18) − 𝛼 𝑑 (cid:19) exp ( 𝑖𝜔𝑡 ) d 𝜔 . (4)he approximations introduced with 𝑛 and 𝛼 allow us to transform our model back to thetime-domain analytically. A more general approach would be to proceed with the frequencydomain model in Eq. (4) in order to account for the frequency dependent material parameters.We leave this as subject to future work and continue by performing the inverse Fourier transformin equation 4. This leads to the time-domain representation of the model 𝐸 ( 𝑑, 𝑡 ) = 𝐸 ref (cid:18) 𝑡 − 𝑛 − 𝑛 𝑐 𝑑 (cid:19) exp (cid:18) − 𝛼 𝑑 (cid:19) . (5)Eq. (5) describes the interaction of a mathematical THz ray with a material of thickness 𝑑 andmaterial parameters 𝑛 and 𝛼 . One can see that the interaction between sample and THz pulseis only expressed by a time-delay of the THz pulse with respect to the reference pulse, and anabsorption term according to Lambert-Beers law.In order to describe the CT measurement procedure (see section 2.3), we introduce the discretesample positions 𝑠 𝑖 and projection angles 𝜃 𝑗 for each point measurement, and the density function 𝑓 ( 𝑥 ) , representing the sample geometry. Analogously to our previous work [17], we nowintroduce the Radon transform 𝑅 [7] which is defined as the line integral along the line 𝐿 𝑖, 𝑗 ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1) = ∫ 𝐿 𝑖, 𝑗 𝑓 ( 𝑥 ) d 𝑥 , (6)where the line 𝐿 𝑖, 𝑗 is characterized by the sample position 𝑠 𝑖 and the projection angle 𝜃 𝑗 .In our specific use-case of extruded plastic profiles, the samples consist of only one materialand air. The density function 𝑓 in Eq. (6) is therefore piecwise constant with values either 0 or 𝛼 .For a given (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1) , the Radon transform of 𝑓 is therefore proportional to the sample thicknessalong the line integral 𝐿 𝑖, 𝑗 ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1) = ¯ 𝛼𝑑 𝑖, 𝑗 . (7)We can use Eq. (7) to rewrite Eq. (5) for the time domain THz-CT model 𝐸 𝑖, 𝑗 ( 𝑡 ) = 𝐸 ref (cid:18) 𝑡 − ¯ 𝑛 − 𝑛 𝑐𝛼 ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1)(cid:19) exp (cid:18) − ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1)(cid:19) , (8)where 𝐸 𝑖, 𝑗 denotes the final measurement signal in time-domain.In the following, we deduce two different approaches for mapping each measurement signal 𝐸 𝑖, 𝑗 ( 𝑡 ) to a scalar sinogram value 𝑔 𝑖, 𝑗 that acts as the basis for the image reconstruction approachesshown in section 4. In the conventional approach of X-CT, information about the sample is retrieved from thetransmitted intensity of X-Ray radiation 𝐼 𝑖, 𝑗 with respect to a reference intensity 𝐼 ref . Thesinogram values for X-CT are given according to 𝑔 𝑋𝑖, 𝑗 = − ln (cid:2) 𝐼 𝑖, 𝑗 / 𝐼 ref (cid:3) .Analogously, the energy of a THz pulse is defined as ∫ +∞−∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) d 𝑡 . By calculating theenergy from Eq. (8) we find 𝑔 𝐸𝑖, 𝑗 : = − ln +∞ ∫ −∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) d 𝑡 +∞ ∫ −∞ 𝐸 ref ( 𝑡 ) d 𝑡 = ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1) , (9)with 𝑔 𝐸𝑖, 𝑗 denoting the pulse energy sinogram values. .2. Pulse Phase Sinogram Mapping
The approach presented above extracts information about the sample based on the signal amplitude.Since the measurement setup is designed to measure the time evolution of the electric field of theTHz pulse, it is crucial to find a way to access the phase information as well. According to ourmodel in Eq. (8), the phase of the THz pulse is given by the time-shift term ¯ 𝑛 (cid:48) − 𝑛 𝑐𝛼 ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1) .Standard signal processing offers a straightforward way to calculate the time-shift, for instanceby means of peak-finder algorithms and thresholding. However, the resulting sinogram isnon-continuous, which can lead to artifacts in the reconstructed image. Furthermore, thisapproach does not offer a robust way of dealing with THz signals containing multiple peaks. Inpractice, these signals typically appear due to the finite cross-section of the THz beam (FWHM:1 mm) travelling through differently thick parts of the sample simultaneously. We have recentlydiscussed this effect in [17] by means of a more complex, non-linear imaging model.In this work, we therefore proceed with our more compact model of the mathematical THz rayand extract the time-shift term in Eq. (8) (cid:52) 𝑡 : = ¯ 𝑛 − 𝑛 𝑐𝛼 ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1) . (10)Calculating the quantity ∫ +∞−∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) 𝑡 d 𝑡 from Eq. (8) by substituting 𝑡 (cid:48) = 𝑡 − (cid:52) 𝑡 results in +∞ ∫ −∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) 𝑡 d 𝑡 = +∞ ∫ −∞ 𝐸 ref ( 𝑡 (cid:48) ) 𝑡 (cid:48) d 𝑡 (cid:48) + (cid:52) 𝑡 +∞ ∫ −∞ 𝐸 ref ( 𝑡 (cid:48) ) d 𝑡 (cid:48) exp (cid:2) ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1)(cid:3) . (11)The remaining integrals in Eq. (11) are independent of 𝑓 and calculated from the referencemeasurement. From Eq. (9) we find that dividing Eq. (11) by ∫ +∞−∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) d 𝑡 and rearrangingterms results in (cid:52) 𝑡 = +∞ ∫ −∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) 𝑡 d 𝑡 +∞ ∫ −∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) d 𝑡 − +∞ ∫ −∞ 𝐸 ref ( 𝑡 ) 𝑡 d 𝑡 +∞ ∫ −∞ 𝐸 ref ( 𝑡 ) d 𝑡 . (12)Reinserting the definition Eq. (10) back into Eq. (12) and rearranging terms leads to the pulsephase sinogram values 𝑔 𝑃𝑖, 𝑗 : = 𝑐𝛼 ¯ 𝑛 − 𝑛 +∞ ∫ −∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) 𝑡 d 𝑡 +∞ ∫ −∞ 𝐸 𝑖, 𝑗 ( 𝑡 ) d 𝑡 − +∞ ∫ −∞ 𝐸 ref ( 𝑡 ) 𝑡 d 𝑡 +∞ ∫ −∞ 𝐸 ref ( 𝑡 ) d 𝑡 = ( 𝑅 𝑓 ) (cid:0) 𝑠 𝑖 , 𝜃 𝑗 (cid:1) . (13)The latter term in Eq. (12) is calculated from the reference measurement and is therefore constant,while the material parameters cause a linear scaling of the sinogram values.The advantage of this procedure over standard peak finding algorithms is its robustness and theabsence of tuning parameters (e.g. thresholds). Furthermore, integrating over the THz signalshas a time-domain averaging effect, leading to a continuous sinogram and reduced SNR. Inaddition, it automatically takes into account the case of multiple peaks appearing in one THzsignal, as discussed above. . Image Reconstruction Approaches Both the pulse energy and the pulse phase mappings in Eq. (9) and (13) amount to a linear system ( 𝑅 𝑓 )( 𝑠 𝑖 , 𝜃 𝑗 ) = 𝑔 𝑖, 𝑗 , (14)for the unknown density 𝑓 , with right hand sides 𝑔 𝑖, 𝑗 = 𝑔 𝐸𝑖, 𝑗 and 𝑔 𝑖, 𝑗 = 𝑔 𝑃𝑖, 𝑗 , respectively. Since Eq.(14) is nothing else than a discretized version of the continuous Radon transform equation [7, 18] ( 𝑅 𝑓 )( 𝑠, 𝜃 ) = 𝑔 ( 𝑠, 𝜃 ) , (15)standard reconstruction methods for this equation can adapted to solve our THz-CT problem.Perhaps the most well-known of these methods is filtered back-projection, which is basedon the classic Radon inversion formula [7, 18]. The filtering step in this method acts as aregularization, which is necessary since the problem is in general ill-posed and unstable withrespect to unavoidable measurement noise in the data. Another popular approach for solvingEq. (15) is Tikhonov regularisation [11], which determines a stable approximation of 𝑓 byminimizing || 𝑅 𝑓 − 𝑔 || + 𝛽 || 𝑓 || , (16)where 𝛽 denotes a regularization parameter. The minimizer can be computed explicitly, in ourdiscrete case as the solution of a linear matrix-vector equation. Alternatively, a common iterativeapproach for solving (15) is given by Landweber iteration [11], which is defined by 𝑓 𝑘 + = 𝑓 𝑘 + 𝛾𝑅 ∗ ( 𝑔 − 𝑅 𝑓 𝑘 ) , where 𝛾 denotes a step-size parameter. Here, the choice of stopping index acts as regularization.For our THz-CT problem we focus on the above three reconstruction methods, i.e., on filteredback-projection, Tikhonov regularization, and Landweber iteration. For further details on thesemethods, as well as for other possible reconstruction methods we refer to [11, 18].
5. Results
In order to compare the outcomes of the two different sinogram mapping methods and reconstruc-tion approaches presented above, we comparatively walk through and analyze the results obtainedfor the sample in Fig. 3. Further on, more reconstructed sample cross-sections, retrieved withthe most promising combination of sinogram mapping and reconstruction approach, are shownbelow and in the additionally provided supplementary material (see Supplement 1).
The sample in Fig. 3 was scanned over a range of 0 - 359 ° and the raw measurement data waspreprocessed according to section 2.3.In section 3.1 and 3.2, two different approaches for extracting information from THz signalsare discussed. For both approaches, the calculated sinograms are shown in Fig. 4 (a, c). Thereconstructed cross-sections in Fig. 4 (b, d) were obtained by applying the classical filteredbackprojection algorithm to the sinograms. For image reconstruction, we have chosen the built-inMATLAB [19] function iradon() with default settings. In this subsection, the different image reconstruction approaches discussed in section 4 areapplied to the pulse phase sinogram in Fig. 4 (c). They were implemented in MATLAB usingthe AIR TOOLS II toolbox [20] for assembling the Radon transform matrix. Fig. 5 shows theresults of Tikhonov regularization (left) and Landweber iteration (right). These were obtained ig. 3. 3D printed sample for THz CT measurements: 3D representation of a circularsample with the letters "THZ" inscribed (left). Cross-section of the sample (right).Cross-hatched areas represent solid walls of PP. a) b)c) d)
Fig. 4. Comparison of sinogram mapping methods. (a) Pulse energy sinogram, (b) FBPof pulse energy sinogram, (c) Pulse phase sinogram, (d) FBP of pulse phase sinogram.The reconstructed image (d) obtained by pulse phase sinogram mapping (c) showsa significantly better representation of the sample geometry than conventional pulseenergy based FBP (b). ith the choice 𝛽 = 𝛾 = . · − for thestep-size in Landweber iteration, which was stopped after 70 iterations. These parameters wereoptimized manually for minimal noise and maximum image sharpness. b)a) Fig. 5. Comparison of image reconstruction approaches: (a) Tikhonov regularization,(b) Landweber iteration. Both reconstruction approaches yield similar results, withLandweber iteration (b) resulting in a more homogeneous representation of the letters
THZ . For quantitative validation of our THz-CT measurements, we designed and measured two furthersamples with more ordinary geometries. Due to the sample design, scattering artifacts in thereconstruction are to be expected only due to 3D-printed inhomogenities (e.g. layer stacks,joints).While the geometry of the first object in Fig. 6 (top) is parameterized by a spiral with linearlyincreasing radius, the second object (bottom) is highly symmetrical and consists of planarelements only.Both images were reconstructed using Landweber iteration with the same parameters as insection 5.2. As a reference measurement for both objects, the nominal wall thickness (2 mm)was verified within minor deviations of ± . )b) Fig. 6. Quantitative evaluation of reconstructed images: Sample cross-sections (left),Quantitative evaluation of a curved sample (a), Quantitative evaluation of a planarsample (b). The blue hatched area in the sample cross-section represents the samplewalls. Blue lines in the reconstructed image represent the calculated sample contourand green lines mark the identified sample boundaries using the Hough transform.Sample a) shows a gradually increasing deviation of the wall thickness towards itscenter of rotation. Overall good agreement within the accuracy of the 3d printer wasachieved for sample b).
6. Discussion
We begin this section with a qualitative comparison of the sinogram mapping methods in Fig. 4.In both cases (b, d), the cylinder around the letters
THZ is well identifiable. Since the cylinder is3D printed from multiple circles next to and on top of each other, an inhomogenious joint/stackof layers is produced. This joint can be identified in both reconstructed images at the position ( 𝑦 | 𝑧 ) ≈ (
55 mm |
45 mm ) .Generally, the pulse energy sinogram (a) and reconstruction (b) yields more pronounced edgesand inhomogenities than the pulse phase sinogram. We presume that this is caused by the stronginfluence of scattering effects on small sample structures and edges, effectively reducing thedetected signal amplitude. On the other hand, the pulse phase mapping approach eliminates theinfluence of the absolute signal amplitude by normalizing each measured signal with respect toits energy (see Eq. 13). This property effectively suppresses the presence of large peaks in thereconstructed image (d), and leads to a better representation of the sample geometry in general.Furthermore, it is worth noting that the thickness of the letters in Fig. 4 (d) is significantly largerhan the cylinders wall thickness, despite both being printed with a thickness of 2 . ± . THZ we observe a much more homogeneous structure andless artifacts for Tikhonov regularization and Landweber iteration. After manual parameteroptimization, both reconstruction approaches seem to outperform the classical approach offiltered backprojection.We have further examined two more samples for quantitative testing of the accuracy of theproposed THz-CT approach. In Fig. 6 (a) we have observed a locally increasing wall thickness ofup to 2.6 mm. The tendency indicates that the amount of surrounding sample structure layersartificially increases the wall thickness in the reconstructed image. Throughout this work, thiseffect was consistently and explicitly observed for all samples with curved surrounding structures(see Supplement 1). Thus, we strongly assume that this artifact is predominantly caused bythe wave propagation effect of refraction. So far, we have not accounted for this effect in ourmodel. Nevertheless, the spiral structure is qualitatively well identifiable and we therefore planto use the current approach as the basis for further improvements in the image reconstruction ofcurved objects. In the industrially more relevant case of a planar profile (see Fig. 6, b), excellentagreement with the nominal wall thickness of 2 . ± .
7. Conclusion
In this work, we have presented two computationally efficient and robust approaches for performingTHz-CT measurements on plastic profiles by means of a fast THz-TDS system. While the firstapproach, in the style of classical X-CT, utilizes the amplitude information of the measuredsignals for image reconstruction, the second approach makes use of the availability of phaseinformation in the THz-TDS measurement. We have experimentally shown that THz-CT based onthe signal amplitude is prone to wave scattering artifacts. Structures with sizes in the dimensionof the wavelength or sharp edges lead to extended peaks in the reconstructed image, potentiallyovershadowing other structures. In contrast, the influence of scattering is mostly eliminated inthe second approach, leading to overall better representations of the sample geometry.We have further compared the performance of multiple image reconstruction methods, namelyfiltered backprojection, Tikhonov regularization, and Landweber iteration. It was empiricallyfound that Tikhonov regularization and Landweber iteration yield more homogeneous imagesthan filtered backprojection.As a last point, the quantitative performance of our THz-CT system was evaluated using theexample of a curved and a planar sample (see Fig. 6). Local deviations were found for thespecial case of the curved sample. The planar sample could be reconstructed correctly withinan accuracy of ± . Funding.
Österreichische Forschungsförderungsgesellschaft (871974); European Commission (777222);ustrian Science Fund (F6805-N36).
Acknowledgments.
Part of this work was performed under the scope of the COMET program within theresearch project "Photonic Sensing for Smarter Processes". This program is promoted by BMK, BMDW,the federal state of Upper Austria and the federal state of Styria, represented by SFG. The research was alsosupported by the strategic program "Innovatives OÖ 2010 plus" by the Upper Austrian Government.
Disclosures.
The authors declare no conflicts of interest.
Data availability.
Data underlying the results presented in this paper are not publicly available at thistime but may be obtained from the authors upon reasonable request.
Supplemental document.
See Supplement 1 for supporting content.
References
1. R. J. B. Dietz, N. Vieweg, T. Puppe, A. Zach, B. Globisch, T. Göbel, P. Leisching, and M. Schell, “All fiber-coupledTHz-TDS system with kHz measurement rate based on electronically controlled optical sampling,” Opt. letters ,6482–5 (2014).2. S. Krimi, J. Klier, J. Jonuscheit, G. von Freymann, R. Urbansky, and R. Beigang, “Highly accurate thicknessmeasurement of multi-layered automotive paints using terahertz technology,” Appl. Phys. Lett. , 021105 (2016).3. K. B. Cooper, R. J. Dengler, G. Chattopadhyay, E. Schlecht, J. Gill, A. Skalare, I. Mehdi, and P. H. Siegel, “Ahigh-resolution imaging radar at 580 ghz,” IEEE Microw. Wirel. Components Lett. , 64–66 (2008).4. B. Recur, J. P. Guillet, L. Bassel, C. Fragnol, I. Manek-Hönninger, J. Delagnes, W. Benharbone, P. Desbarats,J. Domenger, and P. Mounaix, “Terahertz radiation for tomographic inspection,” Opt. Eng. , 1–8 (2012).5. S. Mukherjee and J. Federici, “Study of structural defects inside natural cork by pulsed terahertz tomography,” in (2011).6. J. Neu and C. A. Schmuttenmaer, “Tutorial: An introduction to terahertz time domain spectroscopy (THz-TDS),” J.Appl. Phys. , 231101 (2018).7. F. Natterer, The Mathematics of Computerized Tomography (Society for Industrial and Applied Mathematics, 2001).8. A. Wald and T. Schuster, “Tomographic terahertz imaging using sequential subspace optimization,” in
New Trends inParameter Identification for Mathematical Models,
B. Hofmann, A. Leitão, and J. Zubelli, eds. (Birkhäuser, Cham,2018).9. J. Tepe, T. Schuster, and B. Littau, “A modified algebraic reconstruction technique taking refraction into account withan application in terahertz tomography,” Inverse Probl. Sci. Eng. (2016).10. M. L. Louarn, C. Vérinaud, V. Korkiakoski, N. Hubin, and E. Marchetti, “Adaptive optics simulations for theEuropean Extremely Large Telescope,” in in Advances in Adaptive Optics II no. 6272 in Proc. SPIE, (2006), p.627234.11. H. W. Engl, M. Hanke, and A. Neubauer,
Regularization of inverse problems. (Dordrecht: Kluwer AcademicPublishers, 1996).12. I. Pupeza, R. Wilk, and M. Koch, “Highly accurate optical material parameter determination with thz time-domainspectroscopy,” Opt. Express , 4335–4350 (2007).13. M. Scheller, C. Jansen, and M. Koch, “Analyzing sub-100- µ m samples with transmission terahertz time domainspectroscopy,” Opt. Commun. , 1304–1306 (2009).14. Y. Jin, G. Kim, and S. Jeon, “Terahertz Dielectric Properties of Polymers,” J. Korean Phys. Soc. , 513–517 (2006).15. B. A. Shenoi, Finite Impulse Response Filters (John Wiley & Sons, Ltd, 2005), chap. 5, pp. 249–302.16. D. J. Griffiths,
Introduction to electrodynamics (Pearson, 2014), 4th ed.17. S. Hubmer, A. Ploier, R. Ramlau, P. Fosodeder, and S. van Frank, “A mathematical approach towards thz tomographyfor non-destructive imaging,” submitted, available from: https://arxiv.org/abs/2010.14938 (2021).18. A. K. Louis,
Inverse und schlecht gestellte Probleme , Teubner Studienbücher Mathematik (Vieweg+Teubner Verlag,1989).19. MATLAB, (The MathWorks Inc., Natick, Massachusetts, 2016).20. P. C. Hansen and J. Jorgensen, “Air tools ii: algebraic iterative reconstruction methods, improved implementation,”Numer. Algorithms79