-
超分辨率图像重建的观测模型如图 1所示。图 1中,PSF指点扩散函数(point spread function)。由于设备和环境的原因,成像系统获得的图像通常是高分辨率图像经过几何变形、模糊、降采样等退化过程之后形成的低分辨率图像[21],该过程可以由下式表示:
$ {\mathit{\boldsymbol{I}}_{{\rm{LR}}}} = (\mathit{\boldsymbol{D}} * {\mathit{\boldsymbol{I}}_{{\rm{HR}}}}) \downarrow + \mathit{\boldsymbol{n}} $
(1) 式中,*表示卷积,ILR为LR图像,IHR为HR图像,D为退化矩阵,↓为降采样,n为成像过程中的各类噪声。
超分辨率图像重建是图 1所示图像退化过程的逆过程。该过程是一个不适定问题,需要引入外部约束和先验信息进行优化求解。由LR图像估计HR图像的表达式为:
$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{{\rm{HR}}}} = {\rm{arg}}{\kern 1pt} {\kern 1pt} \mathop {{\rm{min}}}\limits_{{\mathit{\boldsymbol{I}}_{{\rm{HR}}}}} [({{\left\| {{\mathit{\boldsymbol{I}}_{{\rm{HR}}}} * \mathit{\boldsymbol{D}}) \downarrow - {\mathit{\boldsymbol{I}}_{{\rm{LR}}}}} \right\|}^2} + }\\ {\beta {{\left\| {U({\mathit{\boldsymbol{I}}_{{\rm{HR}}}})} \right\|}^2}]} \end{array} $
(2) 式中,IHR为重建后的图像,U(IHR)为正则项,参量β为加权平衡参量。减小β将减少正则项对图像重建的作用, 当β=0时,该式退化为最大似然估计。基于梯度变换的重建方法将理想图像的梯度先验信息作为正则项(即(2)式中的U(IHR))来约束(2)式,达到锐化图像边缘的目的。
-
基于梯度变换的超分辨率重建方法框架如图 2所示。该方法首先对低分辨率太赫兹图像进行插值,初步提高图像质量,消除锯齿状边缘。插值后的太赫兹图像(Iu)的细节信息包含在其梯度图像(▽Iu)当中。然后提取Iu的梯度图像,通过梯度变换使该梯度图像更加尖锐,变换后的梯度图像(▽${\mathit{\boldsymbol{\hat I}}_{{\rm{HR}}}}$)将包含更多的细节信息。利用▽${\mathit{\boldsymbol{\hat I}}_{{\rm{HR}}}}$作为正则项并根据最大后验概率(maximum a posteriori,MAP)方法约束图像的重建过程可以获得高分辨率图像(IHR),达到锐化图像边缘的目的。▽${\mathit{\boldsymbol{\hat I}}_{{\rm{HR}}}}$越尖锐,IHR的边缘也越尖锐,同时也将包含更多的细节信息。受到太赫兹图像对比度较低以及空间噪声的影响,▽Iu的质量无法满足有效进行梯度变换的要求。本文中引入基于空间信息熵的直方图匹配技术和双边滤波器(bilateral filter, BF)对Iu进行处理,优化其梯度图像。
-
受到成像系统效率的限制,太赫兹原始图像中的像素点较少,边缘不连续。从太赫兹原始图像中提取的梯度图像质量较差,无法对该梯度图像进行梯度变换。因此,本文中首先利用有理分形插值对原始图像进行处理,采用具有较高插值精度的有理插值函数和能够较好地描述图像纹理细节的分形插值方法,扩大图像规模的同时,能够更好地保持图像边缘,从而提高梯度图像的质量。分形插值通过一个特定的迭代函数系统生成分形插值函数而实现,假设原始图像所在的平面区域为Ω=I×J=[a, b]×[c, d],其中I,J分别表示原始图像水平方向和垂直方向像素坐标的集合, 原始图像的像素点为{(xi, yj, zi, j), i=1, …, N; j=1, …, M},记区间Ii=[xi, xi+1],Ji=[yj, yj+1],其中i∈Γ={1, …, N-1},分形插值函数的表达式为:
$ \left\{ {\begin{array}{*{20}{l}} {{\phi _i}(x) = {a_i}x + {b_i}}\\ {{\phi _j}(y) = {c_j}y + {d_j}}\\ {{F_{i,j}}(x,y,z) = {s_{i,j}}z + {q_{i,j}}(x,y)} \end{array}} \right. $
(3) 式中,ϕi(x):I→Ii,ϕj(y):J→Jj为插值区间映射,Fi, j:Ω→R为各插值点的灰度值映射,R为插值后图像的区域,si, j为尺度因子。有理分形插值将(3)式中的qi, j(x, y)定义为:
$ {q_{i,j}}(x,y) = {P_{i,j}}({\phi _i}(x),{\phi _j}(y)) - {s_{i,j}}{B_{i,j}}(x,y) $
(4) 式中,Pi, j(ϕi(x), φj(y))和Bi, j(x, y)分别为双变量有理插值函数和有理扰动基函数[22]。
-
受到成像信号强度的影响,太赫兹原始图像的对比度较低,从插值后的太赫兹图像中提取的梯度图像质量较低,对该梯度图像直接进行梯度变换将无法得到理想的结果。因此,本文中利用基于空间信息熵的直方图匹配技术对插值后的太赫兹图像进行增强,提高对比度,优化梯度图像。考虑到图像灰度级的空间分布, 将图像划分为多个子区域,并在每个子区域根据下式进行图像增强处理,可以实现图像的全局增强。
$ \left\{ {\begin{array}{*{20}{l}} {{I_{k,e}} = {\rm{int}} ({F_k}({I_{{\rm{max}}}} - {I_{{\rm{min}}}}) + {I_{{\rm{min}}}})}\\ {{F_k} = \sum\limits_{l = 1}^k {{f_l}} } \end{array}} \right. $
(5) 式中,k为子区域内的灰度级,Ik, e为利用(5)式对灰度级k进行映射后在经过增强的图像中的灰度级,Fk为灰度级映射函数。fl为一个关于l的离散函数,根据图像的空间信息熵计算[23]。Imin和Imax分别为输入图像最小和最大的灰度值。经过全局增强之后,再对图像进行离散余弦变换并对其变换系数根据下式进行加权实现局部增强,进一步提升图像的对比度。
$ \left\{ {\begin{array}{*{20}{l}} {{{\hat D}_{uv}} = \left( {1 + \frac{{\eta - 1}}{{M - 1}}u} \right)\left( {1 + \frac{{\eta - 1}}{{N - 1}}v} \right){D_{uv}}}\\ {\eta = {{\left( {\sum\limits_{k = 1}^K {{f_k}} {\rm{lo}}{{\rm{g}}_2}{f_k}} \right)}^\gamma }} \end{array}} \right. $
(6) 式中,K为子区域内的灰度级数量,参量γ用于平衡全局和局部增强的程度,其取值范围为,本文中该参量取0.25以保证图像的增强效果,同时不会由于过度增强图像而放大其中的噪声。
-
由于成像过程中太赫兹信号存在衰减且不同采样点采集到的信号衰减程度不同,同时受到成像系统其它噪声(如探测器随机噪声、扫描系统电磁噪声等)的影响,太赫兹图像中的空间噪声将对梯度图像的质量产生较大影响。因此,本文中采用双边滤波器来消除空间噪声的影响,该方法的优点在于其能够在保持图像边缘的前提下对图像进行平滑,不会破坏图像的边缘梯度轮廓从而不影响后续的梯度变换,双边滤波器[24]的表达式为:
$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{I}}_{BF}}{|_{i,j}} = \\ \frac{{\sum\limits_{s = - W}^W {\sum\limits_{t = - W}^W {{\rm{exp}}} } \left( { - \frac{{{s^2} + {t^2}}}{{2\sigma _{\rm{q}}^2}}} \right){\rm{exp}}\left[ { - \frac{{{{({\mathit{\boldsymbol{I}}_{i - s,j - t}} - {\mathit{\boldsymbol{I}}_{i,j}})}^2}}}{{2\sigma _{\rm{r}}^2}}} \right]{\mathit{\boldsymbol{I}}_{i - s,j - t}}}}{{\sum\limits_{s = - W}^W {\sum\limits_{t = - W}^W {{\rm{exp}}} } \left( { - \frac{{{s^2} + {t^2}}}{{2\sigma _{\rm{q}}^2}}} \right){\rm{exp}}\left[ { - \frac{{{{({\mathit{\boldsymbol{I}}_{i - s,j - t}} - {\mathit{\boldsymbol{I}}_{i,j}})}^2}}}{{2\sigma _{\rm{r}}^2}}} \right]}} \end{array} $
(7) 式中,W=3σq,i和j为图像中像素点的坐标;σq和σr为几何扩散和广度扩散,由实验确定,分别用于控制空间域和值域权值的衰减[24-25],本文中取σq=3,σr=20。
-
自适应梯度变换(adaptive gradient transform,AGT)利用双边高斯模型拟合图像边缘[20],构造“差函数”形式的梯度变换优化梯度图像,使其更加尖锐,自适应梯度变换表达式为:
$ \left\{ {\begin{array}{*{20}{l}} {▽ {{\mathit{\boldsymbol{\hat I}}}_{{\rm{HR}}}}{|_{i,j}} = {{\left. {\frac{{{\sigma _{{\rm{LR}}}}}}{{{\sigma _{{\rm{HR}}}}}}\left( {▽ {\mathit{\boldsymbol{I}}_{{\rm{LR}}}} - \frac{{{{({\sigma _{{\rm{LR}}}})}^2}}}{{\hat \sigma }}\left| {\frac{{\partial (▽ {\mathit{\boldsymbol{I}}_{{\rm{LR}}}})}}{{\partial n}}} \right|} \right)} \right|}_{i,j}}}\\ {\hat \sigma = \frac{{{{({\sigma _{{\rm{LR}}}})}^2}{{({\sigma _{{\rm{HR}}}})}^2}}}{{{{({\sigma _{{\rm{LR}}}})}^2} - {{({\sigma _{{\rm{HR}}}})}^2}}}} \end{array}} \right. $
(8) 式中,σLR和σHR分别为LR图像的梯度图像和变换后梯度图像中梯度轮廓的锐度(拟合高斯模型的标准差),n为梯度的方向。σLR和σHR之间的关系[18]为:
$ {\sigma _{{\rm{HR}}}} = [1 - {\rm{exp}}( - \mu {\sigma _{{\rm{LR}}}})]{\sigma _{{\rm{LR}}}} $
(9) 式中,参量σLR利用LR图像中每个像素点及其八邻域内的信息计算[20]。
-
根据MAP方法,利用优化后的梯度图像(▽${\mathit{\boldsymbol{\hat I}}_{{\rm{HR}}}}$)作为正则约束条件代入(2)式可以获得高分辨率太赫兹图像(IHR)[10, 18],其表达式为:
$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{{\rm{HR}}}} = {\rm{arg}}{\kern 1pt} \mathop {{\kern 1pt} {\rm{min}}}\limits_{{\mathit{\boldsymbol{I}}_{{\rm{HR}}}}} [{{\left\| {({\mathit{\boldsymbol{I}}_{{\rm{HR}}}} * \mathit{\boldsymbol{D}}) \downarrow - {\mathit{\boldsymbol{I}}_{{\rm{LR}}}}} \right\|}^2} + }\\ {\beta {{\left\| {\nabla {\mathit{\boldsymbol{I}}_{{\rm{HR}}}} - \nabla {{\mathit{\boldsymbol{\hat I}}}_{{\rm{HR}}}}} \right\|}^2}]} \end{array} $
(10) 式中,▽IHR为重建后的图像的梯度图像。本文中,参量β取值为0.15。考虑到太赫兹光束的强度服从高斯分布,因此重建时所采用的退化矩阵D为高斯矩阵。利用梯度下降算法[18-19]优化(10)式可以获得高分辨率太赫兹图像。
-
实验中采用的太赫兹透射式成像系统(Zomega Terahertz Corporation(USA))如图 3所示。其中M1~M10为反射镜,OAP1和OAP2为离轴抛物面镜(off axis parabolic mirror, OAP),L1和L2为聚乙烯透镜。飞秒激光器产生的激光光束被偏振分束器(polarized beam splitter, PBS)分为抽运光和探测光。半波板(half-wave-plate, HWP)用于改变抽运光和探测光的比例。延线由两个反射镜构成,用于调节抽运光和探测光之间的时间延迟。抽运光经过衰减器,同时被反射镜反射后照射到砷化镓(GaAs)光电导天线上产生太赫兹波。探测光被用于探测太赫兹波,氧化铟锡(indium tin oxide, ITO)导电玻璃和平衡探测器(balance detector, BD)共同构成成像系统的探测装置。测量过程中,系统的温度维持在21℃(294K)。成像系统工作频率范围为0.1THz~1.0THz。样品为一个金属制斩波器,其半径为5cm,原始图像的大小为55pixle×55pixle,插值后的太赫兹图像的大小为440pixle×440pixle。样品被放入遮蔽物中并被放置于2维电动平移台上,通过控制程序移动平移台实现对样品的逐点扫描。扫描过程中,探测器对每一个扫描点测量太赫兹光谱,再将提取的太赫兹光谱按扫描次序排列成一幅太赫兹图像。
Figure 3. The THz transmission imaging system a—the illustration of THz imaging system b—the experimental setup of the THz imaging system
探测器在非目标区域采集到的时域信号和频谱分别如图 4a、图 4b所示。在放置遮蔽物后,太赫兹信号会存在一定的时延和衰减。经过傅里叶变换之后,对频域信号进行滤波即可获得单频信号。本文中使用0.25THz,0.50THz,0.75THz图像来分析成像信号频率对重建性能的影响。
-
本文中采用离散梯度算子[-1/2, 0, 1/2]和[-1/2, 0, 1/2]T与图像进行卷积,分别获得每个像素点水平方向和垂直方向的梯度。以0.75THz图像为例,图 5a为太赫兹原始图像的梯度图像,图 5b为插值后的太赫兹图像的梯度图像,图 5c为对插值后的太赫兹图像直接进行增强后(未引入双边滤波)提取的梯度图像,图 5d为对插值后的太赫兹图像进行增强和双边滤波后提取的梯度图像,图 5e为变换后的梯度图像。图 5中白色线条为太赫兹图像的边缘梯度轮廓,白色线条越窄,越明亮,梯度的幅值越大,图像边缘越尖锐。
由图 5a、图 5b可知,经过插值之后,部分图像细节得到恢复,边缘变得连续;由图 5c可知,直接对经过插值的太赫兹图像进行增强虽然可以提高其边缘梯度轮廓的幅值,但也会放大图像中的噪声,从而降低梯度图像的质量;由图 5d可知,引入基于空间信息熵的直方图匹配技术和双边滤波器对经过插值的太赫兹图像进行处理后,其梯度图像的质量大幅提高,从而可以使得梯度变换的效果更加理想;另外,从图 5e中可以看出,梯度变换可以有效压缩梯度轮廓的范围(白色线条更细),提高梯度的幅值。
-
经过插值处理之后,太赫兹图像的分辨率得到提高(如图 6b、图 6g和图 6l所示),部分细节信息更加清晰。相比于0.25THz的太赫兹信号,频率为0.50THz和0.75THz的信号的光斑较小,因此0.50THz和0.75THz图像经过插值后边缘更清晰。总体而言,插值处理虽然恢复了太赫兹图像中的部分细节信息,增强了图像边缘的连续性,但插值并不能解决由于衍射极限的影响而产生的模糊问题。因此,在对图像进行插值之后,还需要对图像进行重建,锐化图像边缘。
利用有理分形插值,Wiener滤波,Lucy-Richardson迭代和基于梯度变换的方法对本实验中所采用的成像系统获得的太赫兹图像进行复原的结果如图 6所示。对经过插值的太赫兹图像采用基于梯度变换的方法重建之后,其边缘更加尖锐, 同时图像中的空间噪声并没有因为边缘锐化而被放大且没有出现Lucy-Richardson迭代算法处理后出现的振铃现象(如图 6d、图 6i、图 6n所示)。其中0.75THz图像边缘最尖锐,这是因为3幅图像中获取0.75THz图像的太赫兹信号频率最高,波长最短,光斑最小。但由于其信号强度较弱,0.75THz图像的对比度不及0.25THz和0.50THz图像。由图 6e可见,0.25THz图像边缘锐化的效果最明显,这是由于受到衍射效应的影响,和经过插值的0.50THz,0.75THz图像相比,插值后的0.25THz图像的边缘较为模糊,因此, 该方法在对0.25THz图像质量的提升更加明显。
本文中采用边缘强度(ICV)和平均梯度(Idef)两种评价标准[26]定量评价图像重建的效果。两种评价标准的计算式分别为:
$ {{I_{{\rm{CV}}}} = \frac{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N | } {\nabla ^2}\mathit{\boldsymbol{I}}|}}{{MN}}} $
(11) $ {{I_{{\rm{def}}}} = \frac{{\sum\limits_{i = 1}^{M - 1} {\sum\limits_{j = 1}^{N - 1} | } \nabla \mathit{\boldsymbol{I}}|}}{{MN}}} $
(12) 式中,▽2为拉普拉斯算子,M和N表示图像的大小。
经过插值的太赫兹图像和重建后的图像的边缘强度和平均梯度计算结果如表 1所示。由表 1中的数据可知,该重建方法可以有效提高太赫兹图像的边缘强度和平均梯度。其中该方法对0.25THz图像的边缘强度的提升效果最为明显。0.50THz图像经过重建之后,边缘强度和平均梯度最高,这是因为和0.25THz图像相比,其成像信号频率较高,光斑较小,边缘较为清晰。和0.75THz图像相比,虽然0.50THz图像成像信号频率较低,但其信号强度更大,因此图像的对比度更好。
Table 1. ICV and Idef of interpolated and reconstructed images
frequency ICV Idef Iu Wiener Lucy-Richardson proposed Iu Wiener Lucy-Richardson proposed 0.25THz 2.99 3.69 5.68 8.04 5.54 6.59 6.18 6.43 0.50THz 5.66 5.60 10.08 12.22 6.95 6.37 8.05 8.91 0.75THz 4.45 7.55 5.75 9.06 5.44 5.65 5.25 6.77 综上所述,利用基于梯度变换的超分辨率重建方法可以有效恢复低分辨率太赫兹图像中的细节信息,锐化图像边缘,重建高分辨率太赫兹图像。由于需要提取较为清晰的梯度图像进行梯度变换,该方法对于空间噪声较小以及图像对比度较高的太赫兹原始图像将有更好的重建效果。
基于梯度变换的太赫兹图像超分辨率重建
Super-resolution reconstruction for terahertz images based on gradient transform
-
摘要: 为了提高太赫兹图像的质量,克服边缘模糊的缺陷,采用有理分形插值和基于梯度变换的图像超分辨率重建算法相结合的方法对0.25THz,0.50THz和0.75THz图像进行超分辨率重建实验,并对实验结果进行了定量分析,利用基于空间信息熵的直方图匹配技术和双边滤波器对重建算法进行了优化,增强了该方法的适用性。结果表明,对经过插值的太赫兹图像采用基于梯度变换的超分辨率重建方法处理之后,0.25THz,0.50THz和0.75THz图像的边缘强度分别提高了169%,116%和104%,平均梯度分别提高了16%,28%和24%;同时,成像信号频率和强度也会对重建性能产生影响。该方法可以有效恢复太赫兹图像当中的细节信息,锐化图像边缘,提高图像质量且不会出现振铃现象,具有较好的实用价值。Abstract: In order to improve the quality of terahertz image and overcome the problem of edge blur of terahertz image, a super-resolution reconstruction method, which combines rational fractal interpolation and gradient field transform, was proposed in this paper for terahertz image reconstruction with frequencies of 0.25THz, 0.50THz, and 0.75THz. Meanwhile, spatial entropy-based image enhancement and bilateral filtering are introduced to optimize the reconstruction. Experimental results illustrate that after processing interpolated terahertz images by super-resolution reconstruction based on gradient field transform, the edge strength of images with frequencies of 0.25THz, 0.50THz, and 0.75THz were respectively improved by 169%, 116%, and 104%, and the average gradient of those images were improved by 16%, 28%, and 24%, respectively. Moreover, signal frequency and intensity would affect the performance of the reconstruction. This method can recover the detail information in terahertz images, sharpen edges of the object, improve the quality of terahertz images without ringing effect, and has practical value.
-
Key words:
- image processing /
- terahertz /
- super-resolution /
- interpolation /
- gradient field transformation /
- edge sharpening
-
a—original image(0.25THz) b—interpolated image(0.25THz) c—reconstructed image(0.25THz, Wiener) d—reconstructed image(0.25THz, Lucy-Richardson) e—reconstructed image(0.25THz, proposed) f—original image(0.50THz) g—interpolated image (0.50THz) h—reconstructed image(0.50THz, Wiener) i—reconstructed image(0.50THz, Lucy-Richardson) j—reconstructed image(0.50THz, proposed) k—original image(0.75THz) l—interpolated image(0.75THz) m—reconstructed image(0.75THz, Wiener) n—reconstructed image(0.75THz, Lucy-Richardson) o—reconstructed image(0.75THz, proposed)
Figure 6. The experimental results of terahertz image reconstruction
Table 1. ICV and Idef of interpolated and reconstructed images
frequency ICV Idef Iu Wiener Lucy-Richardson proposed Iu Wiener Lucy-Richardson proposed 0.25THz 2.99 3.69 5.68 8.04 5.54 6.59 6.18 6.43 0.50THz 5.66 5.60 10.08 12.22 6.95 6.37 8.05 8.91 0.75THz 4.45 7.55 5.75 9.06 5.44 5.65 5.25 6.77 -
[1] MITTLEMAN D M. Twenty years of terahertz imaging[J]. Optics Express, 2018, 26(8): 9417-9431. doi: 10.1364/OE.26.009417 [2] YANG X, ZHAO X, YANG K, et al. Biomedical applications of te-rahertz spectroscopy and imaging[J]. Trends in Biotechnology, 2016, 34(10): 810-824. doi: 10.1016/j.tibtech.2016.04.008 [3] OH S J, KIM S H, JI Y B, et al. Study of freshly excised brain ti-ssues using terahertz imaging[J]. Biomedical Optics Express, 2014, 5(8): 2837-2842. doi: 10.1364/BOE.5.002837 [4] STANTCHEV R I, SUN B, HORNETT S M, et al. Noninvasive, near-field terahertz imaging of hidden objects using a single-pixel detector[J]. Science Advances, 2016, 2(6): e1600190. doi: 10.1126/sciadv.1600190 [5] DONG J, LOCQUET A, CITRIN D. Depth resolution enhancement of terahertz deconvolution by autoregressive spectral extrapolation[J]. Optics Letters, 2017, 42(9): 1828-1831. doi: 10.1364/OL.42.001828 [6] PARK C, PARK J Y, SON J H, et al. Terahertz imaging of excised oral cancer at frozen temperature[J]. Biomedical Optics Express, 2013, 4(8): 1413-1421. doi: 10.1364/BOE.4.001413 [7] XU L M, FAN W H, LIU J. De-noising and enhancement for terahertz imaging[J]. Infrared and Laser Engineering, 2013, 42(10): 2865-2870(in Chinese). [8] CHEN W Y, XIE Q, WU H. Applied research on supper resolution reconstruction method for image enhancement[J].Journal of Naval University of Engineering, 2015, 27(1): 74-78(in Chinese). [9] SHEN H F, LI P X, ZHANG L P, et al. Overview on super resolution image reconstruction[J]. Optical Technique, 2009, 35(2): 194-203(in Chinese). [10] YUE L, SHEN H, LI J, et al. Image super-resolution: The techniques, applications, and future[J]. Signal Processing: Image Communication, 2016, 128: 389-408. [11] XIE Y Y, HU Ch H, SHI B, et al. Super-resolution image reconstruction and its application in terahertz images[J]. System Simulation Technology, 2013, 9(4):306-309(in Chinese). [12] ZHENG X T, YUAN Y, LU X Q. Single image super-resolution re-storation algorithm from external example to internal self-similarity[J].Acta Optica Sinica, 2017, 37(3): 0318006(in Chinese). doi: 10.3788/AOS201737.0318006 [13] XU L M, FAN W H, LIU J. High-resolution reconstruction for terahertz imaging[J]. Applied Optics, 2014, 53(33): 7891-7897. doi: 10.1364/AO.53.007891 [14] SHI J, WANG Y, XU D, et al. Terahertz imaging based on morphological reconstruction[J]. IEEE Journal of Selected Topics in Quantum Electronics, 2017, 23(4): 1-7. [15] AHI K. A method and system for enhancing the resolution of terahertz imaging[J]. Measurement, 2019, 138: 614-619. doi: 10.1016/j.measurement.2018.06.044 [16] WEI M G, LIANG D Ch, GU J Q, et al. Terahertz radar imaging based on time-domain spectroscopy[J]. Journal of Radars, 2015, 4(2): 222-229(in Chinese). [17] FATTAL R. Image upsampling via imposed edge statistics[J]. ACM Transactions on Graphics, 2007, 26(3): 95. doi: 10.1145/1276377.1276496 [18] SUN J, SUN J, XU Z B, et al. Gradient profile prior and its applications in image super-resolution and enhancement[J]. IEEE Transactions on Image Processing, 2011, 20(6): 1529-1542. doi: 10.1109/TIP.2010.2095871 [19] YAN Q, XU Y, YANG X, et al. Single image superresolution based on gradient profile sharpness[J]. IEEE Transactions on Image Processing, 2015, 24(10): 3187-3202. doi: 10.1109/TIP.2015.2414877 [20] QIANG S, XIONG R, DONG L, et al. Fast image super-resolution via local adaptive gradient field sharpening transform[J]. IEEE Transactions on Image Processing, 2018, 27(4): 1966-1980. doi: 10.1109/TIP.2017.2789323 [21] PARK S C, PARK M K, KANG M G. Super-resolution image reconstruction: A technical overview[J]. IEEE Signal Processing Magazine, 2003, 20(3): 21-36. doi: 10.1109/MSP.2003.1203207 [22] ZHANG Y, FAN Q, BAO F, et al. Single-image super-resolution based on rational fractal interpolation[J]. IEEE Transactions on Image Processing, 2018, 27(8): 3782-3797. doi: 10.1109/TIP.2018.2826139 [23] CELIK T. Spatial entropy-based global and local image contrast enhancement[J]. IEEE Transactions on Image Processing, 2014, 23(12): 5298-5308. doi: 10.1109/TIP.2014.2364537 [24] CHAUDHURY K N, DABHADE S D. Fast and provably accurate bilateral filtering[J]. IEEE Transactions on Image Processing, 2016, 25(6): 2519-2528. doi: 10.1109/TIP.2016.2548363 [25] WEI X, LIU Q, GUO Y Zh, et al. Image denoising algorithm based on joint bilateral filter and multi-resolution analysis[J]. Computer Engineering and Design, 2016, 37(12): 3327-3833(in Chinese). [26] ZHANG X, ZHAO Y M, DENG Ch, et al. Study on the passive te-rahertz image target detection[J].Acta Optica Sinica, 2013, 33(2): 0211002(in Chinese). doi: 10.3788/AOS201333.0211002