Fast Terahertz Imaging Model Based on Group Sparsity and Nonlocal Self-Similarity

Abstract

In order to solve the problems of long-term image acquisition time and massive data processing in a terahertz time domain spectroscopy imaging system, a novel fast terahertz imaging model, combined with group sparsity and nonlocal self-similarity (GSNS), is proposed in this paper. In GSNS, the structure similarity and sparsity of image patches in both two-dimensional and three-dimensional space are utilized to obtain high-quality terahertz images. It has the advantages of detail clarity and edge preservation. Furthermore, to overcome the high computational costs of matrix inversion in traditional split Bregman iteration, an acceleration scheme based on conjugate gradient method is proposed to solve the terahertz imaging model more efficiently. Experiments results demonstrate that the proposed approach can lead to better terahertz image reconstruction performance at low sampling rates.

1. Introduction

The information obtained from imaging at the Terahertz (THz) frequencies (0.1–10 THz) has received much attention in recent decades due to its unique properties [1,2,3,4,5,6]. Terahertz time domain spectroscopy (THZ-TDS) imaging is one of the most powerful techniques in the terahertz imaging field. It can reconstruct the terahertz image with high signal-to-noise (SNR) ratio and high resolution by scanning the target pixel by pixel. Therefore, it has great application prospects in fields requiring high terahertz imaging quality [7,8,9]. However, THz-TDS imaging technology requires a long time for terahertz image acquisition. Therefore, a fast imaging scheme that can solve the disadvantage of long-time acquisition while maintaining the advantages of THz-TDS detection will attract great interest for practical application.

Numerous techniques have been proposed to solve the problem of long-time acquisition of THz-TDS technology. Terahertz detector arrays have shown great application potential in shortening long-time acquisition [10,11,12], but they suffer from high equipment and complexity costs. A Fourier based terahertz imaging method was proposed to achieve fast terahertz image reconstruction by using sparsity prior to the image [13]. However, it still needs raster scanning on the Fourier plane. Then, a promising terahertz imaging paradigm, namely single pixel imaging system, was proposed [14,15]. This system uses only one detector and avoids raster scanning of the object. Despite the fact that this system could decrease the number of measurements, the speed of switching from one spatial pattern to another is time consuming. Then, various techniques for effective spatial modulation of terahertz beams were studied, but still required additional hardware. Moreover, a block compressive sensing (CS) based method was proposed for terahertz imaging in the spatial domain [16,17], which requires no additional hardware devices and can shorten imaging time. However, a single sparsity prior is utilized in this method for terahertz imaging, and the reconstructed terahertz image quality is not satisfactory. A dual sparsity constraint based spatial domain terahertz imaging strategy was proposed in [18], but may suffer from over-smoothing when the sampling rate is reduced.

In the past few years, a series of nonlocal self-similarity methods have been proposed to extend image denoising from two-dimensional space to three-dimensional space [19,20,21]. These methods use the nonlocal self-similarity of image patches to improve the image quality. A hybrid sparsity model was proposed for terahertz imaging in [22], which utilized the local sparsity and nonlocal self-similarity of the terahertz image to improve image quality, but it does not use the image structure information. Recently, in order to remove noise and reconstruct the image more effectively, methods based on group sparsity have been widely used in image processing, which can extract sparse information from the image structure [23,24,25]. Inspired by the success of group sparsity and nonlocal self-similarity in image processing, a novel fast terahertz imaging model, combined with group sparsity and nonlocal self-similarity (GSNS) is proposed in this paper. It exploits the structure sparsity and similarity of terahertz image patches in both two-dimensional and three-dimensional space to ensure high-quality terahertz image reconstruction, and has the advantages of detail clarity and edge preservation. Furthermore, to overcome the high computational costs of the proposed terahertz imaging model, an acceleration scheme based on a conjugate gradient method is proposed to solve the model more efficiently. Compared with the other methods, the proposed GSNS algorithm can achieve superior performance in fast terahertz imaging.

The rest of the paper is structured as follows. In Section 2, we present the proposed fast terahertz imaging model based on group sparsity and nonlocal self-similarity, and give a solution algorithm using an acceleration scheme in detail. The experiments are utilized to demonstrate the performance of the proposed method in Section 3. Finally, the conclusion is drawn in Section 4.

2. Materials and Methods

The THZ-TDS system realizes image reconstruction by scanning each pixel of the target sample. The image acquisition time of THZ-TDS is determined by the number of scanned pixels. Therefore, in order to achieve fast terahertz image reconstruction, we need to reduce the number of scanned pixels. Given a target image $x\in R^N$, we randomly select M positions from a total of N pixels of the target image for terahertz scanning detection. Then the fast terahertz image reconstruction model of the THz-TDS system can be represented as

(1)$\mathit{y}=\mathbf{R}\mathit{x}$

where $\mathit{y}\in R^M$ denotes the measured vector of terahertz data, $\mathit{x}\in R^N$ is the target image to be reconstructed, and $\mathbf{R}$ denotes the observation matrix composed of 0 and 1 values. The positions of elements with value 1 in the matrix $\mathbf{R}$ are determined by the scanning positions of the terahertz system, and the values of other elements are 0.

The purpose of terahertz image reconstruction is to obtain the target image x from the sparse measured vector y. Since the amount of measured terahertz data is less than the number of pixels in the target image, problem (1) appears ill conditioned. Since the spectral density of an common terahertz image is usually concentrated in the low frequency band and has strong sparsity, imaging methods based on CS can be utilized to reconstruct the terahertz image to reduce the image acquisition time. Based on the sparsity of the terahertz image in frequency domain, the fast terahertz image reconstruction model can be written as

(2)$\begin{array}{l}\underset{x}{\min }\left\{{‖F\mathit{x}‖}_1\right\}\hspace{1em}\\ s.t.\hspace{1em}\mathit{y}=\mathbf{R}\mathit{x}\end{array}$

where F is the sparse transformation.

2.1. Proposed Fast Terahertz Imaging Model

To further improve the performance of terahertz image reconstruction, a novel fast terahertz imaging model which integrates group sparsity and nonlocal self-similarity constraints is proposed in this section.

Given an image x, we divide it into P patches of size L × L, and the p_th patch can be denoted by

(3)${\mathit{x}}_p={\mathit{C}}_p\mathit{x}$

where C_p represents the operation of dividing the image into small patches. For each patch x_p, we search for its c-1 best matching patches in a T × T search window, and then stack all the best matching patches to form a two-dimensional matrix g_p of size L² × c, where each similar patch is used as a column of g_p. Then the matrix g_p can be expressed as

(4)${\mathit{g}}_p={\mathit{G}}_p{\mathit{x}}_p={\mathit{G}}_p{\mathit{C}}_p\mathit{x}$

where ${\mathit{G}}_p$ is an operator that extracts all patches with similar structures of x_p to construct a patch group. For a given dictionary D_p learned from the group, the patch group g_p can be sparsely expressed as

(5)$\begin{array}{l}{\mathit{\alpha }}_p=\arg \underset{{\mathit{\alpha }}_p}{\min }{‖{\mathit{\alpha }}_p‖}_1\\ s.t.{\mathit{g}}_p={\mathit{D}}_p{\mathit{\alpha }}_p\end{array}$

where ${\mathit{\alpha }}_p$ denotes the sparse coefficient vector. Using the structure similarity and sparsity of the terahertz image, group sparsity could obtain better reconstruction results. After obtaining all the sparse coefficient vectors ${\mathit{\alpha }}_p,p=1,2,\dots P$, the entire image x can be obtained by [25]

(6)$\mathit{x}={\displaystyle \sum _{p=1}^P{\mathit{G}}_p^T}\left({\mathit{D}}_p{\mathit{\alpha }}_p\right)·/{\displaystyle \sum _{p=1}^P{\mathit{G}}_p^T{\mathbf{1}}_{L^2\times c}}$

where ${\mathbf{1}}_{L^2\times c}$ is a $L^2\times c$ matrix and all its elements are 1. ${\mathit{G}}_p^T$ is the transpose of G_p, which can put a group back to the p-th position of the reconstructed image and fill in zeros at other positions. Equation (6) shows that we can recover the image x by averaging all the overlapped groups.

In addition, motivated by the successful application of nonlocal self-similarity in the three-dimensional transform domain for image denoising, we utilize the nonlocal self-similarity of the terahertz image to explore the structure similarity and sparsity in three-dimensional space in order to ensure high-quality terahertz image reconstruction. Similar to the group sparsity, for a patch x_p, we search for its c-1 best matching patches in a T × T search window. However, the difference from group sparsity is that all the best matching patches are stacked to form a three-dimensional matrix z_p with the size L × L$\times$c. Then the matrix z_p can be expressed as

(7)${\mathit{z}}_p={\mathit{S}}_p{\mathit{x}}_p={\mathit{S}}_p{\mathit{C}}_p\mathit{x}$

where ${\mathit{S}}_p$ is an operator that extracts all patches with similar structures of x_p to construct a three-dimensional matrix z_p. Let T_3D be an orthogonal three-dimensional transformation operator, then the transformation coefficients of the image can be expressed as

(8)${\mathbf{\Theta }}_p={\mathit{T}}_{3D}{\mathit{z}}_p={\mathit{T}}_{3D}{\mathit{S}}_p{\mathit{C}}_p\mathit{x}={\mathrm{\Psi }}_p^{3D}\mathit{x}$

where ${\mathrm{\Psi }}_p^{3D}={\mathit{T}}_{3D}{S}_p{\mathit{C}}_p$ is the operator of nonlocal self-similarity.

Therefore, nonlocal self-similarity can be used to measure the structure similarity and sparsity of terahertz images in three-dimensional space, which can be expressed as

(9)${\displaystyle \sum _{p=1}^P{‖{\mathbf{\Theta }}_p‖}_1}={\displaystyle \sum _{p=1}^P{‖{\mathrm{\Psi }}_p^{3D}\mathit{x}‖}_1}$

In summary, combining group sparsity and nonlocal self-similarity, we propose a new fast terahertz imaging model, which can be expressed as

(10)$\underset{x}{\min }{\displaystyle \sum _{p=1}^P\left({‖{\mathit{\alpha }}_p‖}_1+{‖{\mathrm{\Psi }}_p^{3D}\mathit{x}‖}_1\right)+}\frac{\mu }{2}{‖\mathbf{R}\mathit{x}-\mathit{y}‖}_2^2$

where μ is the positive regularization parameter. Obviously, it is extremely difficult to solve the optimization problem (10) directly. Therefore, solving this problem (10) effectively is one of the main contributions of this paper.

2.2. Solvution Algorithm Using Acceleration Scheme

Split Bregman iteration can transform a difficult optimization problem into several simple subproblems. Then the subproblems are solved and updated alternately [26,27,28]. Based on the methodology of split Bregman iteration, we introduce ${\mathrm{\Theta }}_p={\mathrm{\Psi }}_p^{3D}\mathit{x}$, ${\mathit{\alpha }}_p={\mathit{D}}_p^{-1}{\mathit{g}}_p={\mathit{D}}_p^{-1}{\mathit{G}}_p{\mathit{C}}_p\mathit{x}={\mathrm{\Psi }}_p^{2D}\mathit{x}$, Bregman variables ${\mathit{b}}_{\alpha p}^{}$ and ${\mathit{b}}_{\mathrm{\Theta }p}^{}$, convert the problem (10) into the equivalent split Bregman formula, then the problem (10) can be rewritten as

(11)$\underset{x,r_e,{\mathrm{\Theta }}_p}{\min }{\displaystyle \sum _{p=1}^P{‖{\mathit{\alpha }}_p‖}_1}+{\displaystyle \sum _{p=1}^P{‖{\mathbf{\Theta }}_p‖}_1}+\frac{\mu }{2}{‖\mathbf{R}x-y‖}_2^2+\frac{\lambda }{2}{\displaystyle \sum _{p=1}^P{‖{\mathit{\alpha }}_p-{\mathrm{\Psi }}_p^{2D}\mathit{x}-{\mathit{b}}_{\alpha p}‖}_2^2}+\frac{\gamma }{2}{\displaystyle \sum _{p=1}^P{‖{\mathbf{\Theta }}_p-{\mathrm{\Psi }}_p^{3D}\mathit{x}-{\mathit{b}}_{\mathrm{\Theta }p}‖}_2^2}\hspace{1em}$

where λ and γ are the positive regularization parameters. Then the split Bregman iteration formula of the problem (10) is

(12)$\begin{array}{l}({\mathit{x}}^{i+1},{\mathit{\alpha }}_p^{i+1},{\mathbf{\Theta }}_p^{i+1})=\arg \underset{x,{\mathit{\alpha }}_p,{\mathbf{\Theta }}_p}{\min }{\displaystyle \sum _{p=1}^P{‖{\mathit{\alpha }}_p‖}_1}+{\displaystyle \sum _{p=1}^P{‖{\mathbf{\Theta }}_p‖}_1}+\frac{\mu }{2}{‖\mathbf{R}\mathit{x}-\mathit{y}‖}_2^2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{\lambda }{2}{\displaystyle \sum _{p=1}^P{‖{\mathit{\alpha }}_p-{\mathrm{\Psi }}_p^{2D}\mathit{x}-{\mathit{b}}_{\alpha p}^i‖}_2^2}+\frac{\gamma }{2}{\displaystyle \sum _{p=1}^P{‖{\mathbf{\Theta }}_p-{\mathrm{\Psi }}_p^{3D}\mathit{x}-{\mathit{b}}_{\mathrm{\Theta }p}^i‖}_2^2}\end{array}$

and the Bregman updates are

(13)$\left\{\begin{array}{c}{\mathit{b}}_{\alpha p}^{i+1}={\mathit{b}}_{\alpha p}^i+({\mathrm{\Psi }}_p^{2D}{\mathit{x}}^{i+1}-{\mathit{\alpha }}_p^{i+1})\\ {\mathit{b}}_{\mathrm{\Theta }p}^{i+1}={\mathit{b}}_{\mathrm{\Theta }p}^i+({\mathrm{\Psi }}_p^{3D}{\mathit{x}}^{i+1}-{\mathbf{\Theta }}_p^{i+1})\end{array}\right.$

The l₁ norm in the optimization problem (12) is not differentiable, so it is still difficult to solve directly. Based on the strategy of separating variables, we convert the optimization problem (12) into three subproblems, and iteratively update each variable by fixing other variables.

x subproblem: fixing ${\mathit{\alpha }}_p$,${\mathbf{\Theta }}_p$,${\mathit{b}}_{\alpha p}^{}$ and ${\mathit{b}}_{\mathrm{\Theta }p}^{}$, the subproblem of updating x can be given from (12)

(14)$\begin{array}{l}{\mathit{x}}^{i+1}=\underset{x}{\min }{Q}_1(\mathit{x})\\ =\underset{x}{\min }\frac{\mu }{2}{‖\mathbf{R}\mathit{x}-\mathit{y}‖}_2^2+\frac{\lambda }{2}{\displaystyle \sum _{p=1}^P{‖{\mathit{\alpha }}_p^i-{\mathrm{\Psi }}_p^{2D}\mathit{x}-{\mathit{b}}_{\alpha p}^i‖}_2^2}+\frac{\gamma }{2}{\displaystyle \sum _{p=1}^P{‖{\mathbf{\Theta }}_p^i-{\mathrm{\Psi }}_p^{3D}\mathit{x}-{\mathit{b}}_{\mathrm{\Theta }p}^i‖}_2^2}\hspace{1em}\end{array}$

The subproblem (14) is a minimization problem of typical convex function. Using the differential for $Q_1(x)$ and setting the result equal to zero, we can obtain

(15)$⟮\mu {\mathbf{R}}^{\mathrm{T}}\mathbf{R}+\lambda {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{2D})}^{\mathrm{T}}{\mathrm{\Psi }}_p^{2D}}+\gamma {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{3D})}^{\mathrm{T}}{\mathrm{\Psi }}_p^{3D}}⟯{\mathit{x}}^{i+1}=\mu {\mathbf{R}}^{\mathrm{T}}\mathit{y}+\lambda {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{2D})}^{\mathrm{T}}\left({\mathit{\alpha }}_p^i-{\mathit{b}}_{\alpha p}^i\right)}+\gamma {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{3D})}^{\mathrm{T}}\left({\mathbf{\Theta }}_p^i-{\mathit{b}}_{\mathrm{\Theta }p}^i\right)}$

Then, x can be updated as

(16)${\mathit{x}}^{i+1}={⟮\mu {\mathbf{R}}^{\mathrm{T}}\mathbf{R}+\lambda {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{2D})}^{\mathrm{T}}{\mathrm{\Psi }}_p^{2D}}+\gamma {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{3D})}^{\mathrm{T}}{\mathrm{\Psi }}_p^{3D}}⟯}^{-1}{\mathit{z}}^i$

where

(17)${\mathit{z}}^i=\mu {\mathbf{R}}^{\mathrm{T}}\mathit{y}+\lambda {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{2D})}^{\mathrm{T}}\left({\mathit{\alpha }}_p^i-{\mathit{b}}_{\alpha p}^i\right)}+\gamma {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{3D})}^{\mathrm{T}}\left({\mathrm{\Theta }}_p^i-{\mathit{b}}_{\mathrm{\Theta }p}^i\right)}$

${\mathbf{\Theta }}_p$ subproblem: fixing $\mathit{x}$,${\mathit{\alpha }}_p$,${\mathit{b}}_{\alpha p}^{}$ and ${\mathit{b}}_{\mathrm{\Theta }p}^{}$, the subproblem of updating ${\mathbf{\Theta }}_p$ can be given from (12)

(18)${\mathbf{\Theta }}_p^{i+1}=\underset{{\mathbf{\Theta }}_p}{\min }{\displaystyle \sum _{p=1}^P\left({‖{\mathrm{\Theta }}_p‖}_1+\frac{\gamma }{2}{‖{\mathbf{\Theta }}_p-{\mathrm{\Psi }}_p^{3D}{\mathit{x}}^{i+1}-{\mathit{b}}_{\mathrm{\Theta }p}^i‖}_2^2\right)}$

By introducing the soft thresholding algorithm, each ${\mathbf{\Theta }}_p$ can be updated as [29]

(19)${\mathbf{\Theta }}_p^{i+1}=\mathrm{shrink}\left({\mathrm{\Psi }}_p^{3D}{\mathit{x}}^{i+1}+{\mathit{b}}_{\mathrm{\Theta }p}^i,1/\gamma \right)$

where

(20)$\begin{array}{l}\mathrm{shrink}\left(x,\lambda \right)=\mathrm{sgn}(x)\times \max \left(\left|x\right|-\lambda ,0\right)\\ =\left\{\begin{array}{c}x-\lambda ,\\ 0,\\ x+\lambda ,\end{array}\right.\hspace{1em}\begin{array}{c}x\in \left(\lambda ,+\infty \right)\\ x\in \left(-\lambda ,\lambda \right)\\ \hspace{1em}x\in \left(-\infty ,-\lambda \right)\end{array}\end{array}$

${\mathit{\alpha }}_p$ subproblem: Fixing $\mathit{x}$,${\mathbf{\Theta }}_p$,${\mathit{b}}_{\alpha p}^{}$ and ${\mathit{b}}_{\mathrm{\Theta }p}^{}$, the subproblem of updating ${\mathit{\alpha }}_p$ can be given from (12)

(21)${\mathit{\alpha }}_p^{i+1}=\underset{{\mathit{\alpha }}_p}{\min }{\displaystyle \sum _{p=1}^P⟮{‖{\mathit{\alpha }}_p‖}_1+\frac{\lambda }{2}{‖{\mathit{\alpha }}_p-{\mathrm{\Psi }}_p^{2D}\mathit{x}-{\mathit{b}}_{\alpha p}^i‖}_2^2⟯}$

Similar to the optimization problem (18), each ${\mathit{\alpha }}_p$ can be updated by the soft thresholding algorithm.

(22)${\mathit{\alpha }}_p^{i+1}=\mathrm{shrink}\left({\mathrm{\Psi }}_p^{2D}{\mathit{x}}^{i+1}+{\mathit{b}}_{\alpha p}^i,1/\lambda \right)$

Then, these simple subproblems are solved iteratively to obtain the reconstructed terahertz image.

Although the optimization problem (12) can be solved by iterating three subproblems (16), (18) and (21), the computational complexity is still very high. The main computational complexity comes from solving the x subproblem because the cost of computing the matrix inversion is too high, which is O(N³). So, an effective method is highly desirable. Conjugate gradient method is an effective strategy for solving large-scale optimization problems, which has a simple iterative form [30,31]. The basic concept of the conjugate gradient method is to combine the conjugate property with the steepest descent method, construct a set of conjugate directions by using the gradient at known points, search along this set of directions, and find the minimum point of the objective function. Therefore, an acceleration scheme based on the conjugate gradient method is proposed to tackle the x subproblem.

For convenience, we let

(23)$\mathit{A}=\mu {\mathbf{R}}^{\mathrm{T}}\mathbf{R}+\lambda {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{2D})}^{\mathrm{T}}{\mathrm{\Psi }}_p^{2D}}+\gamma {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{3D})}^{\mathrm{T}}{\mathrm{\Psi }}_p^{3D}}$

(24)$\mathit{z}=\mu {\mathbf{R}}^{\mathrm{T}}\mathit{y}+\lambda {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{2D})}^{\mathrm{T}}\left({\mathbf{\alpha }}_p^i-{\mathit{b}}_{\alpha p}^i\right)}+\gamma {\displaystyle \sum _{p=1}^P{({\mathrm{\Psi }}_p^{3D})}^{\mathrm{T}}\left({\mathbf{\Theta }}_p^i-{\mathit{b}}_{\mathrm{\Theta }p}^i\right)}$

(25)$\mathit{u}={\mathit{x}}^{i+1}$

Then, (15) can be rewritten as

(26)$\mathit{Au}=\mathit{z}$

Therefore, the solution of the Equation (26) can be equivalent to

(27)$\underset{u}{\min }f(\mathit{u})=\underset{u}{\min }\frac{1}{2}{\mathit{u}}^{\mathrm{T}}\mathit{Au}-{\mathit{z}}^T\mathit{u}\hspace{1em}$

It can be seen that the optimization problem (27) is a quadratic programming problem, which can be solved by the conjugate gradient method.

Firstly, given an initial point u⁰, calculate the gradient of the objective function f (u) at this point u⁰ and set the first search direction d⁰ as the initial negative gradient direction

(28)${\mathit{d}}^0=-\nabla f({\mathit{u}}^0)=-{\mathit{g}}^0$

Search along direction d⁰ to obtain the next iteration point u¹. Calculate the gradient g¹ of f (u) at u¹. If ||g¹|| ≠ 0, use g¹ and d⁰ to construct the second search direction d¹, and search along d¹ to obtain point u². Generally, if the point u^k and the search direction d^k are known, we search along d^k starting from u^k and get the next iteration point

(29)${\mathit{u}}^{k+1}={\mathit{u}}^k+{\alpha }^k{\mathit{d}}^k\hspace{1em}$

where ${\alpha }^k$ is step size, satisfying

(30)${\alpha }^k=\arg \underset{{\alpha }^k}{\min }f({\mathit{u}}^k+{\alpha }^k{\mathit{d}}^k)$

Solve (30) to obtain

(31)${\alpha }^k=\frac{{({\mathit{g}}^k)}^T{\mathit{d}}^k}{{({\mathit{d}}^k)}^T\mathit{A}{\mathit{d}}^k}\hspace{1em}$

and the gradient g^{k + 1} of the objective function f (u) at the current iteration point u^{k + 1} is

(32)${\mathit{g}}^{k+1}=\nabla f({\mathit{u}}^{k+1})={\mathit{g}}^k-{\alpha }^k\mathit{A}{\mathit{d}}^k\hspace{1em}$

Then, a new direction d^{k + 1} orthogonal to the previous search direction is constructed by using the linear combination of the previous search direction d^k and the gradient g^{k + 1} of the objective function f (u) at the current iteration point u^{k + 1}

(33)${\mathit{d}}^{k+1}={\mathit{g}}^{k+1}+{\beta }^k{\mathit{d}}^k\hspace{1em}$

According to the conjugation of d^{k + 1} and d^k with respect to A, we get

(34)${\beta }^k=\frac{{({\mathit{g}}^{k+1})}^T{\mathit{g}}^{k+1}}{{({\mathit{g}}^k)}^T{\mathit{g}}^k}\hspace{1em}$

When the iteration ends, we get the solution ${\mathit{x}}^{i+1}$ of the x subproblem. The main computational complexity of the conjugate gradient method comes from the multiplication of matrix and vector, which is O(N²). Therefore, the computational complexity can be reduced by using the conjugate gradient method instead of matrix inversion.

To sum up, we introduce the acceleration scheme into the split Bregman iterative framework, and obtain the main steps of the proposed algorithm summarized in Algorithm 1.

Algorithm 1 The proposed algorithm by combining group sparsity and nonlocal self-similarity.

Input:   measured terahertz data y, observation matrix R, sparsity basis ${\mathrm{\Psi }}_p^{2D}$, nonlocal self-similarity operator ${\mathrm{\Psi }}_p^{3D}$.

Initialization:   ${\mathit{\alpha }}_p^0={\mathbf{\Theta }}_p^0={\mathit{b}}_{\alpha p}^0={\mathit{b}}_{\mathrm{\Theta }p}^{\mathbf{0}}=0$, μ, λ, γ, estimate an initial image x0.

Loop: set i = 0 and repeat until $\left|\right|{\mathit{x}}^{i+1}-{\mathit{x}}^i|{|}_2<\delta$  ${\mathit{g}}^0={\mathit{z}}^i-\mathit{A}{\mathit{x}}^i$  ${\mathit{d}}^0=-{\mathit{g}}^0$  $k=0$  while  ${\alpha }^k=\frac{{({\mathit{g}}^k)}^{\mathrm{T}}{\mathit{d}}^k}{{({\mathit{d}}^k)}^{\mathrm{T}}\mathit{A}{\mathit{d}}^k}$  ${\mathit{u}}^{k+1}={\mathit{u}}^k+{\alpha }^k{\mathit{d}}^k$  ${\mathit{g}}^{k+1}={\mathit{g}}^k-{\alpha }^k\mathit{A}{\mathit{d}}^k$  if $\left|\right|{\mathit{g}}^{k+1}\left|\right|<\epsilon$, break   ${\beta }^k=\frac{{({\mathit{g}}^{k+1})}^T{\mathit{g}}^{k+1}}{{({\mathit{g}}^k)}^{\mathrm{T}}{\mathit{g}}^k}$  ${\mathit{d}}^{k+1}={\mathit{g}}^{k+1}+{\beta }^k{\mathit{d}}^k$  $k=k+1$  return ${\mathit{u}}^{k+1}$  ${\mathit{x}}^{i+1}={\mathit{u}}^{k+1}$  ${\mathrm{\Theta }}_p^{i+1}=\mathrm{shrink}\left({\mathrm{\Psi }}_p^{3D}{x}^{i+1}+{\mathit{b}}_{\mathrm{\Theta }p}^i,1/\gamma \right)$  ${\mathit{\alpha }}_p^{i+1}=\mathrm{shrink}\left({\mathrm{\Psi }}_p^{2D}{\mathit{x}}^{i+1}+{\mathit{b}}_{\alpha p}^i,1/\lambda \right)$  ${\mathit{b}}_{\alpha p}^{i+1}={\mathit{b}}_{\alpha p}^i+({\mathrm{\Psi }}_p^{2D}{\mathit{x}}^{i+1}-{\mathit{\alpha }}_p^{i+1})$  ${\mathit{b}}_{\mathrm{\Theta }p}^{i+1}={\mathit{b}}_{\mathrm{\Theta }p}^i+({\mathrm{\Psi }}_p^{3D}{\mathit{x}}^{i+1}-{\mathrm{\Theta }}_p^{i+1})$  $i=i+1$

End Loop

Output:x

3. Experiments and Discussion

The THz-TDS reflection imaging system with a measurement range of 5 × 5 cm was used in the experiment. The terahertz pulse is generated by the laser with a pulse width of 100 fs and a repetition rate of 80 MHz. The terahertz beam generated by pump beam is focused onto the sample; then the beam reflected from the sample is sent to a ZnTe crystal, where it overlaps with the probe beam. The probe beam is modulated by the terahertz field and then is guided to the detector [18]. During the imaging process, the sample is moved in a mechanical grating scan, and the spatial resolution is set to 0.25 mm. By programming the scan positions determined by the observation matrix, we could easily obtain the proposed fast terahertz imaging system from a traditional THz-TDS system.

In the experiment, a glass fragment, a metal screw, a wood chip and a stone embedded in flour with a depth of 10 mm were measured by THz-TDS system at room temperature. Figure 1a shows the sample placed in a 3.5 cm diameter dish, and Figure 1b presents the terahertz image of the sample obtained from the reflection of four objects and flour through full scan imaging process. The terahertz image of the sample is obtained at 0.8 THz and the number of the full scan pixels is 150 × 150. In the experiment, all parameters are set empirically. The size of the patches is 8 × 8. The horizontal or vertical distance between two adjacent patches is 4. The number of best matching patches is set to 10 and the search window is set to 40 × 40. To quantitatively evaluate the imaging performance of terahertz images, peak signal-to-noise ratio (PSNR) [32] and relative l₂ norm error (RLNE) [33] are used as image quality assessment indexes in this paper, which are

(35)$\mathrm{PSNR}=10{\log }_{10}\frac{peakva{l}^2}{\mathrm{MSE}(\mathit{x},\widehat{\mathit{x}})}$

(36)$\mathrm{RLNE}=\frac{{‖\mathit{x}-\widehat{\mathit{x}}‖}_2}{{‖\mathit{x}‖}_2}$

where $peakval$ denotes the peak value of the terahertz image, and $\mathrm{MSE}(x,\widehat{x})$ denotes the mean square error.

To evaluate the proposed GSNS method, we compare the GSNS with the single sparse constraint algorithm (SSC) [17], dual sparsity constraint algorithm (DSC) [18] and the hybrid sparsity model (HSM) [22]. In the experiment, all the methods were compared using the same observation matrix. The reconstructed terahertz images obtained by the four methods at different sampling rates are shown in Figure 2. Figure 2a,b show the reconstructed terahertz images of SSC when the sampling rate (SR) is 20% and 40%, respectively. Figure 2c,d are the terahertz images recovered by DSC with different sampling rates, Figure 2e,f are the terahertz images obtained by HSM, and the reconstructed terahertz images of the proposed GSNS are displayed in Figure 2g,h. When comparing Figure 2b,d,f,h, we can see that HSM and the proposed method can reconstruct the terahertz images with satisfaction at 40% SR. In particular, the terahertz image obtained by the GSNS has fewer recovery errors. However, at the same SR, the reconstructed results of DSC and SSC present some degradation and blurs. When the SR drops to 20%, it can be seen from the areas marked by circles and squares that the high under-sampling rate leads to a serious partial loss of the SSC method. For a close-up comparison, we enlarged the areas marked by squares to evaluate the image quality in Figure 3. It can be seen that the reconstruction results of DSC and HSM have some fuzzy details. Compared with the other three methods, the proposed GSNS performs well in preserving image details and edges. The terahertz image reconstructed by the proposed method is closer to the full scan terahertz image than the other two methods. Overall, GSNS obtains better results compared with the comparative methods.

To further demonstrate the effectiveness of the proposed method, in Figure 4 we present the PSNR and RLNE curves of SSC, DSC, HSM and GSNS with SR from 5% to 50%. As can be seen from Figure 4, for the terahertz image with different sampling rates, the proposed GSNS has higher PSNR and lower RLNE values than the other three methods. The imaging performance of the proposed method is always better than SSC, DSC and HSM under the same SR.

The above results demonstrated that the proposed algorithm can achieve better reconstruction results. Furthermore, in order to verify the computational efficiency of the proposed algorithm, 40% SR data are used to reconstruct the terahertz image. The computation times before and after improvement are 228.1 s and 94.7 s respectively, which shows the effectiveness of the proposed algorithm.

4. Conclusions

In this paper, we present a fast terahertz imaging model which integrates group sparsity and nonlocal self-similarity constraints for better terahertz image reconstruction. The advantage of the proposed method is that it utilizes the structure similarity and sparsity of image patches in both two-dimensional and three-dimensional space to reconstruct the terahertz image. Moreover, an acceleration scheme based on conjugate gradient method is proposed to overcome the high computational costs of the proposed terahertz imaging model. Compared with existing methods, the experimental results show that the proposed method has better performance in detail clarity and edge preservation of the reconstructed terahertz images.

Author Contributions

Conceptualization, X.R. and Y.N.; methodology, X.R. and Y.N.; software, X.R. and Y.J.; validation, Y.B. and Y.J.; formal analysis, X.R. and Y.J.; investigation, Y.B. and Y.J.; resources, Y.B. and Y.J.; data curation, Y.J.; writing—original draft preparation, X.R. and Y.N.; writing—review and editing, Y.B. and Y.J.; project administration, X.R.; funding acquisition, X.R. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the National Natural Science Foundation of China, grant number U1904120, Henan Science and Technology Department Science and Technology Research Program, grant number 182102310759, the Fundamental Research Funds for the Henan Provincial Colleges and Universities, grant number 2018RCJH09.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. The sample and the terahertz image obtained through full scan imaging process. (a) sample; (b) terahertz image.

Figure 2. Reconstruction images of different methods for the sample at different sampling rates, and the areas marked by circles and squares are selected for detail comparison. (a) SSC at SR = 20%; (b) SSC at SR = 40%; (c) DSC at SR = 20%; (d) DSC at SR = 40%; (e) HSM at SR = 20%; (f) HSM at SR = 40%; (g) GSNS at SR = 20%; (h) GSNS at SR = 40%.

Figure 3. Reconstruction images of different methods for the sample at different sampling rates in the selected areas. (a) SSC at SR = 20%; (b) DSC at SR = 20%; (c) HSM at SR = 20%; (d) GSNS at SR = 20%; (e) SSC at SR = 40%; (f) DSC at SR = 40%; (g) HSM at SR = 40%; (h) GSNS at SR = 40%.

Figure 4. Comparison of the PSNR and RLNE curves to different sampling rates. (a) PSNR; (b) RLNE.

References

1. Mittleman, D.M. Twenty years of terahertz imaging. Opt. Express; 2018; 26, pp. 9417-9431. [DOI: https://dx.doi.org/10.1364/OE.26.009417] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29715894]

2. Zanotto, L.; Piccoli, R.; Dong, J.L.; Morandotti, R.; Razzari, L. Single-pixel terahertz imaging: A review. Opto Electron. Adv.; 2020; 3, 200012. [DOI: https://dx.doi.org/10.29026/oea.2020.200012]

3. Totero Gongora, J.S.; Olivieri, L.; Peters, L.; Tunesi, J.; Cecconi, V.; Cutrona, A.; Tucker, R.; Kumar, V.; Pasquazi, A.; Peccianti, M. Route to intelligent imaging reconstruction via terahertz nonlinear ghost imaging. Micromachines; 2020; 11, 521. [DOI: https://dx.doi.org/10.3390/mi11050521] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32443881]

4. Liu, H.B.; Zhong, H.; Karpowicz, N.; Chen, Y.Q.; Zhang, X.C. Real-time terahertz imaging with a single-pixel detector. Proc. IEEE; 2007; 95, pp. 1514-1527. [DOI: https://dx.doi.org/10.1109/JPROC.2007.898903]

5. Stoik, C.; Bohn, M.; Blackshire, J. Nondestructive evaluation of aircraft composites using reflective terahertz time domain spectroscopy. NDT E Int.; 2020; 43, pp. 106-115. [DOI: https://dx.doi.org/10.1016/j.ndteint.2009.09.005]

6. Guerboukha, H.; Nallappan, K.; Skorobogatiy, M. Toward real-time terahertz imaging. Adv. Opt. Photon.; 2018; 10, pp. 843-938. [DOI: https://dx.doi.org/10.1364/AOP.10.000843]

7. Arco, A.D.; Di Fabrizio, M.; Dolci, V.; Petrarca, M.; Lupi, S. THz pulsed imaging in biomedical applications. Condens. Matter; 2020; 5, 25.

8. Xiang, Y.; Xiang, Z.; Ke, Y.; Liu, Y.; Yang, L. Biomedical applications of terahertz spectroscopy and imaging. Trends Biotechnol.; 2016; 34, pp. 810-824.

9. Pickwell, E.; Wallace, V.P. Biomedical applications of terahertz technology. J. Phys. D Appl. Phys.; 2006; 39, R301. [DOI: https://dx.doi.org/10.1088/0022-3727/39/17/R01]

10. Pradarutti, B.; Müller, R.; Freese, W.; Matthäus, G.; Riehemann, S.; Notni, G.; Nolte, S.; Tünnermann, A. Terahertz line detection by a microlens array coupled photoconductive antenna array. Opt. Express; 2008; 16, pp. 18443-18450. [DOI: https://dx.doi.org/10.1364/OE.16.018443]

11. Nallappan, K.; Li, J.; Guerboukha, H.; Markov, A.; Petrov, B.; Morris, D.; Skorobogatiy, M. A dynamically reconfigurable terahertz array antenna for 2D imaging applications. Proceedings of the 2017 Photonics North (PN); Ottawa, ON, Canada, 6–8 June 2017.

12. Malhotra, I.; Jha, K.R.; Singh, G. Design of highly directive lens-less photoconductive dipole antenna array with frequency selective surface for terahertz imaging applications. Optik; 2018; 173, pp. 206-219. [DOI: https://dx.doi.org/10.1016/j.ijleo.2018.08.010]

13. Chan, W.L.; Moravec, M.L.; Baraniuk, R.G.; Mittleman, D.M. Terahertz imaging with compressed sensing and phase retrieval. Opt. Lett.; 2008; 33, pp. 974-976. [DOI: https://dx.doi.org/10.1364/OL.33.000974]

14. Chan, W.L.; Charan, K.; Takhar, D. A single-pixel terahertz imaging system based on compressed sensing. Appl. Phys. Lett.; 2008; 93, 121105. [DOI: https://dx.doi.org/10.1063/1.2989126]

15. Zhang, Z.; Ma, X.; Zhong, J. Single-pixel imaging by means of Fourier spectrum acquisition. Nat. Commun.; 2015; 6, 6225. [DOI: https://dx.doi.org/10.1038/ncomms7225]

16. Cho, S.H.; Lee, S.-H.; Gung, C.N.; Oh, S.J.; Son, J.H.; Park, H.; Ahn, C.B. Fast terahertz reflection tomography using block-based compressed sensing. Opt. Express; 2011; 19, pp. 16401-16409. [DOI: https://dx.doi.org/10.1364/OE.19.016401]

17. Hwang, B.M.; Lee, S.H.; Lim, W.T.; Ahn, C.B.; Son, J.H.; Park, H. A fast spatial-domain terahertz imaging using block-based compressed sensing. J. Infrared Millim. Terahertz Waves; 2011; 32, pp. 1328-1336. [DOI: https://dx.doi.org/10.1007/s10762-011-9822-5]

18. Ren, X.; Jiang, Y. Spatial domain terahertz image reconstruction based on dual sparsity constraints. Sensors; 2021; 21, 4116. [DOI: https://dx.doi.org/10.3390/s21124116]

19. Buades, A.; Coll, B.; Morel, J.-M. A non-local algorithm for image denoising. Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR); San Diego, CA, USA, 20–25 June 2005; pp. 60-65.

20. Lebrun, M. An analysis and implementation of the BM3D image denoising method. Image Process. Line; 2012; 2, 175. [DOI: https://dx.doi.org/10.5201/ipol.2012.l-bm3d]

21. Zha, Z.; Zhang, X.; Wang, Q.; Tang, L.; Liu, X. Group-based sparse representation for image compressive sensing reconstruction with non-convex regularization. Neurocomputing; 2018; 296, pp. 55-63. [DOI: https://dx.doi.org/10.1016/j.neucom.2018.03.027]

22. Ren, X.; Bai, Y.; Jiang, Y. Hybrid spatial model for fast terahertz imaging. Micromachines; 2021; 12, 1181. [DOI: https://dx.doi.org/10.3390/mi12101181]

23. Fan, X.Y.; Lian, Q.S.; Shi, B.S. Compressed sensing MRI based on image decomposition model and group sparsity. Magn. Reson. Imaging; 2019; 60, pp. 101-109. [DOI: https://dx.doi.org/10.1016/j.mri.2019.03.011] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30910695]

24. Hao, H.X.; Wu, L.D. A new complex valued dictionary learning method for group-sparse representation. Optik; 2019; 196, 163150. [DOI: https://dx.doi.org/10.1016/j.ijleo.2019.163150]

25. Zhang, J.; Zhao, D.; Gao, W. Group-based sparse representation for image restoration. IEEE Trans. Image Process.; 2014; 23, pp. 3336-3351. [DOI: https://dx.doi.org/10.1109/TIP.2014.2323127] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24835225]

26. Goldstein, T.; Osher, S. The split bregman method for L1 regularized problems. SIAM J. Imaging Sci.; 2009; 2, pp. 323-343. [DOI: https://dx.doi.org/10.1137/080725891]

27. Zou, J.; Li, H.F.; Liu, G.Q. Split bregman algorithm for structured sparse reconstruction. IEEE Access; 2018; 6, pp. 21560-21569. [DOI: https://dx.doi.org/10.1109/ACCESS.2018.2825323]

28. Corsaro, S.; Simone, V.D.; Marino, Z. Split bregman iteration for multi-period mean variance portfolio optimization. Appl. Math. Comput.; 2021; 392, 125715. [DOI: https://dx.doi.org/10.1016/j.amc.2020.125715]

29. Beck, A.; Teboulle, M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Imaging Sci.; 2009; 2, pp. 183-202. [DOI: https://dx.doi.org/10.1137/080716542]

30. Xiao, Y.; Zhu, H. A conjugate gradient method to solve convex constrained monotone equations with applications in compressive sensing. J. Math. Anal. Appl.; 2013; 405, pp. 310-319. [DOI: https://dx.doi.org/10.1016/j.jmaa.2013.04.017]

31. Xiao, Y.; Wu, C.; Wu, S.Y. Norm descent conjugate gradient methods for solving symmetric nonlinear equations. J. Glob. Optim.; 2015; 62, pp. 751-762. [DOI: https://dx.doi.org/10.1007/s10898-014-0218-7]

32. Turaga, D.S.; Chen, Y.; Caviedes, J. No reference psnr estimation for compressed pictures. Signal Process. Image Commun.; 2004; 19, pp. 173-184. [DOI: https://dx.doi.org/10.1016/j.image.2003.09.001]

33. Qu, X.B.; Hou, Y.K.; Lam, F.; Guo, D.; Zhong, J.H.; Chen, Z. Magnetic resonance image reconstruction from undersampled measurements using a patch-based nonlocal operator. Med. Image Anal.; 2014; 18, pp. 843-856. [DOI: https://dx.doi.org/10.1016/j.media.2013.09.007]