Breaking the Limitations with Sparse Inputs by Variational Frameworks (BLIss) in Terahertz Super-Resolution 3D Reconstruction

Data acquisition, image processing, and image quality are the long-lasting issues for terahertz (THz) 3D reconstructed imaging. Existing methods are primarily designed for 2D scenarios, given the challenges associated with obtaining super-resolution (SR) data and the absence of an efficient SR 3D reconstruction framework in conventional computed tomography (CT). Here, we demonstrate BLIss, a new approach for THz SR 3D reconstruction with sparse 2D data input. BLIss seamlessly integrates conventional CT techniques and variational framework with the core of the adapted Euler-Elastica-based model. The quantitative 3D image evaluation metrics, including the standard deviation of Gaussian, mean curvatures, and the multi-scale structural similarity index measure (MS-SSIM), validate the superior smoothness and fidelity achieved with our variational framework approach compared with conventional THz CT modal. Beyond its contributions to advancing THz SR 3D reconstruction, BLIss demonstrates potential applicability in other imaging modalities, such as X-ray and MRI. This suggests extensive impacts on the broader field of imaging applications.


Introduction
Gaining increasing attention in recent years, terahertz (THz) imaging has emerged its worth as a versatile tool for various applications such as non-destructive inspection [1][2][3], security screening [4][5][6][7], and non-contact testing [8][9][10][11].The unique properties of THz radiation, including its low photon energy, make it suitable for safe use in bio-related scenarios [12,13].Additionally, its capability to penetrate optically opaque materials allows extracting internal geometric and material information of tested objects [14][15][16][17][18]. THz imaging includes a broad spectrum of modalities, with pulsed THz imaging emerging as one of the most prevalent and widely employed techniques in this field.This is often associated with a THz time-domain spectroscopy (THz-TDS) system [13,[19][20][21][22][23], enabling the acquisition of ultrafast responses within the domain of light-matter interactions across temporal, spatial, and spectral dimensions.Recently, it has paved the way for a range of applications, including analysis of conformational behaviors [24], ultrafast molecular dynamics [25], dielectric responses [26], and electrical properties [27].
Leveraging these advantages, pulsed THz imaging finds practical use in visualizing the internal structure of 3D objects.Among the frequently employed pulsed THz imaging techniques, computed tomography (CT) represents an effective computational method.This technique provides cross-sectional images (slices) reconstructed from a sequence of projected images taken at various projection angles by inverse Radon transform R (IRT) [9,28].While X-ray CT has found extensive use in various real-world applications, its ionizing characteristics and high penetration capacity inherently restrict its utility, particularly in the visualization of composite materials, soft dielectrics, or liquids.Furthermore, the limited differentiation of material responses within the X-ray spectrum poses challenges in identifying the chemical components present within objects [9,29].In light of these limitations, the introduction of pulsed THz CT imaging, first proposed in 2002 [30], emerged as a promising supplement to traditional X-ray CT.
While pulsed THz CT imaging exhibits considerable potential as a supplement with nonionizing characteristics and sensitive material differentiation, two primary concerns are (i) the protracted process of data acquisition with limited scanning resolution by THz-TDS; and (ii) Fraunhofer diffraction during the THz wave propagation through the object [30,31].For clarification, the term resolution typically refers to the number of pixels in digital imaging rather than diffraction-limited resolution.This is because all datasets are collected by THz-TDS and processed under the pixel-based scanning resolution.Concerning (i), THz-TDS has historically relied on mechanical time-delay scanning.However, non-mechanical timedomain sampling methods, such as asynchronous optical sampling (ASOPS) introduced in 2005 [32] and electronically controlled optical sampling (ECOPS) introduced in 2010 [33], bring a promising alternative.These techniques can significantly reduce the data acquisition time for each single-cycle THz pulse, shifting from hour/minute scales to milliseconds.Despite advancements, the implementation of ASOPS and ECOPS requires costly or bulky femtosecond lasers.Moreover, projection data acquisition for THz 3D CT imaging still involves mechanical operations, such as raster scanning and object rotation, which constrain the resolution of projection image due to scanning resolution, impacting subsequent 3D reconstruction.Recently, a highthroughput alternative is the use of terahertz focal-plane arrays (THz-FPAs) [34].These detector arrays, designed based on field-effect transistors [35], microbolometers [36], and plasmonic nanoantennas [37,38], offer the potential to bypass mechanical operations.Nevertheless, the shift towards employing THz-FPAs introduces a unique set of challenges.These arrays require scalable high-speed electronics, a considerable increase in computational power, and increased system complexity, all of which depend on the number of THz detectors in use, and notably, this transition could be cost-intensive.Besides, while THz-FPAs can speed up data acquisition, the quality of reconstructed 3D images would be affected by the uniformity of each device and, most importantly, the concern (ii) related to diffraction effects with noise.
Addressing these issues without costly upgrades to the imaging system, algorithm-based methods show promise in directly enhancing low-resolution (LR) images into super-resolution (SR) images, tackling concerns (i) and (ii).For instance, by using prior knowledge of the degradation based on the specific point spread function (PSF) measured [39][40][41] or simulated [42][43][44] from the THz imaging system, a variety of SR algorithms have been explored, including Richardson-Lucy algorithm [45] and PSF deconvolution [46].These methods have demonstrated success in improving spatial resolution and reducing the diffractive limit, although challenges in rendering image edge and texture details remain.Nowadays, deep-learning-based SR methods have been increasingly recognized as a prominent approach.By establishing an end-to-end mapping correlation between LR and HR images, they have shown advantages in comprehensive learning ability and the overall effectiveness of image reconstruction after training.Examples include the convolutional neural network (CNN) with its variants, the very deep super-resolution (VDSR) method, the super-resolution network for multiple degradations (SRMD), and generative adversarial networks (GAN) [14,16,17,[47][48][49][50][51][52][53][54].However, the application of deep-learningbased methods in THz single-image super-resolution (SISR) faces several challenges, including the complex THz image degradation process, limited efficiency in deblurring and denoising procedures, and the issue of gradient vanishing during deep network training.In a nutshell, both algorithm-based and deep-learning-based methods demonstrate their potential in improving the spatial resolution of THz imaging, but their primary focus on 2D image processing often overlooks the global characteristics of the 3D post-reconstruction.Moreover, the need for sufficient datasets and costly workstations with GPUs for model training presents another challenge.
In this work, we focus on improving 3D THz reconstruction by integrating variational frameworks with conventional CT methods.Here, we introduce BLIss (as delineated in Fig. 1) with the aim of facilitating rapid and smooth SR 3D reconstruction, even when working with a limited pool of LR images.In detail, this framework integrates traditional CT techniques, including signal processing (e.g., projections, sinograms, and slices by IRT R), with the underpinning mathematical principles of variational formulations, i.e., Willmore-based formulation and Euler-Elastica-based formulation.To assess the quality and fidelity of the reconstructed results by our BLIss, the conventional evaluation metrics, e.g.mean-squared error (MSE), peak signal-to-noise ratio (PSNR), and structural similarity index (SSIM), which are commonly employed in 2D image assessment, are not sufficient due to the inconsistency of data type in 3D.To address this, we introduce an approach that quantifies the global smoothness stability of 3D surface quality using concepts from discrete geometry.This approach involves calculating the standard deviation of Gaussian and mean curvatures for the triangular meshes of the reconstructed surfaces.In addition, we evaluate the fidelity of reconstruction when using limited input data.To achieve this, we extend our analysis to incorporate the calculation of the multi-scale structural similarity index measure (MS-SSIM), a 3D structural similarity metric that builds upon the foundation of the 2D SSIM.Our proposed variational framework effectively addresses our concerns (i) and (ii), presenting a pathway to rapid and high-precision 3D THz imaging.This advancement holds significant promise for non-destructive sensing and inspection applications.

Experimental Setup for Data Acquisition
The experimental setup is shown in Fig. 2(a).The THz 3D imaging system records time-resolved THz signals for each voxel using an ASOPS THz-TDS system (TERA ASOPS, MenloSystems).
The pulsed THz wave is generated by a THz photoconductive antenna (PCA) emitter and is directed through four plane-convex polypropylene (PP) THz lenses, which focus it onto the sample with a beam size of 1.5 mm before refocusing it onto a THz PCA detector.The photocurrent signal is subsequently amplified by a transimpedance amplifier (TIA) and digitized by an analog-to-digital converter (ADC) processing unit.The data collection process accelerates due to the usage of two asynchronous femtosecond (fs) Er-doped fiber lasers in the ASOPS system.The repetition rate of the pump laser is 100 MHz + 50 Hz (  pump ), the probe laser operates at 100 MHz (  probe ), and the fixed offset frequency (Δ  =  pump −  probe ) is 50 Hz.The sampling rate ( sr ) is 10 MHz, which allows for a data acquisition time of approximately 20 ms per pulse with a temporal resolution of 50 fs.Each THz signal pulse in the dataset corresponds to 5000 sampling points per voxel and spans 250 ps, as instantiated in Fig. 2

Conventional THz CT Imaging
Conventional techniques for THz 3D CT imaging usually involve a series of steps to construct 3D images from the collected THz dataset.The initial step in the THz 3D CT imaging process is the acquisition of projection data.In this step, the Time-MAX projection technique is one of the most fundamental methods in THz 3D CT imaging.Time-MAX involves extracting the maximum field strength from each time-resolved THz signal trace, which is then used to represent the projection of the scanned object at a particular projection angle [14,15].This method effectively reduces the temporal dimensionality of the THz signals to a 2D projection image (see Fig. 3(a)).
The data obtained from Time-MAX is primarily associated with object thickness, surface contour, and material contrast, providing potential information into the spatial distribution of composite materials within the objects.Next, a sinogram (see Fig. 3(b)) is generated by stacking all the 2D projection images at different rotation angles in a 3D array.The sinogram visualizes the internal structure of the object from multiple angles and forms the basis for image reconstruction.
The horizontal axis represents the spatial location along the projection, and the vertical axis corresponds to the projection angle.The value at each pixel in the sinogram indicates the THz signal intensity at the corresponding location and angle.The final step in the conventional THz CT imaging process is the inverse Radon transform (IRT) R [15,55].This mathematical procedure reconstructs an image  (x) from its projections where the Radon transform R  (, ) is the integral of  along the line   ,  = {x : ⟨x, n  ⟩ = } whose distance from the origin is  and whose direction is orthogonal to the unit vector n  = (cos , sin ) with the convolution kernel ℎ satisfied ĥ() = ||.Applied to the THz field, IRT takes the sinogram as input and reconstructs a 2D cross-section (slice) of the object (see Fig. 3(c)).The value of each reconstructed pixel represents the normalized attenuated intensity, rescaled to the interval [0, 1] based on the minimum and maximum across all pixels.
Then, a 3D volume image of the object can be obtained (see Fig. 4(a)) after stacking all 2D slices along the vertical direction and selecting an appropriate threshold.Here, the threshold can be determined based on characteristics of the histogram, such as peaks and valleys.These stages are fundamental to conventional THz 3D CT imaging.While they can yield decent reconstructed THz 3D images, they are still limited by computational cost and complexity, the need for extensive data processing, and potential issues related to artifacts in the reconstructed image.Nevertheless, they continue to serve as the cornerstone methods in THz 3D CT imaging, thanks to their firmly established theoretical foundation and widespread implementation.

Variational Frameworks and Numerical Algorithm with ADMM
In this section, we present a variational approach based on the adapted Euler-Elastica-based formulation and implement it through a numerical algorithm using the alternating direction method of multipliers (ADMM).The variational approach, often referred to as the mathematical approach, and the physical modeling approach are two distinct methodologies used in the field of computational modeling.The physical modeling approach is domain-specific and directly relies on the physical laws governing a particular system, whereas variational frameworks are primarily rooted in mathematical principles and optimization techniques.The variational frameworks operate based on the principle of optimizing a function, often referred to as a functional, which typically represents energy/cost forms associated with a physical or computational system.Through the optimization of this functional, we can infer the optimal parameters or states of the system.This approach is more versatile and applicable to a broad spectrum of problems.
Our proposed variational framework is particularly designed to address the issues of prolonged data acquisition time and resolution constraints by scanning resolution with diffraction effects (referred to as concerns (i) and (ii) in Section 1).This framework has the ability to handle LR input slices, even if they are limited in number, and adeptly reconstruct them into a smooth SR surface.Reconstruction from these  slices {Π  } =1,..., can be framed as a variational model, which is motivated by the works [15,18,[56][57][58].Mathematically, the goal is to find the (local) minimum  * that satisfies our new adapted Euler-Elastica-based formulation arg min Here,  ∈ Ω ⊂ R 3 is the phase-field function to provide a 3D representation of the object, governed by the profile function ;  is the profile function designed to be piecewise, taking a value of 1 within the object, half on the boundary, and 0 outside the object, to maintain continuity and smoothness in the phase-field function; ∇, △ is the gradient and Laplacian operator (note that, △ = div ∇ with divergence operator div);  signifies the phase-field variable and is linked to the thickness parameter;   ≤  ≤   is the linear obstacle restrictions, ensuring that the shape of the input data is preserved where the phase-field profiles,   and   , are derived from the interior {Π   } and exterior region {Π   } through the profile function ; and  () is the double-well potential defined by  () = 1 2  2 (1 − ) 2 , which has two minima related to the profile function  and its first derivative is  ′ () = ( − 1) (2 − 1).Specifically, each pair of coefficients (m, n) is used to control the type of formulation applied.If m = 1, n = 0, the formulation is perimeter-based P  from the perimeter energy, focusing on preserving the edge features of the object.Conversely, if m = 0, n = 1, the Willmore-based formulation W  is applied, which is related to the mean curvature from Willmore energy and promotes smoothness in the output.Remark that, without certain restrictions, this formulation could lead to an output that takes the shape of a sphere, as the sphere is the configuration that minimizes the Willmore energy.Lastly, when m = 1, n = 1, the Euler-Elastica-based formulation E  is used, which is an optimal combination of edge preservation and smoothness.For detailed mathematical preliminaries and clarifications, please refer to Sec. 1 of Supplement 1 [15,18,57,58].
To address our model as represented in (1), we employ the alternating direction method of multiplier (ADMM).This approach breaks down complex optimization problems into smaller, more manageable parts for swifter numerical approximation.We further implement the linear obstacle restriction through orthogonal projection, thereby managing the inequality with improved flexibility, which is achieved by With w = ∇, the augmented Lagrangian functional is then expressed using a penalty parameter  > 0 and a Lagrange multiplier , as follows Hence, the application of the ADMM involves a three-step process: solving two subproblems , w and updating one multiplier Utilizing the numerical Euler semi-implicit discretization scheme with a time step , the optimal solutions for the -and w-subproblems in (4) can be updated as follows: Lastly, recall that the multiplier  will be updated by  +1 =   +  (∇ +1 − w +1 ) .This marks the completion of one iteration of the ADMM method where the comprehensive algorithm is presented in Algorithm 1 (more details of derivations in Sec. 2 of Supplement 1).

Experimental Results with Discussion and Application
In this section, we showcase experimental results obtained using our BLIss.To facilitate our experiments, these objects are fabricated economically using a fused deposition modeling (FDM) 3D printer and high impact polystyrene (HIPS) material (see optical images in Sec. 3 of Supplement 1).The choice of HIPS material is for its low absorption in the THz spectrum [59].This selection helps prevent considerable degradation in retrieved time-resolved THz signals.
Building on the conventional techniques for THz 3D CT imaging outlined in Section 2.2, we generate preliminary LR reconstructions by stacking the designated  slices for three objects, as illustrated in Fig. 4   In our implementation, we preserve the dimensional consistency between inputs and outputs for each object, maintaining isotropic resolution even when using fewer lateral samples.Specifically, our algorithm records and processes all data within the LR 3D matrix space R 3   ×   ×   , ensuring a consistent dimensional framework throughout.All results are then reconstructed using isosurface, a 3D triangular mesh representation well-suited for our framework.This method allows for the precise reconstruction of details, effectively demonstrating the isotropic resolution capabilities of our SR algorithm.However, due to the inconsistency of triangular mesh surfaces reconstructed by two formulations W  and E  in our BLIss (i.e.variations in the number of vertices and faces), a fair quantitative comparison is challenging.To address this, we introduce a metric to evaluate the performance of surfaces reconstructed by both formulations, W  and E  , serving as the indicator of the global smoothness of the surfaces.This metric is to calculate the standard deviation of Gaussian curvatures (GC)  GC and of mean curvatures (MC)  MC based on the theory from discrete geometry (details in Sec.4A of Supplement 1).In Table 1, we provide the values of  GC and  MC for various input slices and gap(s).For identical input slices, both   GC and  MC metrics indicate that results obtained using E  exhibit superior smoothness when compared to W  .This is attributed to the ability of E  to automatically handle topologically complex geometries, which overcomes the limitations of W  .Note that a lower value in these curvature computations signifies better smoothness, indicating less variability in the curvature.Besides, as the number of input slices decreases for the same object, both  GC and  MC values also decrease using the same formulation.This trend implies that a reduced number of input slices leads to the loss of finer details in the reconstruction.Even though our framework visually demonstrates its capacity for reconstructed results using limited input slices, we still require accuracy verification to confirm the fidelity compared to the reconstructed results using full input slices by the same formulation.To achieve this accuracy verification, we employ the multi-scale structural similarity index measure (MS-SSIM) to ascertain the precision of reconstructions from sparse slices (details in Sec.4B of Supplement 1).Then, the result obtained with full input slices serves as the benchmark for these evaluations, which is indicated as one in Table 1.When comparing results using the same formulation with different input slices, the MS-SSIM values consistently fall within the range of 0.97 to 0.99.These high MS-SSIM scores underscore the 3D structural fidelity of our reconstructions, even when using only 1/6 of the original slices in our BLIss.Additionally, we present the average elapsed time for each iteration during the reconstruction of the three objects in Table 1.Here, our framework achieves an average elapsed time per iteration within the scale of seconds to show the efficiency of our proposed BLIss.In expanding the discussion on our findings, we emphasize the real-world applications and contributions of our research to the field of THz imaging.Our results pave the way for advanced diagnostic and analytical techniques in non-invasive medical imaging, material characterization, and security screening [4-7, 12, 13].By enhancing the resolution and fidelity of THz imaging, our approach could significantly improve the detection and analysis of bioinformatics, the multi-scale information inside materials, and the inspection of art and historical artifacts without causing damage.These applications underscore the potential of our findings to revolutionize various industries by providing safer, more accurate, and more detailed imaging capabilities.To illustrate the practical application of this research, we have chosen a peanut as our test subject.The peanut exhibits a textured surface and contains two inner seeds along with an air gap.Refer to the optical image of the halved peanut parts after data acquisition in Fig. 7(a).We first present the results in Fig. 7(b)-7(c), using the Time-MAX information to reconstruct the outer shell.Furthermore, our approach can be adapted to other data types obtained from the time-resolved THz pulse, which contains hyperspectral information in a broad frequency range.One notable aspect is the phase in the frequency domain, which offers distinct advantages, such as depth information, material properties, and better image contrast.Therefore, we derive the phase information in the frequency domain by Φ() = tan −1 (/), from the Fourier amplitude ) −2   d of the time-resolved signals  () by Fourier Transform (FT) F. A representative projection of the unwrapped phase at a chosen frequency of 1.0010 THz is illustrated in Fig. 7(d).The selected frequency is determined by the aim of enhancing material contrast, focusing on fat content for the peanut kernels and fiber for the outer shell, within the THz frequency range.When applying our framework with only 1/3 of the input slices, reducing 3-fold data acquisition time, as shown in Fig. 7(e), we obtain the SR results in Fig. 7(f) using W  and Fig. 7(g) using E  .Notably, results produced by E  exhibit superior smoothness, as indicated by lower values of  GC and  MC .Additionally, the reconstruction process can be completed in seconds using a commonplace laptop, reaffirming the precision and efficiency of our framework with very limited dataset inputs.

Conclusion
In our pursuit of advancing pulsed THz 3D SR reconstruction, we propose a promising variational framework based on the W  and an adapted Euler-Elastica-based formulation E  .With our 3D representation utilizing the implicit phase-field function in conjunction with the designed variational framework, we achieve super-resolution results, surpassing the resolution limitations imposed by coarse input slices due to scanning resolution (see Sec. 5 of Supplement 1).Furthermore, this framework is versatile, capable of processing a reduced number of LR input slices, speeding up both data acquisition and computational processes.The adaptability to a range of data types boosts its value, particularly in the context of pulsed THz imaging.To substantiate our results and provide a quantitative assessment of their quality, we employ two metrics: global smoothness, measured by the standard deviation of Gaussian curvatures  GC and of mean curvatures  MC , and accuracy verification using MS-SSIM.The capabilities of this framework are illustrated through its application to non-destructive THz imaging of peanuts, where it significantly improved the precision and clarity of capturing the internal structure.Our work makes substantial strides in addressing long-standing challenges within THz 3D imaging, which have been intricately linked to wave propagation physics and hardware limitations.This advancement has the potential to transform the landscape of THz imaging, expanding its utility across both research and industry.Moreover, its implications extend to other domains, including X-ray and MRI imaging.By unlocking the capability to achieve higher resolution with sparse inputs, we are paving the way to unveil previously hidden details, offering profound insights into materials and processes that were once enigmatic.
While our findings lay a foundation for improved THz 3D reconstruction, the voyage of discovery is far from complete.One promising avenue involves adapting our framework to a broader range of materials and applications, extending its universal applicability.Beyond this, integrating deep learning approaches holds untapped potential to further optimize the imaging process, making it more adaptive and intuitive.Additionally, as hardware technology advances, it could be seamlessly integrated with our approach, further elevating resolution and image clarity (see Sec. 6 of Supplement 1).In the ever-evolving realm of THz imaging, consider our work as the cornerstone, setting the stage for future layers of transformative innovation.

Fig. 1 .
Fig. 1.Schematic illustration (top) and visualization (bottom) of BLIss in our THz SR 3D reconstruction: conventional methodologies are integrated with proposed variational models to improve quality and resolution (highlighted nodes within the inner color).

Fig. 2 .
Fig. 2. (a) Schematic diagram of the THz 3D imaging system setup using transmission geometry (TIA: transimpedance amplifier).(b) Illustrations of one THz pulse without averaging in the time domain.(c) Illustrations of normalized THz power spectrum obtained through fast Fourier transform (FFT) operation from (b).

Fig. 4 .
Fig. 4. Illustrations of Deer object: (a) conventional LR reconstruction with 218 slices where smooth SR reconstruction (e) by Willmore-based W  and (i) by Euler-Elasticabased model E  ; (b) input 109 slices with 1 gap where reconstruction (f) by W  and (j) by E  ; (c) input 55 slices with 3 gaps where reconstruction (g) by W  and (k) by E  ; (d) input 37 slices with 5 gaps where reconstruction (h) by W  and (l) by E  .

Fig. 5 .
Fig. 5. Illustrations of Polarbear object: (a) conventional LR reconstruction with 258 slices where smooth reconstruction (e) by Willmore-based W  and (i) by Euler-Elasticabased formulation E  ; (b) input 129 slices with 1 gap where reconstruction (f) by W  and (j) by E  ; (c) input 65 slices with 3 gaps where reconstruction (g) by W  and (k) by E  ; (d) input 43 slices with 5 gaps where reconstruction (h) by W  and (l) by E  .

Fig. 6 .
Fig. 6.Illustrations of Skull object: (a) conventional LR reconstruction with 138 slices where smooth SR reconstruction (e) by Willmore-based W  and (i) by Euler-Elasticabased formulation E  ; (b) input 69 slices with 1 gap where reconstruction (f) by W  and (j) by E  ; (c) input 35 slices with 3 gaps where reconstruction (g) by W  and (k) by E  ; (d) input 23 slices with 5 gaps where reconstruction (h) by W  and (l) by E  .

Fig. 7 .
Fig. 7. Illustrations of (a) the optical image of the unveiled peanut kernels.(b) projected Time-MAX THz 2D image.(c) reconstructed THz 3D image of peanuts from (b).(d) representative projection of unwrapped phase at 1.0010 THz.(e) 46 slices extracted from (d) by IRT.(f) reconstruction of peanut kernels by Willmore model W  from (e) with higher computational time (more iterations) and lower smoothness.(g) optimized reconstruction of kernels by Euler-Elastica model E  from (e) with lower computational time (less iterations) and higher smoothness.

Table 1 .
Comparison metrics for three reconstructed objects Deer, Polarbear, and Skull with various input slices and gap(s) by Willmore-based W  and Euler-Elastica-based E  formulations: the standard deviation of Gaussian curvatures  GC and of mean curvatures  MC ; the average elapsed time of each iteration (seconds/iteration or s/iter); and the multi-scale structural similarity index measure (MS-SSIM) where ↑ (↓): higher (lower) is better.