CN106097291A - 一种基于径向梯度全变差的art火焰切片重构算法 - Google Patents

Abstract

本发明公开了一种基于径向梯度全变差的ART火焰切片重构算法，对火焰图像进行预处理；生成投影图像；根据ART准则求迭代权重矩阵W；根据ART代数迭代法重构图像；求解每次迭代图像的径向梯度和径向梯度全变差；根据径向TV方法调整图像；输出重构图像。与现有技术相比，本发明能够有效提高火焰图像的切片重构质量和重构速度；尤其在投影角度较少的情况下，该算法的去相关性更好，径向梯度图像更稀疏，且火焰重构过程稳定性好，算法鲁棒性强。

Description

一种基于径向梯度全变差的ART火焰切片重构算法

技术领域

本发明涉及火焰图像重构领域，特别是涉及一种内燃机燃烧火焰切片重构算法。

背景技术

近年来汽车尾气排放造成的空气污染日趋严重，已成为雾霾天气形成的重要原因之一，内燃机的节能减排已引起全球关注。燃烧火焰作为内燃机缸内燃烧过程的重要表征，对研究火焰燃烧机理、降低燃烧排污和提高内燃机性能具有重要意义。现代光学方法和计算机技术的发展，使得利用光学可视化技术直接对内燃机火焰燃烧过程进行观测成为可能。目前对内燃机燃烧火焰研究，大多利用从活塞底部或顶部获取的缸内火焰单幅二维投影图像进行二维特征分析，但由于火焰具有半透明性，利用该方法拍摄的火焰图像是整个缸内燃烧火焰在特定投影面的二维叠加结果，不能表征火焰的三维空间特征。因此有必要对火焰进行三维形态重构，为火焰的定量分析提供足够信息。

国内外科研人员对火焰三维形态重构进行了相应研究，主要利用高速同步摄像机拍摄火焰在多个方向的二维投影图像，并基于计算机断层扫描(Computed Tomography，简称CT)技术对火焰二维切片进行重构以形成三维立体图像。多伦多大学的Samuel借鉴层析成像技术和CT技术图像重构的思想，采用分层重构的方法对自然界火焰进行三维重构。J.Floyd等基于5个CCD相机及镜像***对火焰进行拍摄，并利用CT技术对获取的10幅火焰视图进行重构，其重构切片效果比较粗糙，有待于进一步完善。日本的Tadashi Ito等引入一个外部辐射源的开关来控制火焰热辐射投影的获取并使用CT技术来重构火焰的温度分布图像。中国科技大学的周怀春和浙江大学的岑可法院士等对炉膛内火焰三维重构方面进行了大量研究，采用数字化摄像装置从炉膛燃烧空间提取燃烧二维火焰辐射图像，并利用CT技术进行火焰三维温度场的重构。因此，对火焰三维形态重构的关键是对火焰二维投影图像进行切片重构，其切片重构精度决定了火焰三维重构的效果，但上述研究方法主要利用CT技术直接对火焰切片进行重构，没有有效利用火焰的稀疏特性，导致其在投影角度较小的情况下重构效果较差。

发明内容

针对以上问题，本发明提出了一种基于径向TV的ART火焰切片重构算法，该方法是对火焰二维切片重构算法的优化，依据燃烧火焰的径向扩散特征，考虑到火焰图像在全变差(Total Variation，TV)域的稀疏性，将传统TV从笛卡尔坐标系扩展到极坐标系模型中，引入到代数重构法(Algebraic reconstruction technique,简称ART)火焰二维切片重构中。

本发明提出了基于径向TV的ART火焰切片重构算法，该方法包括以下步骤：

步骤1、对火焰切片图像进行预处理：包括灰度变换和平滑去噪操作；

步骤2、求得火焰切片图像重建的迭代权重矩阵W并获得相应的投影图像矩阵P，具体步骤如下：

将二维火焰切片离散为N＝n×n个独立的体元，每个体元都具有特定的代表该单元亮度大小的值x_j(j＝1,2,…,N)，二维投影图像上的每个像素点都对应一条穿过该切片体元的射线L_k(k＝1,2,…,M)，M是所有投影的射线总数，等于投影角度个数与某角度投影射线数的乘积；而该像素点的值P_k(k＝1,2,…,M)为这条射线上所有体元对该射线亮度贡献的叠加；设ω_jk为体元j对射线L_k的权重因子，代表了体元与射线间的相关性，体元j对射线L_k的贡献为：

P_jk＝ω_jk·x_j

切片上所有体元对L_k亮度贡献和表示为：

$P_k=\underset{j=1}{\overset{N}{\Σ}}{P}_{jk}=\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}\·x_j,k=1,2,...,M$

用矩阵的形式表示如下：

Wx＝P

P为M维投影矩阵，表示二维投影的像素值；x为N维待重构切片数据矢量；W为M×N维投影权重矩阵；

步骤(3)、设定初始迭代次数i＝1和迭代次数的最大值Maxcount＝15，并设置迭代的初始向量f⁽⁰⁾＝{0,0,...,0}；

步骤4、计算重构图像，方程组公式如下：

$f^{\left(m\right)}=f^{(m-1)}+{\mathrm{\ω}}_{jk}(P_k-\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}{f}^{(m-1)})/\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}^2$

其中m＝1,2,…,M代表每一次完整迭代的方程组的方程个数；根据上述方程组公式，计算M×N方程组；将初始向量f⁽⁰⁾＝{0,0,...,0}代入方程组的第一个方程求出f⁽¹⁾，由第二个方程求出f⁽²⁾，依次类推，最终由第M个方程求得f^(M)，这样就完成了一次完整迭代；

步骤5、求解f^(M)的径向梯度公式如下：

$\▿f^{\left(M\right)}=(D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}},D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}})$

其中D_ρ和D_θ分别是沿极径和极角方向的离散微分算子，有如下公式计算：

$D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{-},\mathrm{\θ}}$

$D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ},\mathrm{\θ}}_{-}}$

其中ρ_-和ρ₊分别代表沿极径减少和极径增加方向的变化，θ_-和θ₊则分别代表沿极角顺时针和极角逆时针方向的变化；

步骤6、求解f^(M)径向梯度全变差的偏导函数并对f^(M)进行径向梯度全变差优化，如下式：

$f^{\left(M\right)}\′=f^{\left(M\right)}-q\·\▿{||f^{\left(M\right)}||}_{R-TV}$

其中，q为下降步长因子，本发明中取值为0.1/k，并用f^(M)'代换原来的f^(M)；

步骤7、如果i≤Maxcount，将f^(M)'根据步骤4的方程组进行下一次完整迭代，否则，停止迭代，输出重构图像。

本发明的基于径向TV的ART火焰切片重构算法，能够有效提高火焰图像的切片重构质量和重构速度；尤其在投影角度较少的情况下，该算法的去相关性更好，径向梯度图像更稀疏，且火焰重构过程稳定性好，算法鲁棒性强。

附图说明

图1为x_ρ,θ的相邻像素点的寻找过程示意图；附图标记为：(1-1)、寻找(1-2)、寻找x_ρ+,θ；(1-3)寻找(1-4)、寻找

图2为本发明的基于径向TV的ART火焰切片重构算法流程图；

图3为两幅模拟火焰切片的灰度图像，附图标记为：(3-1)、图像1；(3-2)、图像2；

图4为两个图像不同投影角度下重构图像的PSNR随着迭代次数的变化曲线，附图标记为：(4-1)、图像1在18角下；(4-2)、图像1在36角下；(4-3)、图像1在72角下；(4-4)、图像2在18角下；(4-5)、图像2在36角下；(4-6)、图像2在72角下；

图5为重构结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

压缩感知稀疏重构和代数迭代重构法都具有从少量投影数据重构出原始完全图像数据的特点，且汽油机燃烧火焰图像呈现亮度聚集块的稀疏结构形式，其在TV域呈现较好的稀疏特性，因此可以将全变差优化引入到ART算法中来对汽油机火焰进行切片图像重构，以提高火焰的切片的重构效果。

首先，建立径向梯度全变差模型：将火焰图像的中心设为原点，建立极坐标系，计算出火焰图像的径向梯度和径向梯度全变差；定义原点处梯度为0，对于像素点x_ρ,θ在极坐标系中的位置用M(ρ,θ)表示，其局部梯度表示为：

$\▿x_{\mathrm{\ρ},\mathrm{\θ}}=(D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}},D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}})---\left(1\right)$

其中D_ρ和D_θ分别是沿极径和极角方向的离散微分算子，有如下公式计算：

$\begin{array}{c}{D}_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{-},\mathrm{\θ}}\\ D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{-}}\end{array}---\left(2\right)$

由径向梯度求解径向梯度全变差，需要对极径ρ和极角θ两个方向图像像素点梯度向量的模进行相加，即梯度向量的模为图像像素沿极径ρ和极角θ两个方向梯度的分量的平方和的开方求和，计算公式如下：

$\begin{array}{c}{||x||}_{R-TV}=\underset{\mathrm{\ρ},\mathrm{\θ}}{\Σ}\sqrt{D_{\mathrm{\ρ}}{x_{\mathrm{\ρ},\mathrm{\θ}}}^2+D_{\mathrm{\θ}}{x_{\mathrm{\ρ},\mathrm{\θ}}}^2}\\ =\underset{\mathrm{\ρ},\mathrm{\θ}}{\Σ}\sqrt{{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{-},\mathrm{\θ}})}^2+{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{-}})}^2}\end{array}---\left(3\right)$

程序实现中，该算法将Sidky的导数逼近的全变差优化算法作为参考，然后将径向梯度TV调整用导数逼近思想按如下计算公式进行了实现，至此建立了径向梯度全变差模型。

$\begin{array}{c}\▿{||x||}_{R-TV}=\frac{\∂{||x||}_{R-TV}}{\∂x_{\mathrm{\ρ},\mathrm{\θ}}}\≈\frac{2x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{\_},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{\_}}}{\sqrt{{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{\_},\mathrm{\θ}})}^2+{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{\_}})}^2+\mathrm{\ϵ}}}\\ +\frac{x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{+},\mathrm{\θ}}}{\sqrt{{(x_{{\mathrm{\ρ}}_{+},\mathrm{\θ}}-x_{\mathrm{\ρ},\mathrm{\θ}})}^2+{(x_{{\mathrm{\ρ}}_{+},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{\_}})}^2+\mathrm{\ϵ}}}\\ +\frac{x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{+}}}{\sqrt{{(x_{\mathrm{\ρ},{\mathrm{\θ}}_{+}}-x_{\mathrm{\ρ},\mathrm{\θ}})}^2+{(x_{\mathrm{\ρ},{\mathrm{\θ}}_{+}}-x_{{\mathrm{\ρ}}_{-},{\mathrm{\θ}}_{+}})}^2+\mathrm{\ϵ}}}\end{array}---\left(4\right)$

径向梯度全变差模型的径向梯度计算过程中，对于像素点x_ρ,θ的相邻像素点的寻找过程如下：对于极径方向的离散微分，首先寻找极径减小1个单位的像素点集合集合中的点位于以坐标原点为圆心，半径为ρ-1的圆上。在这些像素点中，寻找像素向量与原像素点向量夹角最小的像素点即为需要的极径减少方向的梯度分量相邻像素点，表示为用上述相似的方法寻找极径增加方向的梯度分量相邻像素点，表示为对于极角方向的离散微分，首先寻找在极径不变的像素点集合集合中的点位于以坐标原点为圆心，半径为ρ的圆上。在这些像素点中，以原像素点x_ρ,θ为起点，顺时针寻找像素向量与原像素点向量夹角最小的像素点，即为需要的极角减少方向的梯度分量相邻像素点，表示为用上述相似的方法逆时针寻找极角增加方向的梯度分量相邻像素点，表示为按照上述过程及式(1)即可求解所需要径向梯度。

本发明的基于径向梯度全变差的ART火焰切片重构算法的具体的实现步骤如下：

步骤1：对火焰切片图像进行预处理：包括灰度变换和平滑去噪操作；

步骤2：求得火焰切片图像重建的迭代权重矩阵W并获得相应的投影图像矩阵P；

步骤3：设定初始迭代次数i＝1和迭代次数的最大值Maxcount＝15，并设置迭代的初始向量f⁽⁰⁾＝{0,0,...,0}；

步骤4：依据ART重构算法计算重构图像，一般采用Kaczmarz提出的松弛法来求解，选取松弛因子为恒值1，

$f^{\left(m\right)}=f^{(m-1)}+{\mathrm{\ω}}_{jk}(P_k-\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}{f}^{(m-1)})/\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}^2---\left(5\right)$

其中m＝1,2,…,M代表每一次完整迭代的方程组的方程个数；根据式(5)的方程组，计算M×N方程组；将初始向量f⁽⁰⁾＝{0,0,...,0}代入式(4)方程组的第一个方程求出f⁽¹⁾，由第二个方程求出f⁽²⁾，依次类推，最终由第M个方程求得f^(M)，这样就完成了一次完整迭代；

步骤5：根据上述径向梯度全变差模型，求解f^(M)的径向梯度

步骤6：按照导数逼近思想，求径向梯度全变差的全微分并对f^(M)进行径向梯度全变差优化，即采用梯度下降法进行调整如下式，其中q为下降步长因子，本发明中取值为0.1/k，并用f^(M)'代换原来的f^(M)

$f^{\left(M\right)}\′=f^{\left(M\right)}-q\·\▿{||f^{\left(M\right)}||}_{R-TV}---\left(6\right)$

步骤7：如果i≤Maxcount，将f^(M)'根据式(5)进行下一次完整迭代，否则停止迭代输出重构图像。

径向梯度全变差模型中径向梯度计算时，对于像素点x_ρ,θ的相邻像素点的寻找过程如图1所示。

本发明的具体实施例详细说明如下：

1)、拟选用两幅不同的模拟火焰切片图像如图3所示：

2)、拟选取初始迭代次数i＝1和迭代次数的最大值Maxcount＝15，下降步长因子q为0.1/k；

3)、然后按照上述技术方案的算法流程分别计算出火焰切片图像的径向梯度和径向梯度全变差，并对ART重构出来的重构图像进行优化调整，完成火焰切片重构过程；

4)性能测试

测试中本发明的ART-R-TV算法和传统的ART-TV算法的重构过程均是在18个、36个和72个投影角度下进行的。测试中得到的重构图像的MSE和PSNR的15次迭代的平均值的结果分别如表1所示。并对两种算法的PSNR在迭代过程中的变化曲线进行了比较如图4所示。

表1、图像1和图像2在三个投影角度下的重构结果的MSE和PSNR对比

由表1可以看出，两幅图像在重构过程中利用本发明的算法的重构质量均高于传统算法。在18、36和72个投影角度下迭代15次时，图像1的重构图像的PSNR值分别平均提高了2.93db、3.36db和3.68db，图像2的重构图像的PSNR值分别平均提高了3.77db、4.14db和4.34db。同时，两幅图像的重构图像的MSE对比得到本发明的算法相较于传统算法的MSE平均值更小且重构精度更高。

由图4的PSNR变化情况的对比，可以看出随着投影角度的增多，两种算法的重构图像的PSNR值都会有所提高。而传统ART-TV算法的重构图像的PSNR随着迭代次数的增加提高的有所起伏并不稳定，增长曲线也多为折线。本发明的ART-R-TV算法的重构图像的PSNR随迭代次数增加平滑增长更加稳定，保证了重构过程稳定在一个平均水平，使算法鲁棒性更强。同时，在相同投影角度数下达到相同PSNR时,本发明的算法所需要的迭代次数要少于传统的ART-TV算法。

表2、两个图像分别在18、36和72投影角度下的重构图像

Claims (1)

1.一种基于径向梯度全变差的ART火焰切片重构算法，其特征在于，该算法包括以下步骤：

步骤(1)、对火焰切片图像进行预处理：包括灰度变换和平滑去噪操作；

步骤(2)、求得火焰切片图像重建的迭代权重矩阵W并获得相应的投影图像矩阵P，具体步骤如下：

将二维火焰切片离散为N＝n×n个独立的体元，每个体元都具有特定的代表该单元亮度大小的值x_j(j＝1,2,…,N)，二维投影图像上的每个像素点都对应一条穿过该切片体元的射线L_k(k＝1,2,…,M)，M是所有投影的射线总数，等于投影角度个数与某角度投影射线数的乘积；而该像素点的值P_k(k＝1,2,…,M)为这条射线上所有体元对该射线亮度贡献的叠加；设ω_jk为体元j对射线L_k的权重因子，代表了体元与射线间的相关性，体元j对射线L_k的贡献为：

P_jk＝ω_jk·x_j

切片上所有体元对L_k亮度贡献和表示为：

$P_k=\underset{j=1}{\overset{N}{\Σ}}{P}_{jk}=\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}\·x_j,k=1,2,...,M$

用矩阵的形式表示如下：

Wx＝P

P为M维投影矩阵，表示二维投影的像素值；x为N维待重构切片数据矢量；W为M×N维投影权重矩阵；

步骤(3)、设定初始迭代次数i＝1和迭代次数的最大值Maxcount＝15，并设置迭代的初始向量f⁽⁰⁾＝{0,0,...,0}；

步骤(4)、计算重构图像，方程组公式如下：

$f^{\left(m\right)}=f^{(m-1)}+{\mathrm{\ω}}_{jk}(P_k-\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}{f}^{(m-1)})/\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}^2$

其中m＝1,2,…,M代表每一次完整迭代的方程组的方程个数；根据上述方程组公式，计算M×N方程组；将初始向量f⁽⁰⁾＝{0,0,...,0}代入方程组的第一个方程求出f⁽¹⁾，由第二个方程求出f⁽²⁾，依次类推，最终由第M个方程求得f^(M)，这样就完成了一次完整迭代；

步骤(5)、求解f^(M)的径向梯度公式如下：

$\▿f^{\left(M\right)}=(D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}},D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}})$

其中D_ρ和D_θ分别是沿极径和极角方向的离散微分算子，有如下公式计算：

D_ρx_ρ,θ＝x_ρ,θ-x_ρ-,θ

D_θx_ρ,θ＝x_ρ,θ-x_ρ,θ-

其中ρ_-和ρ₊分别代表沿极径减少和极径增加方向的变化，θ_-和θ₊则分别代表沿极角顺时针和极角逆时针方向的变化；

步骤(6)、求解f^(M)径向梯度全变差的偏导函数并对f^(M)进行径向梯度全变差优化，如下式：

$f^{\left(M\right)}\′=f^{\left(M\right)}-q\·\▿{||f^{\left(M\right)}||}_{R-TV}$

其中，q为下降步长因子，本发明中取值为0.1/k，并用f^(M)'代换原来的f^(M)；

步骤(7)、如果i≤Maxcount，将f^(M)'根据步骤(4)的方程组进行下一次完整迭代，否则，停止迭代，输出重构图像。