CN104112294B - 基于稀疏约束的强度关联成像高速三维重构系统及方法 - Google Patents

Abstract

一种基于稀疏约束的强度关联成像高速三维重构系统，包括数据采集单元、信号处理机、高速存储器、GPU并行计算控制器、GPU计算卡组、三维处理机、三维成像显示器。数据采集单元由CCD和PMT构成，GPU计算卡组由多块GPU计算卡并联构成。本发明可实现准实时强度关联成像三维图像重构，具有计算速度快、信噪比高、成本低和功耗小的特点，可应用于可见光及X光强度关联成像。

Description

基于稀疏约束的强度关联成像高速三维重构系统及方法

技术领域

本发明涉及三维信息处理领域，特别是一种基于稀疏约束的强度关联成像高速三维重构系统及方法，它可应用于可见光及X光强度关联成像领域。

背景技术

强度关联成像不同于常规成像，探测器拍摄到的并非物体图像本身，而是记录热光场的涨落信息，再通过后续信息处理来获取目标物体图像。因此，在强度关联成像应用中，一个重要的问题是重构目标物体图像。采用关联算法能快速得到重构图像，但重构图像信噪比低，只适合于噪声影响较低的实验室成像。根据压缩感知(Compressive Sensing)理论，当长度为N的信号X为稀疏且测量矩阵A满足特定条件时，可以通过CS算法由远少于N的M次测量数据恢复出信号。对于噪声影响较大的远距离三维成像实用系统，可以利用CS方法实现强度关联图像重构。

CS重构图像信噪比高，但数据计算量大，所需重构时间较长。对于二维强度关联成像，只需重构单帧图像，重构算法所需时间尚可接受。但对于三维成像，是对二维切片图像分别进行单帧重构，再由二维切片组合得到物体的三维信息，重构时间由图像序列的帧数决定，要实现高分辨率的三维成像，需要对数千帧图像进行重构，重构时间难以满足应用需求。借助大型计算机或者集群虽然能够实现快速并行计算，但所需成本高，且无法实现户外移动计算。

浙江大学的刘旭等提出了一种基于压缩感知的三维物体成像的简单计算方法，只需要两次CS计算就可以获得目标物体的三维信息，但需要反射镜和空间光调制器DMD等，这些器件只对可见光发生作用，对于高穿透性的X光，没有对应器件，无法应用于X光强度关联成像。

基于稀疏约束的强度关联成像(GISC)方法是中科院上海光机所龚文林等提出的，通过关联计算实现了遥感三维成像，但成像信噪比低，图像模糊，而CS重构方法理论上能获得较高信噪比的图像。

强度关联成像进行m次测量，采集得到的数据为：(1)二维参考图像构成的参考序列{Ir₁,…,Ir_i,…Ir_m}，每帧二维参考图像Ir_i的像素数为n，用矩阵A来表示参考序列，一帧二维参考图像构成A的一行，则矩阵A为m行n列；(2)一维信号数据构成的信号序列{It₁,…,It_i,…It_m}，每帧一维信号数据It_i包含k个元素，用矩阵Y来表示信号序列，一帧一维信号数据构成Y的一行，则矩阵Y为m行k列矩阵。

待求三维重构目标为二维重构图像构成的图像序列{I₁,…,I_j,…I_k}，传统CS算法是通过依次求解每帧二维重构图像实现三维重构。在三维并行CS算法中，每帧二维重构图像I_j的像素数为n，图像序列的帧数为k，用矩阵X来表示三维重构目标，一帧二维重构图像构成X的一列，则矩阵X为n行k列矩阵，三维CS重构过程就是由已知的Y和A求解X的过程，可表示为：

AX＝Y，

其中，AX_j＝Y_j，X_j为n元素向量，Y_j为m元素向量，X＝[X₁,…,X_j,…X_k]，Y＝[Y₁,…,Y_j,…Y_k]。

根据CS理论，由(其中φ代表稀疏变换)可知，稀疏约束条件下强度关联成像三维重构就是求解欠定方程组可转化为求解线性规划问题：

发明内容

本发明要解决的技术问题在于克服上述在先技术的缺陷，提供一种基于稀疏约束的强度关联成像高速三维重构系统及方法。该三维重构系统及方法实现准实时强度关联成像三维图像重构，具有计算速度快、信噪比高、成本低和功耗小的特点，可应用于可见光及X光强度关联成像。

本发明的技术解决方案如下：

一种基于稀疏约束的强度关联成像高速三维重构系统，其特征在于其构成包括数据采集单元、信号处理机、高速存储器、GPU并行计算控制器、GPU计算卡组、三维处理机、三维成像显示器；

所述的数据采集单元由CCD和PMT构成，所述的GPU计算卡组由多块GPU计算卡并联构成；所述的信号处理机、GPU并行计算控制器、三维处理机与高速存储器之间采用高速光缆连接，所述的GPU计算卡组与GPU并行计算控制器通过PCI-E总线连接，所述的三维成像显示器与三维处理机相连接；

所述的数据采集单元将采集到的数据输入信号处理机，信号处理机对采集数据进行噪声过滤处理后，将有效数据传递到高速存储器进行存储；

所述的GPU并行计算控制器从高速存储器读取有效数据进行并行计算预处理后，传递到GPU计算卡组的多个GPU计算卡，GPU计算卡组完成并行计算后将计算结果返回到GPU并行计算控制器，GPU并行计算控制器进行并行计算后处理之后将结果数据发送到高速存储器进行存储；

所述的三维处理机从高速存储器读取结果数据后进行三维降噪合成，并通过三维成像显示器完成结果显示。

利用本发明基于稀疏约束的强度关联成像高速三维重构系统进行高速三维重构的方法，包括以下步骤：

<1>获取强度关联成像三维重构数据矩阵：

强度关联成像中，进行m次测量，CCD采集到的数据为二维参考图像构成的参考序列{Ir₁,…,Ir_i,…Ir_m}，PMT采集到的数据为一维信号数据构成的信号序列{It₁,…,It_i,…It_m}，其中i＝1,2,3,…,m；

对于参考序列{Ir₁,…,Ir_i,…Ir_m}，每帧二维参考图像Ir_i的像素数为n，用矩阵A来表示参考序列，一帧二维参考图像构成A的一行，则矩阵A为m行n列；

对于信号序列{It₁,…,It_i,…It_m}，每帧一维信号数据It_i包含k个元素，用矩阵Y来表示信号序列，一帧一维信号数据构成Y的一行，则矩阵Y为m行k列矩阵；

将参考数据矩阵A和信号数据矩阵Y输入信号处理机。

<2>进行噪声过滤处理：

求信号数据矩阵Y每一列的平均值，得到1行k列的向量，将此向量归一化后，得到1行k列的向量Ys，对Ys的k个元素逐个判断，若小于门限阈值Ts则视为噪声，Ts∈(0,1)，剔除矩阵Y对应的列，得到新的矩阵Y为m行K列，其中K<k；

求信号数据矩阵Y每一行的平均值，得到m行1列的向量Ym，将矩阵Y的每一列Y_j(j＝1,2,3,…,K)与Ym做运算|Y_j-Ym|./Ym，得到m行K列的矩阵Ya，对Ya的m×K个元素逐个判断，若小于门限阈值Ta则视为噪声，Ta∈(0,0.5)，矩阵Y对应的元素置0；

对矩阵Y的每一行作判断，剔除矩阵Y的全0行，及矩阵A对应的行，得到有效数据：M行K列的信号数据矩阵Y’和M行n列的参考数据矩阵A’，其中M<m；

将参考数据矩阵A’和信号数据矩阵Y’存储到高速存储器。

<3>进行并行计算预处理：

由GPU并行计算控制器从高速存储器读取参考数据矩阵A’和信号数据矩阵Y’，根据A’和Y’计算参数μ，计算公式为：

其中μ₀∈(0,1)；

对M行n列的参考数据矩阵A’进行数据分割，GPU计算卡的数目为N，则将参考数据矩阵A’分割为N块，每块对应一个GPU计算卡，每个分块矩阵A’_GPU的大小为M/N行n列；

将参数μ和信号数据矩阵Y’传递到每个GPU计算卡，将分块矩阵A’_GPU传递到对应的GPU计算卡。

<4>进行并行计算：

在多个GPU计算卡上同时对分块矩阵A’_GPU进行稀疏变换，其变换公式为：

其中，φ代表稀疏变换，可以是离散余弦变换(DCT)、小波变换等等；

在多个GPU计算卡上同时根据A_new、Y’和μ求解X’,即求解线性规划问题：

具体算法为：

令作为目标函数，X’初值设为0，沿方向搜索，即其中

α₀＝1，循环终止条件为

将每个GPU计算卡求得的X’矩阵(n/N行K列)返回到GPU并行计算控制器。

<5>进行并行计算后处理：

对多个GPU计算卡返回的矩阵X’进行数据合并，得到完整的n行K列的结果数据矩阵X’；

对X’进行反稀疏变换获取三维重构目标的实空间表达矩阵X，其变换公式为：

其中，φ^-1代表逆稀疏变换，可以是逆离散余弦变换(IDCT)、逆小波变换等等；

将结果数据矩阵X(n行K列)存储到高速存储器。

<6>进行三维降噪合成：

利用所述的三维处理机从高速存储器读取结果数据矩阵X(n行K列)，对矩阵X的每一行搜索最大值(此处为信号，其余为噪声)，记录下此位置和此位置的值，得到位置向量P(n行1列)和值向量V(n行1列)，将P向量和V向量归一化后以hsv真彩格式表示，P向量对应色度，V向量对应亮度，饱和度设为1，即可得第三维以颜色表示的三维目标真彩图。

本发明与在先技术相比，主要有以下几方面的优点：

1、适用于多种三维成像。本发明强度关联成像三维重构系统及方法，与成像光源无关，可用于X光强度关联成像，而部分在先技术仅适用于可见光强度关联成像。本发明强度关联成像三维重构系统及方法，与成像光路无关，可用于三维距离强度关联成像，三维多光谱强度关联成像，三维多点探测强度关联成像等等，而部分在先技术仅适用于三维距离强度关联成像。

2、高信噪比。本发明基于稀疏约束的强度关联成像高速三维重构系统及方法，其重构方法是基于稀疏变换和并行CS算法。根据压缩感知理论，信号越是稀疏，重构所需的测量次数越少，在同样的测量次数时，信号越是稀疏，重构图像质量越高。对于每个二维切片1024*1024像素来说，需要100万次测量才能得到全部待求信息，而采用稀疏变换和CS算法后只需要几千次测量就可高质量重构待求图像。因此，本发明理论上具有比传统算法更高的信噪比。

3、计算速度快。本发明基于稀疏约束的强度关联成像高速三维重构系统及方法，从三个方面提高计算速度。第一，对信号数据进行过滤，选择高变化量信号数据及其对应参考数据，从而降低部分计算量，提高计算速度；第二，采用并行CS算法，速度相对于传统CS算法提高10倍左右；第三，采用GPU设备进行数据级并行加速，最终速度相对于传统串行重构方案提高100倍以上。

4、成本低，功耗小，可移动。现有的常规CS并行计算是采用大型计算机或者集群实现，虽然能够快速计算，但所需成本高(>100万)，功耗大，无法实现户外移动计算。本发明基于稀疏约束的强度关联成像高速三维重构系统采用GPU并行计算方式，8块GPU并行计算卡价格约16万，整套系统价格不超过25万，且功耗小，采用移动供电系统即可。

附图说明

图1是本发明基于稀疏约束的强度关联成像高速三维重构系统及方法的实施例结构示意图，图中：

1：数据采集单元，101：CCD，102：PMT，2：信号处理机，3：高速存储器，4：GPU并行计算控制器，5：GPU计算卡组，6：三维处理机，7：三维成像显示器。

图2是本发明基于稀疏约束的强度关联成像高速三维重构系统及方法的实施例三维重构数据模型示意图

具体实施方式

本发明是一种基于稀疏约束的强度关联成像高速三维重构系统，其构成包括数据采集单元1、信号处理机2、高速存储器3、GPU并行计算控制器4、GPU计算卡组5、三维处理机6、三维成像显示器7。

数据采集单元1由CCD101和PMT102构成，GPU计算卡组5由多块GPU计算卡并联构成；信号处理机2、GPU并行计算控制器4、三维处理机6与高速存储器3之间采用高速光缆连接，GPU计算卡组5与GPU并行计算控制器4通过PCI-E总线连接，三维成像显示器7与三维处理机6相连接；

数据采集单元1将采集到强度关联成像数据Y和A输入信号处理机2，信号处理机2对采集数据进行噪声过滤处理后，将有效数据Y’和A’传递到高速存储器3进行存储；

GPU并行计算控制器4从高速存储器2读取有效数据Y’和A’进行并行计算预处理后得到参数μ和分块矩阵A’_GPU，将μ与Y’传递到GPU计算卡组5的每个GPU计算卡，将分块矩阵A’_GPU传递到对应的GPU计算卡，GPU计算卡组5完成并行计算后将计算结果X’返回到GPU并行计算控制器4，GPU并行计算控制器4进行并行计算后处理之后将结果数据X发送到高速存储器3进行存储；

三维处理机6从高速存储器3读取结果数据X后进行三维降噪合成，并通过三维成像显示器7完成结果显示。

本实施例中，GPU计算卡组为8块GPU计算卡并联，三维强度关联遥感成像采样次数m＝5000，待测目标图像像素数为n＝256*256＝65536，信号数据长度为k＝1024，成像距离为800米至1000米，其三维CS重构数据模型如图2所示。完成三维并行CS重构时间约1分钟，而传统CS重构依次求解每帧二维重构图像，完成重构计算需要4小时以上。

本实施例基于稀疏约束的强度关联成像高速三维重构系统进行高速三维重构的方法，包括以下步骤：

<1>获取强度关联成像三维重构数据矩阵：

进行5000次测量，CCD采集到的数据为二维参考图像构成的参考序列{Ir₁,…,Ir_i,…Ir_m}，PMT采集到的数据为一维信号数据构成的信号序列{It₁,…,It_i,…It_m}，其中i＝1,2,3,…,5000；

对于参考序列{Ir₁,…,Ir_i,…Ir_m}，每帧二维参考图像Ir_i的像素数为65536，用矩阵A来表示参考序列，一帧二维参考图像构成A的一行，则矩阵A为5000行65536列；

对于信号序列{It₁,…,It_i,…It_m}，每帧一维信号数据It_i包含1024个元素，用矩阵Y来表示信号序列，一帧一维信号数据构成Y的一行，则矩阵Y为5000行1024列矩阵；

将参考数据矩阵A和信号数据矩阵Y输入信号处理机2。

<2>进行噪声过滤处理：

求信号数据矩阵Y每一列的平均值，得到1行1024列的向量，将此向量归一化后，得到1行1024列的向量Ys，对Ys的1024个元素逐个判断，若小于门限阈值Ts则视为噪声(Ts＝0.05)，剔除矩阵Y对应的列，得到新的矩阵Y为5000行963列；

求信号数据矩阵Y每一行的平均值，得到5000行1列的向量Ym，将矩阵Y的每一列Y_j(j＝1,2,3,…,963)与Ym做运算|Y_j-Ym|./Ym，得到5000行963列的矩阵Ya，对Ya的5000×963个元素逐个判断，若小于门限阈值Ta则视为噪声(Ta＝0.1)，矩阵Y对应的元素置0；

对矩阵Y的每一行作判断，剔除矩阵Y的全0行(共剔除976行)，同时剔除矩阵A对应的行，得到有效数据：4024行963列的信号数据矩阵Y’和4024行65536列的参考数据矩阵A’；

将参考数据矩阵A’和信号数据矩阵Y’存储到高速存储器3。

<3>进行并行计算预处理：

由GPU并行计算控制器从高速存储器读取参考数据矩阵A’和信号数据矩阵Y’，根据A’和Y’计算参数μ，计算公式为：

其中μ₀＝0.5；

对4024行65536列的参考数据矩阵A’进行数据分割，GPU计算卡的数目为8，则将参考数据矩阵A’分割为8块，每块对应一个GPU计算卡，每个分块矩阵A’_GPU的大小为503行65536列；

将参数μ和信号数据矩阵Y’传递到每个GPU计算卡，将分块矩阵A’_GPU传递到对应的GPU计算卡。

<4>进行并行计算：

在8个GPU计算卡上同时对分块矩阵A’_GPU进行稀疏变换，其变换公式为：

其中，φ代表稀疏变换(离散余弦变换)；

在8个GPU计算卡上同时根据A_new、Y’和μ求解X’,即求解线性规划问题：

具体算法为：

令作为目标函数，X’初值设为0，沿方向搜索，即

其中α₀＝1，循环终止条件为

将8个GPU计算卡求得的X’矩阵块(8192行963列)返回到GPU并行计算控制器4。

<5>进行并行计算后处理：

对8个GPU计算卡返回的X’矩阵块进行数据合并，得到完整的65536行963列的结果数据矩阵X’；

对X’进行反稀疏变换获取三维重构目标的实空间表达矩阵X，其变换公式为：

其中，φ^-1代表逆稀疏变换(逆离散余弦变换)；

将结果数据矩阵X(65536行963列)存储到高速存储器3。

<6>进行三维降噪合成：

从高速存储器3读取结果数据矩阵X(65536行963列)，对矩阵X的每一行搜索最大值(此处为信号，其余为噪声)，记录下此位置和此位置的值，得到位置向量P(65536行1列)和值向量V(65536行1列)，将P向量和V向量归一化后以hsv真彩格式表示，P向量对应色度，V向量对应亮度，饱和度设为1，即可得第三维以颜色表示的三维目标真彩图。

实验表明，本发明基于稀疏约束的强度关联成像高速三维重构系统及方法，可实现准实时高速三维图像重构，具有计算速度快、信噪比高的特点，成本低，功耗小，可广泛应用于强度关联成像领域。

Claims (1)

1.利用基于稀疏约束的强度关联成像高速三维重构系统进行高速三维重构的方法，该系统构成包括数据采集单元、信号处理机、高速存储器、GPU并行计算控制器、GPU计算卡组、三维处理机、三维成像显示器；

所述的数据采集单元由CCD和光电倍增管PMT构成，所述的GPU计算卡组由多块GPU计算卡并联构成；所述的信号处理机、GPU并行计算控制器、三维处理机与高速存储器之间采用高速光缆连接，所述的GPU计算卡组与GPU并行计算控制器通过PCI-E总线连接，所述的三维成像显示器与三维处理机相连接；所述的数据采集单元将采集到的数据输入信号处理机，信号处理机对采集数据进行噪声过滤处理后，将有效数据传递到高速存储器进行存储；所述的GPU并行计算控制器从高速存储器读取有效数据进行并行计算预处理后，传递到GPU计算卡组的多个GPU计算卡，GPU计算卡组完成并行计算后将计算结果返回到GPU并行计算控制器，GPU并行计算控制器进行并行计算后处理之后将结果数据发送到高速存储器进行存储；所述的三维处理机从高速存储器读取结果数据后进行三维降噪合成，并通过三维成像显示器完成结果显示；其特征在于，该方法包括以下步骤：

<1>获取强度关联成像三维重构数据矩阵：

强度关联成像中，进行m次测量，CCD采集到的数据为二维参考图像构成的参考序列{Ir₁,…,Ir_i,…Ir_m}，光电倍增管PMT采集到的数据为一维信号数据构成的信号序列{It₁,…,It_i,…It_m}，其中i＝1,2,3,…,m；

对于参考序列{Ir₁,…,Ir_i,…Ir_m}，每帧二维参考图像Ir_i的像素数为n，用矩阵A来表示参考序列，一帧二维参考图像构成A的一行，则矩阵A为m行n列；

对于信号序列{It₁,…,It_i,…It_m}，每帧一维信号数据It_i包含k个元素，用矩阵Y来表示信号序列，一帧一维信号数据构成矩阵Y的一行，则矩阵Y为m行k列矩阵；

将参考数据矩阵A和信号数据矩阵Y输入信号处理机；

<2>进行噪声过滤处理：

求信号数据矩阵Y每一列的平均值，得到1行k列的向量，将此向量归一化后，得到1行k列的向量Ys，对Ys的k个元素逐个判断，若小于门限阈值Ts则视为噪声，Ts∈(0,1)，剔除矩阵Y对应的列，得到新的矩阵Y_new为m行K列，其中K<k；

求矩阵Y_new每一行的平均值，得到m行1列的向量Ym，将矩阵Y_new的每一列Y_j(j＝1,2,3,…,K)与Ym做运算|Y_j-Ym|./Ym，得到m行K列的矩阵Ya，对Ya的m×K个元素逐个判断，若小于门限阈值Ta则视为噪声，Ta∈(0,0.5)，矩阵Y_new对应的元素置0；

对矩阵Y_new的每一行作判断，剔除矩阵Y_new的全0行和矩阵A对应的行，得到有效数据：M行K列的信号数据矩阵Y’和M行n列的参考数据矩阵A’，其中M<m；

将参考数据矩阵A’和信号数据矩阵Y’存储到高速存储器；

<3>并行计算预处理：

由GPU并行计算控制器从高速存储器读取参考数据矩阵A’和信号数据矩阵Y’，根据A’和Y’按下列公式计算参数μ：

其中μ₀∈(0,1)；

对M行n列的参考数据矩阵A’进行数据分割，GPU计算卡的数目为N，则将参考数据矩阵A’分割为N块，每块对应一个GPU计算卡，每个分块矩阵A’_GPU的大小为M/N行n列；

将参数μ和信号数据矩阵Y’传递到每个GPU计算卡，将分块矩阵A’_GPU传递到对应的GPU计算卡；

<4>进行并行计算：

在多个GPU计算卡上同时对分块矩阵A’_GPU进行稀疏变换，其变换公式为：

其中，代表稀疏变换；

在多个GPU计算卡上同时根据A_new、Y’和μ求解X’：

$X^{\&prime;}=\min \frac{1}{2}\left|\right|Y^{\&prime;}-A_{new}{X}^{\&prime;}|{|}_{l2}^2+\mathrm{\&mu;}|\left|X^{\&prime;}\right|{|}_{l1},$

具体算法为：

令作为目标函数，X’初值设为0，沿方向搜索，即其中α₀＝1，循环终止条件为：

将每个GPU计算卡求得的n/N行K列的矩阵X’返回到GPU并行计算控制器；

<5>进行并行计算后处理：

对多个GPU计算卡返回的矩阵X’进行数据合并，得到完整的n行K列的结果数据矩阵X’；

对完整的n行K列的结果数据矩阵X’进行反稀疏变换获取三维重构目标的实空间表达矩阵X，其变换公式为：

其中，代表逆稀疏变换；

将结果数据n行K列的矩阵X存储到高速存储器；

<6>进行三维降噪合成：

利用所述的三维处理机从高速存储器读取结果数据矩阵X，对矩阵X的每一行搜索最大值，记录下此位置和此位置的值，得到n行1列的位置向量P和n行1列的值向量V，将P向量和V向量归一化后以hsv真彩格式表示，P向量对应色度，V向量对应亮度，饱和度设为1，即得第三维以颜色表示的三维目标真彩图。