Terahertz (THz) waves have high penetration through non-metallic and non-polar materials, and they have low photon energy, so are safer for biological targets and human operators. Because of these special properties, both the two-dimensional (2D) and three-dimensional (3D) THz imaging techniques have gained widespread attention in biomedical fields^[1] and nondestructive testing^[2], among others^[3–7]. The 3D THz imaging can reveal the inner structures of the object, and, especially, the refractive index (RI) distributions usually offer higher contrast and more plentiful information.

A number of methods for 3D imaging with THz radiation have also been proposed and demonstrated. One of main modalities is the THz computed tomography (CT)^[5], which reconstructs the 3D absorption coefficient and RI distributions of the sample from the measured intensity and phase distributions of various projection angles. THz CT is based on a geometrical straight-line model of the propagation of radiation, which assumes that diffraction effects and Fresnel losses can be neglected. The Fourier slice theorem is the central theorem, which leads to the well-known 3D reconstruction algorithm, the filtered backprojection algorithm (FBPJ). The other iterative reconstruction algorithms, such as the simultaneous algebraic reconstruction technique (SART)^[6], were developed. The improved THz CT systems with high scanning speed design, extended depth of field^[7], and RI matching method^[8] have been proposed. However, the THz CT is not proper for 3D imaging of the sample with a fine structure, like high-accuracy 3D RI measurement, and the 3D reconstruction results would be seriously degraded^[9,10] because of the diffraction effect caused by the large wavelength of the THz wave. For conventional CT systems, the resolution is related to the size of the focused beam, which is larger than 1 mm in THz waves. Also, it is hard to improve the imaging lateral resolution.

Diffraction tomography (DT) is one of the common tomography methods that accounts for the diffraction effect^[11,12] and can determine simultaneously the 3D absorption coefficient and RI distributions of the sample with high accuracy by using the measurement of the diffracted wave field. In the visible light domain, presently, the optical DT (ODT) has become an important 3D imaging method for non-invasive observation of biomedical samples^[13,14] and structural inspection of micro-optical elements^[15,16]. In the THz domain, Wang et al.^[17,18] in 2003 reported the first implementation of DT by employing a THz time domain spectroscopy (TDS) system. In their work, the sample was an object composed of three rectangular polyethylene cylinders. The 2D THz pattern was formed on an electro-optic ZnTe crystal, which was transferred by the visible probe light and then recorded by a CCD detector. The Fourier diffraction theorem (FDT) and interpolation in the spatial frequency domain were used to reconstruct the spatial object function. This 3D reconstruction algorithm is denoted as a frequency domain interpolation (FDI) algorithm, which is often used. However, the signal complex interpolation for the 3D spatial spectrum data led to artifacts in the reconstructed RI distributions. Thus, the reconstructed results were not very good with noise background. Besides, the system was rather complicated, and it needed to precisely control the CCD image acquisition and THz beam on–off time sequence. Moreover, the above work also assumes that the object is dispersionless. Therefore, in order to meet the practical requirements for 3D imaging, the optimized THz DT with a simple configuration still needs further development, including the direct intensity recording of a THz wave by the 2D array detectors, the use of single-frequency continuous-wave (CW) THz sources, and more effective 3D reconstruction algorithms. Recently, the progress on the wide-field THz coherent imaging techniques makes the exploration very possible, such as THz digital holography^[19,20] and ptychography^[21], which can give, simultaneously, the 2D quantitative amplitude and phase distributions of the sample.

In this Letter, we proposed a CW-THz DT system by using a coherent single-frequency THz laser and the array detector to directly record the off-axis digital hologram. The setup is simple and easy to operate. By rotating the object, the 2D scattered fields of the object at each rotation angle can be obtained by the digital holographic reconstruction method. With regard to the 3D reconstruction algorithm of the THz DT, to the best of our knowledge, the filtered backpropagation (FBPP) method is first introduced to produce the volume data in the space domain and to achieve high-quality 3D RI measurement. Compared with the regular FDI algorithm, there is no need to do complex interpolation in the Fourier domain, and the experimental results demonstrate the validity of the method with high fidelity.

On the assumption of the scalar wave diffraction theorem, consider a monochromatic plane wave $u_0(\stackrel{\rightharpoonup }{r})=\exp ({\stackrel{\rightharpoonup }{K}}_m·\stackrel{\rightharpoonup }{r})$ incident upon a scattering object, whose physical quantity is termed the scattering potential $f(\stackrel{\rightharpoonup }{r})$. The resultant total field $u(\stackrel{\rightharpoonup }{r})$ [incident field $u_0(\stackrel{\rightharpoonup }{r})$ plus scattered field $u_s(\stackrel{\rightharpoonup }{r})$] satisfies the inhomogeneous Helmholtz equation as^[11]$$({\nabla }^2+k_0^2)u_2(\stackrel{\rightharpoonup }{r})=-f(\stackrel{\rightharpoonup }{r})u(\stackrel{\rightharpoonup }{r}),$$(1)and the scattering potential $f(\stackrel{\rightharpoonup }{r})$ is defined by$$f(\stackrel{\rightharpoonup }{r})=k_0^2\left[(\frac{n(\stackrel{\rightharpoonup }{r})}{n_0}{)}^2-1\right],$$(2)where $\stackrel{\rightharpoonup }{r}=(x,y,z)$ is the 3D spatial vector at the object space, ${\stackrel{\rightharpoonup }{K}}_m=(k_{mx},k_{my},k_{mz})$ is the spatial frequency vector of the incident fields, $k_0=|{\stackrel{\rightharpoonup }{K}}_m|$ is the wave number in a vacuum, $n_0\approx 1$ is the RI value of air, and $n(\stackrel{\rightharpoonup }{r})$ is the RI distributions of the object.

In order to solve Eq. (1), we can make use of Green’s function and apply the Rytov approximation, so the scattered field can be obtained as^[22]$$u_s(\stackrel{\rightharpoonup }{r})\approx u_{\mathrm{Rytov}}(\stackrel{\rightharpoonup }{r})=u_0(\stackrel{\rightharpoonup }{r})\ln [u(\stackrel{\rightharpoonup }{r})/u_0(\stackrel{\rightharpoonup }{r})].$$(3)

Then, a direct relation between the Fourier transform (FT) of the scattered field and the FT of the scattering potential of the object can be deduced, which is commonly called the FDT^[23]:$$\tilde{F}({\stackrel{\rightharpoonup }{K}}_D-{\stackrel{\rightharpoonup }{K}}_m)=-\sqrt{\frac{2}{\pi }}i{k}_{Dz}{\tilde{U}}_s(k_{Dx},k_{Dy};{\varphi }_0),$$(4)where $\tilde{F}({\stackrel{\rightharpoonup }{K}}_D-{\stackrel{\rightharpoonup }{K}}_m)$ is the 3D FT of the scattering potential $f(\stackrel{\rightharpoonup }{r})$, ${\tilde{U}}_s(k_{Dx},k_{Dy};{\varphi }_0)$ is the 2D FT of the scattered wave $u_s(x,y;{\varphi }_0)$ with the rotation angle ${\varphi }_0$, and ${\stackrel{\rightharpoonup }{K}}_D=(k_{Dx},k_{Dy},k_{Dz})$ is the spatial frequency vector of scattered field, here, $k_{Dx}^2+k_{Dy}^2+k_{Dz}^2=k_0^2$. It states that the 2D FT of the filtered data function is equal to the 3D FT of $f(\stackrel{\rightharpoonup }{r})$ evaluated over a semispherical surface oriented at an angle ${\varphi }_0$ in the 3D Fourier space of $f(\stackrel{\rightharpoonup }{r})$.

To reconstruct exact $f(\stackrel{\rightharpoonup }{r})$ directly in the space domain^[23–25], Devaney developed the FBPP algorithm based on Eq. (4) when the frequency domain is represented, by the polar coordinate, and the rotation coordinate transformation is performed. For simplicity, we only consider the rotation of the sample along one axis, which is denoted as the $y$ axis here. The scattering potential can be reconstructed as$$\begin{gathered} f(\stackrel{\rightharpoonup }{r})=\frac{-i{k}_0}{4{\pi }^2}{\int }_0^{2\pi }{\int }_{-k_0}^{k_0}{\int }_{-k_0}^{k_0}|k_{D_x}|{\tilde{U}}_s(k_{D_x},k_{Dy};{\varphi }_0) \\ \times \exp \left[i{z}_{{\varphi }_0}(\sqrt{k_0^2-k_{Dx}^2-k_{Dy}^2}-k_0)\right] \\ \times \exp [i(k_{Dx}{x}_{{\varphi }_0}+k_{Dy}{y}_{{\varphi }_0})]\mathrm{d}{k}_{Dx}\mathrm{d}{k}_{Dy}\mathrm{d}{\varphi }_0, \end{gathered}$$(5)where $$x_{{\varphi }_0}=x\cos {\varphi }_0+z\sin {\varphi }_0,$$(6a)$$y_{{\varphi }_0}=y,$$(6b)and $$z_{{\varphi }_0}=-x\sin {\varphi }_0+z\cos {\varphi }_0.$$(6c)

In Eq. (5), there are two exponential terms inside the integral formula. The first one is a transfer function, which is depth-dependent due to the parameter $z_{\varphi 0}$. The second one is a factor to make the whole integral form a 2D inverse FT over $k_{Dx}$ and $k_{Dy}$. The whole integral formula can be interpreted as propagating the scattered field along the $z_{\varphi 0}$ direction, which is similar to the diffraction propagation. In order to implement the numerical calculation of Eq. (5), the scattered field can be propagated into $N$ sections in the space containing the object, and the interval between adjacent sections is denoted as $\mathrm{\Delta }{z}_{\varphi 0}$, which can be flexibly set. Thus, volume data in the space domain is produced, and the scale of the depth is $N·\mathrm{\Delta }{z}_{\varphi 0}$. Finally, the transformation of coordinates, illustrated as Eq. (6), is performed for 3D volumes, which are then added to the reconstruction volume in the sum over all rotation angles ${\varphi }_0$. Correspondingly, the lateral and axial resolutions of the DT system are related to the resolution of the complex amplitude distribution of a single projection and the number of projections^[26]. Under the proper number of projections, the former factor usually plays an important role.

It is noted that the FBPP algorithm generally requires scattered data measured from view angles in [$0,2\pi$] for exact reconstruction of a complex-valued object function. Compared to the direct linear inversion (the FDI algorithm) based on Eq. (4), there is no need to do the complex interpolation processing in the frequency domain, which will cause large computational errors and produce artifacts. The error-prone frequency interpolation can be avoided by applying the FBPP algorithm, which is theoretically equivalent to the FDI.

Practically, in order to achieve the complex diffracted fields (including the amplitude and phase distributions) of the object for one illumination angle ${\varphi }_0$, the off-axis digital holographic or the phase-shifting interferometry method can be employed. The off-axis digital holographic method is used in this Letter. The filtering frequency spectrum is applied to the digital holograms to obtain $u(x,y;{\varphi }_0)$ and $u_0(x,y;{\varphi }_0)$, which are with and without the object. The scattered field can be calculated by Eq. (3).

An experimental CW THz DT system based on off-axis Fresnel digital holography was built as depicted in Fig. 1 to collect the scattered field of the sample with different rotation angles. The THz source was an optically pumped far-IR gas laser (295-FIRL; Edinburgh Instruments Ltd., UK) with a central wavelength of 118.83 µm (2.52 THz) and a maximum power of 500 mW. The emitted THz beam was expanded and collimated by two off-axis parabolic mirrors (PM1, $f_1=50.8\mathrm{mm}$; PM2, $f_2=152.4\mathrm{mm}$) and then divided into the object beam and the reference beam by a beam splitter (BS). The object was attached to a metal holder for fixing on a rotational stage (PRMTZ8/M; Thorlabs Inc., USA). The interference fringes generated by the two beams formed the hologram, which was recorded by a pyroelectric detector (Pyrocam IV, pixel size $320\times 320$, pixel pitch $80\mathrm{\mu m}\times 80\mathrm{\mu m}$). The distance from the object to the detector was $d=62.9\mathrm{mm}$. To enhance the contrast of the interference fringes, 500 frames were recorded at a chopping frequency of 50 Hz and accumulated via Gaussian fitting. In this experimental setup, the theoretical resolutions are about 0.3 mm.

To validate the proposed method, we first used a single polystyrene (PS) foam sphere with a diameter of 7.44 mm. The RI of the PS foam was $\sim 1.0169$ at 2.52 THz, as measured by the THz-TDS system (TAS7400SU; Advantest, Japan). In the experiments, the sample was rotated gradually by 360° at angular intervals of 3°, and thus 120 holograms with the object were recorded successively, followed by one background hologram without the object. Figures 2(a1) and 2(b1) show the holograms with and without the object, respectively, at the rotation angle of 0°, and Figs. 2(c1) and 2(d1) show the reconstructed amplitude and wrapped phase images of Fig. 2(a1). The least-squares method is used to obtain the unwrapped phase image^[27], then a double exposure is used to compensate for phase aberration^[28], and the result is shown in Fig. 2(a2). Figures 2(a2)–2(d2) show the reconstructed phase images at 0°, 45°, 90°, and 180°, respectively, and Fig. 2(e1) shows the phase profiles of the dotted lines in Figs. 2(a2)–2(d2). As can be seen, the maximum phase values of the different angles are shifted along the X axis. This misalignment exists because it is difficult to precisely adjust the gravity center of the sample to the rotational axis of the stage. To compensate for this mechanical error, a cross-correlation coefficient was used to align the phase profiles. After alignment, the maximum phase values are nearly in the same position along the X axis, as shown in Fig. 2(e2). The mechanical error mentioned above also results in axial shifting of the sample when rotating, so the reconstructed complex images of digital holography are defocused to different extents. To improve this, the auto-focusing method based on the Brenner function as the criterion function is applied to make them with a sharp boundary at each rotation angle^[29].

The reconstructed complex amplitudes by digital holography with various angles were then processed to obtain the scattered field. The scattering potential distributions of the sample were reconstructed by using Eq. (5), where $\mathrm{\Delta }{z}_{\varphi 0}$ is set to be 80 µm, and N is equal to 320. The RI distributions were finally achieved by using Eq. (2). In our experiment, the sample was supported by a metal rod whose large RI led to incorrect reconstructed values in these areas, so the support-rod part has been removed from the reconstructed 3D RI results below.

Figures 3(a1)–3(c1) are the reconstructed results by the FDI method, showing the obtained RI distributions at cross sections $x–y$, $y–z$, and $x–z$, respectively. It is seen that there is some artifacts error inside, especially in the background region, and the fluctuation of the value is relatively serious. The reconstructed RI distributions by the FBPP method are shown in Figs. 3(a2)–3(c2), which show much-improved fidelity of the reconstructed RI tomograms. Figure 3(d) shows the differences of the RI profiles of the PS foam sphere between the results obtained by the FBPP method and the FDI method. The red line represents the ideal RI value of the PS foam material, which is measured by the THz-TDS system. As can be seen, the RI value obtained by the FBPP method is more accurate, where the average value is $\sim 1.0153$, and the error is only 0.16%, thereby verifying the effectiveness of the proposed method. Furthermore, the average diameter of the reconstructed foam sphere is $\sim 7.44\mathrm{mm}$, which corresponds to a reconstruction error of 0.6% compared with the true value. We can also analyze the reconstruction quality of the FBPP and FDI methods from the view of 3D Fourier space, as shown in Fig. 4. Figures 4(a) and 4(b) are the amplitude of the spectrum distribution of FDI and FBPP on the logarithmic scale along the $f_x–f_y,f_x–f_z$, and $f_y–f_z$ cross sections, respectively. For the FDI method, the number of missing spectrum points produced by the reconstruction algorithm is more than that of FBPP, as indicated by the green arrows. It leads to the reconstructed RI distribution being worse, as shown in Figs. 3(a1)–3(c1). Meanwhile, both FBPP and FDI have the same missing spectrum points caused by the “missing apple core” problem, as indicated by the red arrows.

To demonstrate further the applicability of THz DT on non-axisymmetric samples, two glued PS foam spheres were placed horizontally on the rotary stage, and then the sample was rotated around the gravity center of one sphere to obtain the scattered field with various rotation angles. The reconstructed RI distributions of the foam spheres at cross sections $x–y,y–z$, and $x–z$ and the 3D distributions are shown in Figs. 5(a)–5(d), respectively. From Visualization 1, the 3D RI distributions are presented visually. The average RI value of the two reconstructed foam spheres is $\sim 1.0126$, and the difference is 0.42% compared with the value obtained by the THz-TDS system. Note that the RI distributions are not as uniform compared with the previous single foam sphere. There are some fringe patterns. Because multiple scattering exists when the THz beam propagates through the two foam spheres, the Rytov approximation is not well satisfied, and the reconstruction error is increased. To solve this, one possible way is to build the forward propagation model including the multiple scattering effects and apply an iterative DT reconstruction algorithm to improve the quality of 3D reconstructed RI distributions^[30,31].

In summary, we have realized CW THz DT combined with digital holography, and the digital holograms are recorded directly by the array detector at various angles through rotating the samples. The configuration is simple and easy to operate. With regard to the 3D reconstruction algorithm, the FBPP algorithm is adopted to achieve 3D RI distributions of the PS foam spheres. The average RI value has only 0.16% difference from the RI value measured by the THz-TDS system. The reconstructed results have high fidelity compared with results obtained by the FDI algorithm. It verifies the feasibility of the proposed method.

Note that in the reconstructed RI tomograms in Figs. 3 and 5, some values are quite different from their surroundings along the Y axis, which is the result of the “missing apple core” problem. It can be improved by the iterative algorithm with non-negative constraints. Furthermore, to promote THz DT further, RI matching or the multiple scattering non-linear model must be applied to reconstruct high-RI samples. Our view is that THz DT can be an effective method for non-destructive testing and quantitative measurement of the complex RI of complex samples in the future, which can be combined with other THz wide-field phase-contrast imaging methods.