CN107525588B - 一种基于gpu的双相机光谱成像***的快速重构方法 - Google Patents
一种基于gpu的双相机光谱成像***的快速重构方法 Download PDFInfo
- Publication number
- CN107525588B CN107525588B CN201710699197.6A CN201710699197A CN107525588B CN 107525588 B CN107525588 B CN 107525588B CN 201710699197 A CN201710699197 A CN 201710699197A CN 107525588 B CN107525588 B CN 107525588B
- Authority
- CN
- China
- Prior art keywords
- image
- camera
- reconstruction
- gpu
- dual
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 61
- 238000000701 chemical imaging Methods 0.000 title claims abstract description 49
- 238000005457 optimization Methods 0.000 claims abstract description 44
- 238000001228 spectrum Methods 0.000 claims abstract description 20
- 238000011478 gradient descent method Methods 0.000 claims abstract description 7
- 230000006870 function Effects 0.000 claims description 72
- 230000003595 spectral effect Effects 0.000 claims description 53
- 238000004422 calculation algorithm Methods 0.000 claims description 46
- 239000011159 matrix material Substances 0.000 claims description 41
- 238000005070 sampling Methods 0.000 claims description 25
- 239000006185 dispersion Substances 0.000 claims description 18
- 230000004044 response Effects 0.000 claims description 16
- 238000004364 calculation method Methods 0.000 claims description 13
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 claims description 9
- 230000009977 dual effect Effects 0.000 claims description 9
- 238000012545 processing Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 8
- 238000013178 mathematical model Methods 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 5
- 230000003190 augmentative effect Effects 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 3
- 230000008602 contraction Effects 0.000 claims description 3
- 238000011423 initialization method Methods 0.000 claims description 3
- 238000005316 response function Methods 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 2
- 238000003384 imaging method Methods 0.000 abstract description 8
- 238000011160 research Methods 0.000 abstract description 5
- 238000012360 testing method Methods 0.000 description 6
- 230000001133 acceleration Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 239000011324 bead Substances 0.000 description 4
- 230000006835 compression Effects 0.000 description 4
- 238000007906 compression Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 230000006978 adaptation Effects 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000002612 dispersion medium Substances 0.000 description 1
- 125000001495 ethyl group Chemical group [H]C([H])([H])C([H])([H])* 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- ONCZDRURRATYFI-QTCHDTBASA-N methyl (2z)-2-methoxyimino-2-[2-[[(e)-1-[3-(trifluoromethyl)phenyl]ethylideneamino]oxymethyl]phenyl]acetate Chemical compound CO\N=C(/C(=O)OC)C1=CC=CC=C1CO\N=C(/C)C1=CC=CC(C(F)(F)F)=C1 ONCZDRURRATYFI-QTCHDTBASA-N 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/2823—Imaging spectrometer
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/2803—Investigating the spectrum using photoelectric array detector
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/2803—Investigating the spectrum using photoelectric array detector
- G01J2003/282—Modified CCD or like
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J2003/283—Investigating the spectrum computer-interfaced
Landscapes
- Physics & Mathematics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Physics & Mathematics (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Investigating, Analyzing Materials By Fluorescence Or Luminescence (AREA)
Abstract
本发明公开的一种基于GPU的双相机光谱成像***的快速重构方法,涉及能够快速获取高分辨率高光谱图像的方法,属于计算摄像学领域。本发明应用于基于编码孔径快照光谱成像+灰度相机的双相机光谱成像***,将高光谱图像重构问题转化为多个子优化问题,并使用GPU完成各个子问题的求解:使用cuBLAS库和共轭梯度下降法更新高光谱图像;使用软阈值函数更新辅助变量;重复迭代完成高光谱图像的重构。本发明能够高质量地完成双相机光谱成像***的高光谱图像重构,在保证重建结果具备高空间分辨率和高光谱保真性的同时,大幅度提高高光谱图像重建的效率,扩展高光谱图像的应用范围。本发明可用于载人航天、地质勘测和植被研究等多个领域。
Description
技术领域
本发明专利涉及一种用于双相机光谱成像的高光谱图像重构方法,尤其是能够快速获取高分辨率高光谱图像的方法,属于计算摄像学领域。
背景技术
高光谱成像技术是将空间成像技术和频谱成像技术相结合,对目标场景在频谱上进行几十甚至上百个波段的连续测量。该技术得到的图像包含目标场景的二维空间信息和一维光谱信息,被称为数据立方体。相比传统彩色成像,高光谱成像技术能够获得内容更为丰富,细节更为显著的有用信息。该技术已经被应用于植被研究、大气检测,载人航天、医学诊断等多个领域。
近年来,基于压缩感知的计算光谱成像技术被广泛运用到高光谱成像。与传统光谱成像***相比,计算光谱成像能够获得更高空间分辨率、光谱分辨率和时间分辨率的高光谱图像,具有更为广阔的应用前景。Ashwin Wagadarikar等人提出的编码孔径快照光谱成像仪(Coded Aperture Snapshot Spectral Imager,CASSI)采用编码孔径和色散介质来对目标场景进行调制,通过探测器获取三维图谱数据的二维压缩投影。双相机光谱成像***则是在CASSI***的基础上,加入了灰度相机通道。该***在获取目标场景的二维压缩光谱投影的同时,还获取目标场景的二维灰度投影。与最初的CASSI***相比,双相机***能够最大化入射光利用率,大幅度提高高光谱图像的成像质量,因此已成为该领域的重点研究对象。
目前,从双相机***重构出原始的三维高光谱图像的方法主要有两类:基于全变差最小化的两步迭代收缩/阈值算法(Two-Step Iterative Shrinkage/ThresholdingAlgorithms,TwIST)和基于自适应字典学习的重构算法。TwIST算法建立在迭代收缩阈值算法和迭代加权收缩算法的基础上,其核心是使用反向投影函数去噪,并利用前两次的结果进行迭代更新,从而完成高光谱图像的重建。由于其计算复杂度很高,TwIST算法重构高光谱图像时间很长,完全不能满足快速重建的要求。同时,重构的高光谱图像分辨率也有待提高。
基于自适应字典学习的重构算法利用光谱图像和灰度图像在空间维的相似性,采用字典学习和稀疏表示的方法完成高光谱图像的重构。这种方法虽然能够获得较好的重构质量,但是每次重构都需要重新学习字典,导致重构时间比TwIST算法更长,因此也不能满足快速重建的要求,极大地限制了高光谱图像的应用和发展。
CUDA(Compute Unified Device Architecture),是显卡厂商NVIDIA推出的通用并行计算架构,该架构使用图像处理器GPU作为计算核心完成复杂的计算难题。相比传统CPU处理器,GPU拥有更多的计算单元,能够快速完成运算量大的计算任务,且已经在深度学习和人工智能等多个领域得到应用。
发明内容
针对现有重构算法存在的重构时间长、图像质量低等问题。本发明要解决的技术问题是提供一种基于GPU的双相机光谱成像***的快速重构方法,能够快速完成双相机光谱成像***的重构,具有重构速度快,成像质量高等优点。
为达到以上目的,本发明采用以下技术方案:
本发明公开的一种基于GPU的双相机光谱成像***的快速重构方法,应用于基于编码孔径快照光谱成像+灰度相加的双相机光谱成像***,将高光谱图像重构问题转化为多个子优化问题,并使用图形处理器GPU完成各个子问题的求解:使用共轭梯度下降法更新高光谱图像;使用软阈值函数更新中间变量;重复迭代完成高光谱图像的重构。
本发明公开的一种基于GPU的双相机光谱成像***的快速重构方法,包括以下步骤:
步骤101:输入双相机光谱成像仪的采样图像Y、标定后的前向响应矩阵H、主函数最大迭代次数Maxiters和正则化系数τ。
步骤101中所述双相机光谱成像仪为基于编码孔径快照光谱成像仪(CodedAperture Snapshot Spectral Imager,CASSI)+灰度相机的双相机采样光谱成像仪。双相机光谱成像仪主要由分光镜、物镜、编码模版、中继镜、色散棱镜和灰度相机等部件构成。目标场景的高光谱图像X大小为M×N×Ω,高光谱图像X上任意一点的像素值为x(i,j,λ),1≤i≤M,1≤j≤N,1≤λ≤Ω。其中,M×N表示高光谱图像的空间分辨率,Ω表示高光谱图像的频段数。入射光首先到达分光镜将其能量一分为二,一半进入编码孔径快照光谱成像仪CASSI分支,一半进入灰度相机分支。进入编码孔径快照光谱成像仪CASSI分支的光会到达编码模版进行0-1编码。经编码后的图像到达色散棱镜后,不同频段的图像会沿着竖直方向进行偏移。最后所有频段的图像到达灰度相机后进行叠加,得到压缩的二维混叠光谱图像。编码孔径快照光谱成像仪CASSI成像仪的数学模型为:
其中ω(λ)表示CCD相机的频谱响应函数,Cu(i,j)表示编码模版函数,φ(λ)表示色散棱镜的波段偏移函数,yc(i,j)为二维压缩光谱采样图像,v1(i,j)表示高斯白噪声。
进入灰度相机分支的入射光会直接到达灰度相机得到目标场景的二维灰度投影,数学模型为:
将式(1)和式(2)联立,并写成矩阵的形式为:
Y=HX+V (3)
其中Y=[Yc;Yp]表示二维压缩光谱采样图像+二维灰度图像。X表示三维数据立方体,V=[Vc;Vp]表示高斯白噪声。H表示对编码孔径快照光谱成像仪CASSI+灰度相机标定后的前向响应矩阵,为编码模版函数、色散棱镜偏移函数和CCD灰度相机频谱响应的综合作用。
步骤102:申请主机CPU内存,将前向矩阵H和二维压缩光谱采样图像Y读入到主机CPU内存。
步骤102所述申请主机CPU内存优选C语言的malloc内存分配函数,申请的主机CPU内存为一维线性内存,方便后续处理中将CPU内存读入图像处理器GPU显存的操作。
步骤102所述将矩阵H和Y读入主机CPU内存的方式是根据编码孔径光谱仪的色散棱镜偏移方向确定:如果色散棱镜是沿竖直方向偏移,则以二维采样图像的水平方向为基准进行读入,即逐行读入;如果色散棱镜是沿水平方向偏移,则以二维采样图像的竖直方向为基准进行读入,即逐列读入。目的是方便在后续处理中编码孔径快照光谱成像仪CASSI***前向模型的计算,提高计算效率。
步骤103:为输入输出、以及后续计算过程所需要的数据申请GPU显存,将二维压缩光谱图像Y和前向响应矩阵H拷贝到GPU显存中。设置CUDA并行计算时需用的线程格/块数。
步骤103所述申请图形处理器GPU显存的方法优选通用并行计算架构CUDA的cudaMalloc函数,拷贝到图形处理器GPU显存的方法优选通用并行计算架构CUDA的cudaMemcpy函数。
步骤104:初始化重构高光谱图像X0,优化目标函数的初始值f0,以及辅助矩阵S和U,辅助因子ρ。初始化当前迭代次数t=0。
步骤104所述重构高光谱图像X0的初始化方式为:
X0=HTY (4)
其中HT表示编码孔径快照光谱成像仪CASSI***前向响应H的转置,即从二维压缩光谱图像+二维灰度图像反演到三维数据立方体的过程。
步骤104所述优化目标函数为全局优化目标函数。根据自然图像的平滑特性,将高光谱重构问题转化为基于全变差约束的最优化问题,从而得到全局优化目标函数为:
||DX||1=Σi,j,λ|x(i,j+1,λ)-x(i,j,λ)|+|x(i+1,j,λ)-x(i,j,λ)| (6)
矩阵D为差分矩阵,它和图像X的作用结果如下:
公式(7)的结果为两个大小和X相同的矩阵,分别表示图像X在水平方向和竖直方向上的差分值。
步骤104所述辅助矩阵S和U是指求解优化方程(4)时所需要的辅助变量,大小均为2×M×N×Ω,初始化为全零矩阵。公式(5)无法直接求解,将公式(5)转化为多个子优化问题的求解,原理如下:
引入辅助矩阵S=DX,得到如下带约束的优化方程:
其增广拉格朗日方程为:
其中ρ为辅助因子,优选2.048。当式(9)最小时,变量X、S的值即为优化问题式(8)的解,因此得到两个子优化问题,分别为:
接下来,对X、S和U进行交替更新、循环迭代即可完成高光谱图像的重构。
步骤105:使用共轭梯度下降算法更新高光谱图像Xt+1。
针对优化问题式(10),X的最小二乘解为:
Xt+1=(HTH+ρDTD)-1(HTY+DT(Ut+ρSt)) (12)
由于矩阵H和矩阵D的规模很大,无法直接求其解析解,因此需利用共轭梯度下降法求高光谱图像Xt+1的近似解,从而完成高光谱图像Xt+1的更新。
同时,为了加速程序的执行效率,借助NVIDIA提供的矩阵线性运算库cuBLAS来实现共轭梯度下降法,从而完成高光谱图像Xt+1的更新。
步骤106:更新辅助矩阵S和U。
根据优化问题式(11),S的最小二乘解为:
辅助矩阵U的更新公式如下:
Ut+1=Ut+ρ(St-DXt) (14)
Ut+1的更新方式优选矩阵线性运算库cuBLAS库,可以提高程序效率。
步骤107:使用更新后的高光谱图像Xt+1计算全局优化目标函数值ft+1。
步骤107所述全局优化目标函数与步骤103的相同,即:
步骤108:根据步骤107计算的结果执行迭代选择策略,完成双相机光谱成像***的快速重构。
步骤108所述迭代选择策略如下:如果收敛,即步骤107的目标函数值ft+1小于上一个目标函数值ft,则计算目标函数值的相对变化量:
若不满足迭代停止条件,即Tol大于预设阈值或者当前迭代次数t小于最大迭代次数Maxiters,则更新辅助因子ρ=1.05×ρ,当前迭代次数t=t+1,目标函数值ft=ft+1,并转至步骤105进行迭代;否则,停止迭代并输出最后一次更新的高光谱图像,从而完成双相机光谱成像***的快速重构。
如果不收敛,即步骤107的目标函数值ft+1大于上一个目标函数值ft,则更新辅助因子ρ=2.0×ρ。判断更新后的ρ值,若ρ小于预设阈值,则转至步骤105进行迭代;否则,停止迭代并输出最后一次更新的高光谱图像,从而完成双相机光谱成像***的快速重构。
有益效果:
1、本发明公开的一种基于GPU的双相机光谱成像***的快速重构方法,将高光谱重建问题转化为多个子优化问题进行求解,收敛速度和重构效率得到很大提升。
2、本发明公开的一种基于GPU的双相机光谱成像***的快速重构方法,使用交替迭代的方式求解优化问题,每一个更新步骤都能够达到最优解,因此能够提高空间分辨率和光谱保真度。
3、本发明公开的一种基于GPU的双相机光谱成像***的快速重构方法,使用GPU完成高光谱图像的重建过程,并利用cuBLAS库来提高GPU的利用率,大幅度降低双相机光谱成像仪的重建时间。
4、本发明公开的一种基于GPU的双相机光谱成像***的快速重构方法,使用预先统一分配的方式管理GPU显存,最小化计算所需显存的大小,对GPU的要求低,使用入门级显卡就能够获得非常显著的加速效果。
附图说明
图1是本发明中用于双相机光谱成像的***结构图;
图2是本发明公开的基于GPU的双相机光谱成像***的快速重构方法的总流程图;
图3是本发明和对比算法对测试图片beads进行仿真重构后波长为440nm时的对比图,其中:图3-a为本发明重构的结果,图3-b为对比算法重构的结果。
图4是本发明和对比算法对测试图片beads进行仿真重构后波长为540nm时的对比图,其中:图4-a为本发明重构的结果,图4-b为对比算法重构的结果。
图5是本发明和对比算法对测试图片beads进行仿真重构后波长为640nm时的对比图,其中:图5-a为本发明重构的结果,图5-b为对比算法重构的结果。
具体实施方式
为了更好地说明本发明的目的和优点,下面结合附图和实例对发明内容做进一步说明。
实施例1:
本实施例公开的一种基于GPU的双相机光谱成像***的快速重构方法,应用于基于编码孔径快照光谱成像+灰度相机的双相机光谱成像***。编码孔径快照光谱成像***,即CASSI***,最早由Ashwin Wagadarikar等人提出(详见Wagadarikar A,John R,WillettR,Brady D.Single disperser design for coded aperture snapshot spectralimaging[J].Applied optics.2008,47(10):B44-B51.)。CASSI***利用编码模版将光谱信息在空间维上进行随机调制,并使用色散棱镜在光谱维上进行压缩,从而实现了由三维数据立方体到二维数据的压缩。双相机光谱成像***则是在CASSI***的基础上,加入了灰度相机通道,其***结构如图1所示(详见Wang L,Xiong Z,Gao D,et al.Dual-cameradesign for coded aperture snapshot spectral imaging[J].Applied Optics,2015,54(4):848-58.)。双相机光谱成像***在获取目标场景的二维压缩光谱投影的同时,还获取目标场景的二维灰度投影。与CASSI***相比,双相机***能够最大化入射光利用率,同时大幅度提高高光谱图像的成像质量,因此已成为该领域的重点研究对象。
目前,从双相机***重构出原始的三维高光谱图像的方法主要有两类:基于全变差最小化的两步迭代收缩/阈值算法(Two-Step Iterative Shrinkage/ThresholdingAlgorithms,TwIST)和基于自适应字典学习的重构算法。TwIST算法(详见Bioucas-DiasJM,Figueiredo MAT.A New TwIST:Two-Step Iterative Shrinkage/ThresholdingAlgorithms for Image Restoration[J].IEEE Transactions on ImageProcessing.2007,16(12):2992-3004.)建立在迭代收缩阈值算法和迭代加权收缩算法的基础上,其核心是使用反向投影函数去噪(详见Chambolle A.An Algorithm for TotalVariation Minimization and Applications[M].Kluwer Academic Publishers,2004.),并利用前两次的结果进行迭代更新,从而完成高光谱图像的重建。但是TwIST算法的计算复杂度较高,重构时间很长;同时,其重构质量也有待提高。
基于自适应字典学习的重构算法利用光谱图像和灰度图像在空间维的相似性,采用字典学习和稀疏表示的方法完成高光谱图像的重构(详见Wang L,Xiong Z,Gao D,etal.High-speed hyperspectral video acquisition with a dual-camera architecture[J].2015:4942-4950.)。这种方法虽然能够获得较好的重构质量,但是每次重构都需要重新学习字典,导致重构时间比TwIST算法更长,因此也不能满足快速重建的要求,极大地限制了高光谱图像的应用和发展。
针对现有重构算法存在的重构质量低、速度慢等问题,本实施例将高光谱图像重构问题转化为多个子优化问题,采用交替更新、重复迭代的方式完成高光谱图像的重构。同时,本实施例利用NVIDIA公司提出的通用并行计算架构(Compute Unified DeviceArchitecture,CUDA)来完成重构过程,极大地提高高光谱图像重构的效率。本实施例流程图如图2所示。
本实施例公开的一种基于GPU的双相机光谱成像***的快速重构方法,包含以下步骤:
步骤101:输入双相机光谱成像仪的采样图像Y、标定后的前向响应矩阵H、主函数最大迭代次数Maxiters和正则化系数τ。
步骤101中所述双相机光谱成像仪为基于编码孔径快照光谱成像仪(CodedAperture Snapshot Spectral Imager,CASSI)+灰度相机的双通道采样光谱成像仪。双相机光谱成像仪主要由分光镜、物镜、编码模版、中继镜、色散棱镜和灰度相机等部件构成。目标场景的高光谱图像X大小为M×N×Ω,高光谱图像X上任意一点的像素值为x(i,j,λ),1≤i≤M,1≤j≤N,1≤λ≤Ω。其中,M×N表示高光谱图像的空间分辨率,Ω表示高光谱图像的频段数,入射光首先到达分光镜将其能量一分为二,一半进入编码孔径快照光谱成像仪CASSI分支,一半进入灰度相机分支。进入编码孔径快照光谱成像仪CASSI分支的光会到达编码模版进行0-1编码。得到编码后的图像到达色散棱镜后,不同频段的图像会沿着竖直方向进行偏移。最后所有频段的图像到达CCD灰度相机后进行叠加,得到压缩的二维混叠光谱图像。CASSI成像仪的数学模型为:
其中ω(λ)表示CCD相机的频谱响应函数。Cu(i,j)表示编码模版函数。φ(λ)表示色散棱镜的波段偏移函数。yc(i,j)为二维压缩光谱采样图像,v1(i,j)表示高斯白噪声。
进入灰度相机分支的入射光会直接到达灰度相机得到目标场景的二维灰度投影,数学模型为:
将式(1)和式(2)联立,并写成矩阵的形式为:
Y=HX+V (2)
其中Y=[Yc;Yp]表示二维压缩光谱采样图像+二维灰度图像。X表示三维数据立方体,V=[Vc;Vp]表示高斯白噪声。H表示对CASSI成像仪+灰度相机标定后的前向响应矩阵,为编码模版函数、CCD相机频谱响应和色散棱镜偏移函数的综合作用。
步骤102:使用malloc函数申请主机CPU内存,将前向矩阵H和二维压缩光谱采样图像Y读入到主机CPU内存。
步骤102所述malloc函数为C语言内存分配函数,申请的主机CPU内存为一维线性内存,方便后续处理中将CPU内存读入图形处理器GPU显存的操作。
步骤102所述将矩阵H和Y读入主机CPU内存的方式是根据编码孔径光谱仪的色散棱镜偏移方向确定:如果色散棱镜是在竖直方向偏移,则以二维采样图像的水平方向为基准进行读入,即逐行读入;如果色散棱镜是在水平方向偏移,则以二维采样图像的竖直方向为基准进行读入,即逐列读入。目的是方便在后续处理中CASSI***前向模型的计算,提高计算效率。
步骤103:使用cudaMalloc函数为输入输出、以及后续计算过程所需要的数据申请GPU显存,并使用cudaMemcpy函数将二维压缩光谱图像Y和前向响应矩阵H拷贝到GPU显存中。设置通用并行计算架构CUDA计算时需用的线程格/块数。
步骤103所述线程格/块数的设置方法为:线程格/块均为二维布局,单位线程块中的线程个数为(16,16),单位线程格中线程块的个数则是根据高光谱图像的空间分辨率确定为((M+15)/16,(N+15)/16)。
步骤104:初始化重构高光谱图像X0,优化目标函数的初始值f0,以及辅助矩阵S和U,辅助因子ρ。初始化当前迭代次数t=0。
步骤104所述重构高光谱图像X0的初始化方式为:
X0=HTY (4)
其中HT表示CASSI***前向响应H的转置,即从二维压缩光谱图像反演到三维数据立方体的过程。
步骤104所述优化目标函数为全局优化目标函数。根据自然图像的平滑特性,将高光谱重构问题转化为基于全变差约束的最优化问题,从而得到全局优化目标函数为:
||DX||1=∑i,j,λ|x(i,j+1,λ)-x(i,j,λ)|+|x(i+1,j,λ)-x(i,j,λ)| (6)
矩阵D为差分矩阵,它和图像X的作用结果如下:
公式(7)的结果为两个大小和X相同的矩阵,分别表示图像X在水平方向和竖直方向上的差分值。
步骤104所述辅助矩阵S和U是指求解优化方程(4)时所需要的辅助变量,大小均为2×M×N×Ω,且都是用cudaMemset函数将其初始化为全零矩阵。公式(5)无法直接求解,将公式(5)转化为多个子优化问题的求解,原理如下:
引入辅助矩阵S=DX,得到如下带约束的优化方程:
其增广拉格朗日方程为:
其中ρ为辅助因子,本实施例取2.048。当式(9)最小时,变量X、S的值即为优化问题式(8)的解,因此得到两个子优化问题,分别为:
接下来,对X、S和U进行交替更新、循环迭代即可完成高光谱图像的重构。
步骤105:使用cuBLAS和共轭梯度下降算法更新高光谱图像Xt+1。
根据优化问题式(10),X的最小二乘解为:
Xt+1=(HTH+ρDTD)-1(HTY+DT(Ut+ρSt)) (12)
由于矩阵H和矩阵D的规模很大,无法直接求其解析解,因此需利用共轭梯度下降法求高光谱图像Xt+1的近似解(详见Hestenes M R,Stiefel E.Methods of ConjugateGradients for Solving Linear Systems[J].Journal of Research of the NationalBureau of Standards,1952,49(6):409-436.)。同时,为了加速程序的执行效率,本实施例借助NVIDIA提供的矩阵线性运算库cuBLAS来实现共轭梯度下降法,从而完成高光谱图像Xt +1的更新。
步骤106:更新辅助矩阵S和U。
根据优化问题式(11),S的最小二乘解为:
辅助矩阵U的更新公式如下:
Ut+1=Ut+ρ(St-DXt) (14)
使用cuBLAS库,Ut+1能够快速求得。
步骤107:使用更新后的高光谱图像Xt+1计算全局优化目标函数值ft+1。
步骤107所述全局优化目标函数与步骤103的相同,即:
步骤108:根据步骤107计算的结果执行迭代选择策略,完成双相机光谱成像***的快速重构。
步骤108所述迭代选择策略如下:如果收敛,即步骤107的目标函数值ft+1小于上一个目标函数值ft,则计算目标函数值的相对变化量:
若不满足迭代停止条件,即Tol大于预设阈值或者当前迭代次数t小于最大迭代次数Maxiters,则更新辅助因子ρ=1.05×ρ,当前迭代次数t=t+1,目标函数值ft=ft+1,并转至步骤105进行迭代;否则,停止迭代并输出最后一次更新的高光谱图像,从而完成双相机光谱成像***的快速重构。
如果不收敛,即步骤107的目标函数值ft+1大于上一个目标函数值ft,则更新辅助因子ρ=2.0×ρ。判断更新后的ρ值,若ρ小于预设阈值,则转至步骤105进行迭代;否则,停止迭代并输出最后一次更新的高光谱图像,从而完成双相机光谱成像***的快速重构。
为说明本发明的效果,本实施例将在实验条件相同的情况下对两种方法进行对比。
1.实验条件
本实验的硬件测试条件为:Inter i5 3210,内存8G,Matlab 2015b。GPU为GT630M,显存2G,CUDA 7.0。测试所用高光谱图片来自于哥伦比亚大学的Cave数据集。输入的CASSI分支压缩光谱采样图像大小为512×542;灰度相机分支的灰度采样图像大小为512x512;重构以后得到的高光谱图像大小为512×512×31。编码孔径模版为贝努利随机矩阵。主循环最大迭代次数为10次。辅助因子ρ=2.048。共轭梯度下降算法的迭代次数为25次,迭代停止门限值为1.0×10-4。对比算法为基于全变差约束的TwIST算法,其最大迭代次数为200次。两种算法使用的正则化系数均为τ=0.1,预设迭代收敛阈值均为1.0×10-6。
2.实验结果
为了验证本发明的可行性,测试在不同噪声方差sigma下,本发明公开的算法和TwIST算法的重构性能。
首先是两种算法的重构PSNR,单位为dB,结果如表1所示:
表1 PSNR/dB
从表1的结果可以看出,本发明公开的算法迭代10次就能够达到非常好的重构效果,在不同强度的噪声影响下其PSNR均高于TwIST算法。
图3-a、图4-a和图5-a分别为测试图片beads使用本发明仿真重构后,波长分别为440nm、540nm和640nm时的结果;图3-b、图4-b和图5-b分别为使用TwIST算法仿真重构后对应波长的结果。可以看出,本发明公开的方法重构的结果更加清楚,在视觉上的效果要优于TwIST算法。
为了对比两种算法的光谱保真度,其均方根误差RMSE结果如表3所示:
表2 RMSE
从表2的结果可以看出,本发明公开的方法得到的均方根误差更小,表明其在频谱维上的保真度更高,其谱线更接近于真实材料的谱线。
接下来比较两种算法的重构时间,结果如表3所示:
表3重构时间/秒
从表3的结果可以看出,本发明公开的方法重构所用时间非常短,重建时间只需要13秒左右,而TwIST算法重建时间长达1833秒,因此重建速度得到非常大的提升。
此外,为了进一步说明本发明公开的方法在重建效率上的提升,本实施例测试了三种不同高光谱图像分辨率下,两种算法的平均重构时间和加速比,结果如下:
表4不同分辨率下平均重构时间对比
从表4可以看出,在不同分辨率下,本发明公开的方法都能够获得非常高的加速比。同时,随着分辨率的提升,重构计算量成倍增大,本发明公开的方法的加速比也在增加。以分辨率为1340×1018×33的重构图像为例,TwIST算法重建的时间长达2小时多,而使用本发明公开的方法只需要1分钟,其加速效果非常显著。更为重要的,本实施例使用的GPU为GT 630M,属于入门级显卡,如果使用更好的显卡将能够获得更高的加速比。
以上所述的具体描述,对发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (6)
1.一种基于GPU的双相机光谱成像***的快速重构方法,其特征在于:包括以下步骤,
步骤101:输入双相机光谱成像仪的采样图像Y、标定后的前向响应矩阵H、主函数最大迭代次数Maxiters和正则化系数τ;
步骤101中所述双相机光谱成像仪为基于编码孔径快照光谱成像仪(Coded ApertureSnapshot Spectral Imager,CASSI)+灰度相机的双相机采样光谱成像仪;双相机光谱成像仪主要由分光镜、物镜、编码模版、中继镜、色散棱镜和灰度相机部件构成;目标场景的高光谱图像X大小为M×N×Ω,高光谱图像X上任意一点的像素值为x(i,j,λ),1≤i≤M,1≤j≤N,1≤λ≤Ω;其中,M×N表示高光谱图像的空间分辨率,Ω表示高光谱图像的谱段数;入射光首先到达分光镜将其能量一分为二,一半进入编码孔径快照光谱成像仪CASSI分支,一半进入灰度相机分支;进入编码孔径快照光谱成像仪CASSI分支的光会到达编码模版进行0-1编码;经编码后的图像到达色散棱镜后,不同频段的图像会沿着竖直方向进行偏移;最后所有频段的图像到达灰度相机后进行叠加,得到压缩的二维混叠光谱图像;编码孔径快照光谱成像仪CASSI成像仪的数学模型为:
其中ω(λ)表示CCD相机的频谱响应函数,Cu(i,j)表示编码模版函数,φ(λ)表示色散棱镜的波段偏移函数,yc(i,j)为二维压缩光谱采样图像,v1(i,j)表示高斯白噪声;
进入灰度相机分支的入射光会直接到达灰度相机得到目标场景的二维灰度投影,数学模型为:
将式(1)和式(2)联立,并写成矩阵的形式为:
Y=HX+V (3)
其中Y=[Yc;Yp]表示二维压缩光谱采样图像+二维灰度图像;X表示高光谱图像,V=[Vc;Vp]表示高斯白噪声;H表示对编码孔径快照光谱成像仪CASSI+灰度相机标定后的前向响应矩阵,为编码模版函数、色散棱镜偏移函数和CCD灰度相机频谱响应的综合作用;
步骤102:申请主机CPU内存,将前向响应矩阵H和采样图像Y读入到主机CPU内存;
步骤103:为输入输出、以及后续计算过程所需要的数据申请GPU显存,将采样图像Y和前向响应矩阵H拷贝到GPU显存中;设置CUDA并行计算时需用的线程格/块数;
步骤104:初始化重构高光谱图像X0,优化目标函数的初始值f0,以及辅助矩阵S和U,辅助因子ρ;初始化当前迭代次数t=0;
步骤104所述重构高光谱图像X0的初始化方式为:
X0=HTY (4)
其中HT表示编码孔径快照光谱成像仪CASSI***前向响应矩阵H的转置,即从二维压缩光谱图像+二维灰度图像反演到三维数据立方体的过程;
步骤104所述优化目标函数为全局优化目标函数;根据自然图像的平滑特性,将高光谱重构问题转化为基于全变差约束的最优化问题,从而得到全局优化目标函数为:
||DX||1=∑i,j,λ|x(i,j+1,λ)-x(i,j,λ)|+|x(i+1,j,λ)-x(i,j,λ)| (6)
矩阵D为差分矩阵,矩阵D和图像X的作用结果如下:
公式(7)的结果为两个大小和X相同的矩阵,分别表示图像X在水平方向和竖直方向上的差分值;
步骤104所述辅助矩阵S和U是指求解优化方程(4)时所需要的辅助变量,大小均为2×M×N×Ω,初始化为全零矩阵;公式(5)无法直接求解,将公式(5)转化为多个子优化问题的求解,原理如下:
引入辅助矩阵S=DX,得到如下带约束的优化方程:
其增广拉格朗日方程为:
其中ρ为辅助因子,优选2.048;当式(9)最小时,变量X、S的值即为优化问题式(8)的解,因此得到两个子优化问题,分别为:
接下来,对X、S和U进行交替更新、循环迭代即可完成高光谱图像的重构;
步骤105:使用共轭梯度下降算法更新高光谱图像Xt+1;
针对优化问题式(10),X的最小二乘解为:
Xt+1=(HTH+ρDTD)-1(HTY+DT(Ut+ρSt)) (12)
由于矩阵H和矩阵D的规模很大,无法直接求其解析解,因此需利用共轭梯度下降法求高光谱图像Xt+1的近似解,从而完成高光谱图像Xt+1的更新;
步骤106:更新辅助矩阵S和U;
根据优化问题式(11),S的最小二乘解为:
辅助矩阵U的更新公式如下:
Ut+1=Ut+ρ(St-DXt) (14)
步骤107:使用更新后的高光谱图像Xt+1计算全局优化目标函数值ft+1;
步骤107所述全局优化目标函数与步骤104的相同,即:
步骤108:根据步骤107计算的结果执行迭代选择策略,完成双相机光谱成像***的快速重构。
2.如权利要求1所述的一种基于GPU的双相机光谱成像***的快速重构方法,其特征在于:步骤102所述申请主机CPU内存优选C语言的malloc内存分配函数,申请的主机CPU内存为一维线性内存,方便后续处理中将CPU内存读入图像处理器GPU显存的操作。
3.如权利要求2所述的一种基于GPU的双相机光谱成像***的快速重构方法,其特征在于:步骤102所述将矩阵H和Y读入主机CPU内存的方式是根据编码孔径光谱仪的色散棱镜偏移方向确定:如果色散棱镜是沿竖直方向偏移,则以二维采样图像的水平方向为基准进行读入,即逐行读入;如果色散棱镜是沿水平方向偏移,则以二维采样图像的竖直方向为基准进行读入,即逐列读入;目的是方便在后续处理中编码孔径快照光谱成像仪CASSI***前向模型的计算,提高计算效率。
4.如权利要求3所述的一种基于GPU的双相机光谱成像***的快速重构方法,其特征在于:步骤103所述申请图形处理器GPU显存的方法选通用并行计算架构CUDA的cudaMalloc函数,拷贝到图形处理器GPU显存的方法选通用并行计算架构CUDA的cudaMemcpy函数。
5.如权利要求4所述的一种基于GPU的双相机光谱成像***的快速重构方法,其特征在于:步骤106中Ut+1的更新方式选矩阵线性运算库cuBLAS库,能够提高程序效率。
6.如权利要求1所述的一种基于GPU的双相机光谱成像***的快速重构方法,其特征在于:步骤108所述迭代选择策略如下:如果收敛,即步骤107的目标函数值ft+1小于上一个目标函数值ft,则计算目标函数值的相对变化量:
若不满足迭代停止条件,即Tol大于预设阈值或者当前迭代次数t小于最大迭代次数Maxiters,则更新辅助因子ρ=1.05×ρ,当前迭代次数t=t+1,目标函数值ft=ft+1,并转至步骤105进行迭代;否则,停止迭代并输出最后一次更新的高光谱图像,从而完成双相机光谱成像***的快速重构;
如果不收敛,即步骤107的目标函数值ft+1大于上一个目标函数值ft,则更新辅助因子ρ=2.0×ρ;判断更新后的ρ值,若ρ小于预设阈值,则转至步骤105进行迭代;否则,停止迭代并输出最后一次更新的高光谱图像,从而完成双相机光谱成像***的快速重构。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710699197.6A CN107525588B (zh) | 2017-08-16 | 2017-08-16 | 一种基于gpu的双相机光谱成像***的快速重构方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710699197.6A CN107525588B (zh) | 2017-08-16 | 2017-08-16 | 一种基于gpu的双相机光谱成像***的快速重构方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107525588A CN107525588A (zh) | 2017-12-29 |
CN107525588B true CN107525588B (zh) | 2020-04-17 |
Family
ID=60681287
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710699197.6A Expired - Fee Related CN107525588B (zh) | 2017-08-16 | 2017-08-16 | 一种基于gpu的双相机光谱成像***的快速重构方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107525588B (zh) |
Families Citing this family (14)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108805975B (zh) * | 2018-05-29 | 2021-03-16 | 常熟理工学院 | 一种基于改进迭代收缩阈值算法的微观3d重建方法 |
CN109146787B (zh) * | 2018-08-15 | 2022-09-06 | 北京理工大学 | 一种基于插值的双相机光谱成像***的实时重建方法 |
CN109358379B (zh) * | 2018-10-30 | 2020-04-21 | 西安石油大学 | 修正总变分模型约束下基于泛函重构的地球物理反演方法 |
CN109447890A (zh) * | 2019-01-09 | 2019-03-08 | 北京理工大学 | 一种基于卷积神经网络的光谱成像***的编码优化方法 |
CN109741407A (zh) * | 2019-01-09 | 2019-05-10 | 北京理工大学 | 一种基于卷积神经网络的光谱成像***的高质量重构方法 |
CN109886898B (zh) * | 2019-03-05 | 2020-10-02 | 北京理工大学 | 基于优化启发的神经网络的光谱成像***的成像方法 |
CN110717947B (zh) * | 2019-09-25 | 2021-04-27 | 北京理工大学 | 一种基于外部和内部训练的高质量光谱重构方法 |
CN110926611A (zh) * | 2019-12-12 | 2020-03-27 | 北京理工大学 | 一种应用于压缩感知光谱成像***的噪声抑制方法 |
CN111353584B (zh) * | 2020-02-20 | 2023-04-07 | 中山大学 | 一种基于时间序列分析的深度学习训练任务行为预测方法 |
CN111397733B (zh) * | 2020-04-23 | 2021-03-02 | 湖南大学 | 一种单/多帧快照式光谱成像方法、***及介质 |
CN111609931A (zh) * | 2020-05-17 | 2020-09-01 | 北京安洲科技有限公司 | 一种参数可编程的实时高光谱采集***和方法 |
CN112556848B (zh) * | 2020-11-28 | 2022-09-06 | 南京理工大学 | 双相机压缩测量协同张量表示的融合计算成像方法 |
CN113984206A (zh) * | 2021-10-21 | 2022-01-28 | 北京航天创智科技有限公司 | 一种基于多台阶衍射滤光器的快照式光谱成像*** |
CN114485942B (zh) * | 2022-02-16 | 2024-05-28 | 南京大学 | 一种高光谱配准方法及其成像*** |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103400350A (zh) * | 2013-07-22 | 2013-11-20 | 中国科学院西安光学精密机械研究所 | 一种编码孔径光谱成像仪的光谱图像复原方法 |
CN105139352A (zh) * | 2015-08-13 | 2015-12-09 | 西安电子科技大学 | 基于gpu的混叠采样的光谱视频快速重构方法 |
-
2017
- 2017-08-16 CN CN201710699197.6A patent/CN107525588B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103400350A (zh) * | 2013-07-22 | 2013-11-20 | 中国科学院西安光学精密机械研究所 | 一种编码孔径光谱成像仪的光谱图像复原方法 |
CN105139352A (zh) * | 2015-08-13 | 2015-12-09 | 西安电子科技大学 | 基于gpu的混叠采样的光谱视频快速重构方法 |
Non-Patent Citations (1)
Title |
---|
Fast Parallel Implementation of Dual-Camera;Shipeng Zhang et.al;《IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY》;20181231;全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN107525588A (zh) | 2017-12-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107525588B (zh) | 一种基于gpu的双相机光谱成像***的快速重构方法 | |
US10985777B2 (en) | Signal recovery via deep convolutional networks | |
CN110501072B (zh) | 一种基于张量低秩约束的快照式光谱成像***的重构方法 | |
CN109146787B (zh) | 一种基于插值的双相机光谱成像***的实时重建方法 | |
EP3716198A1 (en) | Image reconstruction method and device | |
Chen et al. | Hyperspectral image compressive sensing reconstruction using subspace-based nonlocal tensor ring decomposition | |
CN112001914A (zh) | 深度图像补全的方法和装置 | |
CN107451956B (zh) | 一种编码孔径光谱成像***的重构方法 | |
CN109697697B (zh) | 基于优化启发的神经网络的光谱成像***的重构方法 | |
Lai et al. | Deep plug-and-play prior for hyperspectral image restoration | |
Ye et al. | CSformer: Bridging convolution and transformer for compressive sensing | |
CN109886898B (zh) | 基于优化启发的神经网络的光谱成像***的成像方法 | |
CN110148103A (zh) | 基于联合优化的高光谱和多光谱图像融合方法、计算机可读存储介质、电子设备 | |
CN106097278A (zh) | 一种多维信号的稀疏模型、重建方法和字典训练方法 | |
Ramirez et al. | LADMM-Net: An unrolled deep network for spectral image fusion from compressive data | |
CN113962858A (zh) | 一种多视角深度获取方法 | |
Bernabé et al. | Parallel hyperspectral coded aperture for compressive sensing on gpus | |
CN116402679B (zh) | 一种轻量级红外超分辨率自适应重建方法 | |
Han et al. | Deep learning simulations of the microwave sky | |
CN109948462B (zh) | 基于多gpu协同交互数据流组织的高光谱图像快速分类法 | |
Cherian et al. | A Novel AlphaSRGAN for Underwater Image Super Resolution. | |
Polasek et al. | Vision UFormer: Long-range monocular absolute depth estimation | |
Luo et al. | Piecewise linear regression-based single image super-resolution via Hadamard transform | |
Liu et al. | Fine-grained MRI reconstruction using attentive selection generative adversarial networks | |
CN115760670B (zh) | 基于网络隐式先验的无监督高光谱融合方法及装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200417 Termination date: 20200816 |