基于光偏折层析的脉动火焰四维测量方法

摘要

本发明涉及一种基于光偏折层析的脉动火焰四维测量方法，所述方法构建4D-CT多视角投影有序阵列偏折成像系统，针对被测流场的结构特点，确定最佳的采样角度，将激光源发出的激光束进行扩束准直和空间滤波器后，使用光束分束器进行二次分束，产生多条投射光束，进而生成莫尔偏折条纹阵列；使用基于波前恢复的条纹图处理技术，提取投影条纹的光偏折信息；构建基于偏折角修正重建的压缩传感重建算法对偏折投影信息进行反演计算，得到不同条纹图截面的二维参量分布；使用单阈值分割法对重建后的二维温度分布进行阈值分割，使用基于Marching Cubes面绘制算法和Ray Casting体绘制算法的Visualization Tool Kit技术，实现脉动燃烧流场及火焰结构的三维体数据动态显示。

说明书

技术领域

本发明属于激光技术应用领域，具体涉及一种基于光偏折层析的脉动火焰 四维测量方法。

背景技术

脉动燃烧是当代先进的非常规燃烧技术，它是指在燃烧和声脉动相互激励 和反馈的条件下，温度、压力、气流速度及热释放率等状态参数随时间周期性 波动的一种特殊的不稳定燃烧过程。脉动燃烧具有燃烧强度高、传热效率高、 低NOx排放等常规稳燃无可比拟的优势，对于节约能源、控制污染物排放、研 究高效动力机械等具有重要意义，因而受到了美国、日本、中国和欧洲多国的 高度重视。然而其燃烧机理、燃烧特性、应用技术的研究同样受制于流场结 构、参量分布、燃烧过程的定量测量和精确描述。

近几十年来，结合光学与激光技术、计算机技术、电子技术、信息处理技 术而产生的粒子图像测速技术(PIV)、激光诱导荧光技术(LIF)、激光全息干涉 技术(LIH)、激光－超声技术(LUM)、相干反斯托克斯喇曼散射(CARS)等新一 代流场显示与测量方法蓬勃发展，从而使复杂流动的研究取得突破性进展，这 其中具备非接触测量、多参量分析、定量计算等特质的光学计算层析 (Computed Tomography(CT))技术已成为复杂流场测试领域的有力工具。

光学CT(干涉CT、光偏折CT等)利用激光光波作为信息敏感器，获取光线 通过被测流场多方向的投影数据反演二维乃至三维分布。相比于干涉CT，光偏 折CT具有对光源要求不高、抗干扰能力强、光路简单易调等优点。更为重要 的是，对于高温、高速、非稳态相位流场，光偏折成像的低敏感性使它具备了 宽广的动态测试范围，成为其测量复杂流场的显著优势。

现有技术方案

光偏折CT理论发展和实际应用的关键技术主要有两方面：其一是复杂流 场多方向投影数据的获取，其二是由投影信息重建场分布。现有光偏折CT通 过投影采样、投影信息提取、重建计算等流程可以实现流场的二维乃至三维分 布。

(1)光偏折CT的成像采样

燃烧流场多方向光偏折成像的获取是实现流场三维结构呈现和参量定量测 量的首要条件。基于光线偏折原理测量的成像方法主要有纹影法和叠栅法。在 最新的研究中，基于背景导向纹影法(BOS)的光偏折CT流场测量不断取得突 破，Masanori等人分别使用色带图背景纹影(SPBOS)CT和彩色网格背景纹影 (CGBOS)CT实现了超音速流场的定量密度测量，赵玉新和田立丰等人也分别使 用BOS进行了复杂气动流场测量。然而，Miranda等人通过分析和实验证明，在 流场测量中条纹偏折法(FD)比BOS具有更高的测量精度，并且对应更少的数学 处理。因此对于复杂的非对称湍流流场测量，叠栅层析是更可胜任的方法。

在叠栅投影采样中，旋转扫描法光路简单、采样范围可控，并可实现全场 视角采样，我们之前已将其应用于发热元件生成的非对称温度场和风洞中高速 激波密度场测量。另一种分光路法能够实现多方向投影同时采集，可以对复杂 的瞬态流场进行记录和测量。Floyd等人和高益庆教授的团队基于该思想分别 开展了化学发光CT(CTC)和发射光谱层析(EST)流场测量。

(2)光偏折CT的重建算法

由投影信息反演被测流场一直是CT理论的核心和难点问题，光偏折成像的 一阶偏导数积分投影原理更直接增加了重建难度。在将投影数据进行各种处理 和转化后，以代数重建技术(ART)为代表的迭代类算法，因其在各种采样方 案、非完全数据、噪声控制等方面的鲁棒性决定了它比解析(变换类)算法更适 合解决实际的CT重建问题。近年来涌现的改进的序列子集联合代数重建技术 (ISO-SART)]、基于Newton-Raphson算法的加速代数重建技术(AART)、定向代 数重建技术(DART)、自适应联合代数重建技术(ASART)、改进的自适应代数重 建技术(MAART)、自适应统计迭代重建技术(ASIRT)等研究不断改良着迭代类 算法的重建性能。然而这些未考虑CT成像原理和特点的通用CT重建算法，无 法体现光偏折CT本质(成像原理)上的优势。鉴于此，贺安之教授的团队从基本 的偏折公式出发，推导出偏折角修正重建技术(DARRT)，并进行了多种复杂流 场的光偏折CT测量。

另一方面，燃烧非稳态流场的光偏折投影获取方法决定了CT重建只能在少 数方向投影采样的条件下进行，这进一步增加了重建算法的难度。此类欠采样 的非完全数据重建是近年来各种CT理论研究领域极为活跃的课题，Applied Physics Letters和Optics Express等刊物上都发表了使用有限视场数据反演层析图 像的最新成果。Liu等人使用结合最速下降策略的总变差最小化算法(TVM-SD) 和结合软阈值滤波的全差分最小化算法(TDM-STF)重建有限角投影CT图像，并 通过实验比较了两者在有限角投影重建图像中的效能。Ambartsoumian和 Nquyen等人分别使用改进的最大期望值算法重建非完全数据CT图像。高益庆 等人将最大熵、惩罚函数等技术与迭代算法相结合并应用于非完全数据光学CT 重建。

对于脉动燃烧等具有固定频率的振荡燃烧过程，现有光偏折CT二维乃至三 维的稳态和瞬态测量，无法实现脉动燃烧湍流结构的实时显示和参量分布周期 变化的定量描述，光偏折CT的关键方法亟待发展。

(1)光偏折CT的成像采样

该方案采用固定被测流场旋转光路系统(或固定光路系统旋转被测流场) 的方法获取不同视角的投影条纹，但此方法的缺陷是无法真正实现对同一流场 过程的同时采样，因此只适用于过程近似稳定且可重复的流场测量。分光路法 的能够实现多方向投影同时采集，但将该原理应用于光偏折CT时，不仅需要 多台CCD摄像装置并增加相应的成像元件，而且很难保证光路严格复制，从而 增大测量误差。因此，现有的光偏折CT投影获取方法无法满足诸如脉动燃烧 等非稳态复杂流场的采样要求，直接影响了光偏折CT在复杂流场测量领域的 实用化。另一方面，脉动燃烧具有周期性的湍流结构和参量分布变化，这一重 要的燃烧过程特性也对光偏折CT的投影采样提出了更高的要求。

(2)光偏折CT的重建算法

尽管针对于计算层析的代数重建技术不断发展，然而这些未考虑CT成像 原理和特点的通用CT重建算法，无法体现光偏折CT本质(成像原理)上的优 势。鉴于此，贺安之教授的团队从基本的偏折公式出发，推导出偏折角修正重 建技术(DARRT)，并进行了多种复杂流场的光偏折CT测量。另一方面，燃烧 非稳态流场的光偏折投影获取方法决定了CT重建只能在少数方向投影采样的 条件下进行，这进一步增加了重建算法的难度，此类欠采样的非完全数据重建 是近年来各种CT理论研究领域极为活跃的课题。前期工作中，我们建立了光 偏折CT重建的Tikhonov正则化方法，并讨论了该方法对有限角投影数据的重 建效能。然而对极度欠采样条件下的重建仍不能满足复杂结构湍流燃烧高精度 测量的要求。

发明内容

本发明克服现有技术的不足，提出了一种基于光偏折层析的脉动火焰四维 测量方法，所述方法克服了现有分光路系统需要多组成像元件和多台CCD摄像 装置、难于保证光路严格复制等缺陷，在CCD单幅图像上完成即时6方向投影 有序阵列可视化，实现多视角、同时、同光路条件采样，对全周期非稳态脉动 燃烧流场进行显示和动态采样。

本发明的技术方案为：

一种基于光偏折层析的脉动火焰四维测量方法，所述方法首先构建4D-CT 多视角投影有序阵列偏折成像系统，4D-CT多视角投影有序阵列偏折成像系统 包括反射镜、空间滤波器、激光源、透镜、光束分束器、光栅、毛玻璃屏、高 速摄像系统，然后依据下列步骤重建燃烧流场和火焰结构：

1)针对被测流场的结构特点，确定最佳的采样角度，将激光源发出的激光 束进行扩束准直和空间滤波器后，使用光束分束器进行二次分束，产生多条投 射光束，改变各光束方向，形成精确可控的多视角探测光束，经过被测流场； 设计可调的反射镜阵列，将空间调制后的被测流场多视角光束形成同方向有序 光束阵列，平行通过叠栅，形成莫尔偏折条纹阵列；使用高速摄像系统实现莫 尔条纹阵列的可视化与采样；

2)使用基于波前恢复的条纹图处理技术，提取投影的光偏折信息；

3)构建基于偏折角修正重建的压缩传感重建算法对偏折投影信息进行反演 计算，得到不同条纹图截面的二维参量分布；

4)使用Sobel算子对重建后的二维温度分布进行单阈值分割，确定燃烧火 焰边界，多幅阈值分割后的二维温度分布按位置顺序排列，使用基于Marching Cubes面绘制算法和Ray Casting体绘制算法的Visualization Tool Kit技术，实现 脉动燃烧流场及火焰结构的三维体数据动态显示。

本发明的有益效果是：

(1)本发明克服了现有分光路系统需要多组成像元件和多台CCD摄像装 置、难于保证光路严格复制等缺陷，通过投影阵列偏折仪，在CCD单幅图像上 完成即时6方向投影有序阵列可视化，实现多视角、同时、同光路条件采样， 对全周期非稳态脉动燃烧流场进行显示和动态采样。

(2)本发明将压缩传感理论与偏折角修正迭代技术相结合解决光偏折层析 极少方向的投影重建；

(3)本发明使用折射率(或温度、密度)阈值确定火焰边界，使用基于面 绘制及体绘制方法的VTK技术，实现三维火焰表面结构和参量数据场显示，从 而在周期内形成脉动火焰的四维结构和参量显示。

附图说明

下面结合附图对本发明进一步说明。

图1是本发明结构示意图；

图2是本发明光线通过流场产生偏折及坐标系统示意图。

图中：1、反射镜；2、空间滤波器；3、激光源；4、透镜；5、光束分束器；6、 光栅；7、被测流场：8、毛玻璃屏：9、高速摄像系统；10、光束。

具体实施方式：

参见图1和图2所示，本发明构建4D-CT多视角投影有序阵列偏折成像系 统，针对被测流场的结构特点，确定最佳的采样角度，将激光源发出的激光束 进行扩束准直和空间滤波器后，使用光束分束器进行二次分束，产生多条投射 光束10，进而生成莫尔偏折条纹阵列；使用基于波前恢复的条纹图处理技术， 提取投影条纹的光偏折信息；构建基于偏折角修正重建的压缩传感重建算法对 偏折投影信息进行反演计算，得到不同条纹图截面的二维参量分布；使用单阈 值分割法对重建后的二维温度分布进行阈值分割，使用基于Marching Cubes面 绘制算法和Ray Casting体绘制算法的Visualization Tool Kit技术，实现脉动燃烧 流场及火焰结构的三维体数据动态显示。

具体为：首先构建4D-CT多视角投影有序阵列偏折系统，4D-CT多视角投 影有序阵列偏折系统包括反射镜1、空间滤波器2、激光源3、透镜4、光束分 束器5、光栅6、毛玻璃屏8、高速摄像系统9和光束10，然后依据下列步骤提 取光偏折投影信息并进行流场重建。

1)针对被测流场7的结构特点，确定最佳的采样角度，使用光束分束器5 将扩束准直的由激光源发出的激光束经过空间滤波器后，进行二次分束，产生 多条投射光束10，改变各光束方向，形成精确可控的多视角探测光束，经过被 测流场7；设计光学系统如图1所示，波长为632.8nm的氦氖激光源3发出激光， 经过空间滤波器2扩束、滤波，经过透镜4后，成为直径3～5cm的平行光束。 调整六个光束分束器5，光束分束器5为半透半反镜，调整三个平面反射镜1， 将平行光束分成角度为10°、30°、55°、100°、136°和152°的六条光束通过被测 流场7，并发生微小偏折，再由10个平面反射镜1将六条光束调整成为3(列) ×2(行)的同轴光束，通过2片尺寸为178×178mm平行放置的Ronchi光栅6， 在光栅后的毛玻璃屏8上，将形成3(列)×2(行)的莫尔条纹阵列，使用分辨 率为1024×768的CCD摄像机(C)，可以采集到动态的光偏折投影(莫尔条纹) 阵列图像。

可调的反射镜阵列，将空间调制后的被测流场多视角地发射光束形成同方 向有序光束阵列，经过缩束后，通过叠栅；使用高速摄像系统9实现投影条纹 有序图阵的可视化与采样，使用基于波前恢复的条纹图序阵批处理技术，平行 光通过相位物体后，形成的畸变莫尔条纹分布g(y,z)可以表示成空间相位调制信 号：

$g(y,z)=r(y,z)\underset{-\infty }{\overset{\infty }{\Sigma }}{A}_n\exp \left\{i\right[2\mathrm{\pi n}{f}_{0y}y+2\mathrm{\pi n}{f}_{0z}+\mathrm{n\phi }(y,z)\left]\right\}$

其中r(y,z)表示照射光束通过相位物体后的光强分布，φ(y,z)为莫尔条纹畸变引 起的空间相位改变量，f_0y、f_0z为f₀的两个分量。对上式做二维傅立叶变换得到 傅立叶频谱，用滤波函数提取+1级频谱，并通过平移到频率面原点去除参考频 率f₀，再进行逆傅立叶变换得到：

g₁(y,z)＝A₁r(y,z)exp[iφ(y,z)]

从而可以求得：

${\phi }_w(y,z)\arctan \frac{\mathrm{Im}\left[g_1(y,z)\right]}{\mathrm{Re}\left[g_1(y,z)\right]}$

其中，Im[]和Re[]分别表示对g₁(y,z)做虚部和实部运算。此时得到的φ_w(y,z)被 包裹在-π:π之间，它是真实波面φ(y,z)的2π的模。接下来沿y,z两个方向对 φ_w(y,z)积分，即在包裹的相位图上对发生2π→0相位跃变处后的各点相位增加 2π和对发生0→2π相位跃变处后的各点相位减少2π，这样得到的φ(y,z)是偏折 的莫尔条纹相对于理想直条纹的相位改变量。但在实验中要得到严格的直条纹 是不可能的，因为它要求照射光为理想的平行光以及两光栅严格和入射光垂直。 因此在没有相位物体时，初始莫尔条纹已经存在相位φ₀(y,z)。去除初始相位，

△φ(y,z)＝φ(y,z)-φ₀(y,z)

得到变形波面的真实相位△φ(y,z)。将△φ(y,z)代入下式即可求得偏折角

式中θ为两光栅夹角，d_c为光栅节距。

波前恢复方法即相位场二维投影数据的相位展开方法，它克服了基于条纹 跟踪的莫尔图形态学处理方法只能得到几个孤立层面一维偏折数据的缺陷。在 二维投影数据的基础上，充分利用空间相位场三维连续和相关的特点，即各层 析面上的场分布在轴向是连续且光滑过渡的，减少了采样的投影数据间的不相 容性，提高了相位场定量测试的精度。

2)构建基于偏折角修正重建的压缩传感重建算法对偏折投影信息进行反演 计算，得到不同条纹图截面的二维参量分布。

光偏折CT技术是基于光线在非均匀折射率场中传播发生方向改变而产生光 程差的复杂流场测试技术。光偏折CT是利用激光光波作为信息敏感器，获取光 线通过被测场多方向的偏折投影数据反演二维乃至三维场分布。光线通过流场 产生偏折的基本原理如图2，图中二维相位场O的折射率分布为n(x,y)，环境折射 率n₀，R.O.表示入射光线,x'-y'坐标轴相对x-y轴旋转了角度θ。在角度为θ的观 测方向，偏折角可以由偏折仪显示和分析得到。

式(1)中n是被测折射率分布，n₀是环境折射率。将折射率场划分成I＝w×w个 足够小的相等网格，网格边长为d。设每一网格内的折射率变化△n为常数，若有 l的投影方向，每个投影方向有m条射线，则该平面内任意光线的偏折角可以写 成

式中A_ij是第j条光线通过第i个网格的长度，是第j光线的偏折角。

根据压缩传感(CS)理论和偏折角CT原理，结合物理量(温度、密度) 分布的全变差的稀疏性，提出了一种偏折角压缩传感修正重建算法，并推导得 出算法的约束最小化的数学描述

$\underset{n}{\min }{\left|\right|n\left|\right|}_{\mathrm{TV}}s.t.\mathrm{Dn}=\mathrm{Cn}+q-h---\left(3\right)$

式(3)中

${\left|\right|n\left|\right|}_{\mathrm{TV}}=\underset{i,j}{\Sigma }\sqrt{|n(s+1,t)-n(s,t)|+|n(s,t+1)-n(s,t)|}---\left(4\right)$

n是折射率分布的矩阵，s,t是矩阵n的坐标序号，h,q,C,D的元素是光线入射 斜率k和A的函数，A是光线通过网格的长度矩阵，q和h分别代表迭代过程中的 边界修正量和各采样偏折角。计算过程分为以下三步。

第一步，对式(3)线性转化的偏折角投影方程使用偏折角修正思想进行迭代 重建。将光偏折角公式折射率偏导数进行离散、合并相关同类项等一系列处理， 推导得出偏折角修正重建算法(DARRT)的迭代表达式：

$n_i^{m+1}=n_i^m+\omega \frac{(h_j-q_j+\underset{i=0}{\overset{I-1}{\Sigma }}{D}_{\mathrm{ij}}{n_i}^m)-\underset{i=0}{\overset{I-1}{\Sigma }}{C}_{\mathrm{ij}}{n_i}^m}{{\left(\underset{i=0}{\overset{I-1}{\Sigma }}{C}_{\mathrm{ij}}\right)}^2}{C}_{\mathrm{ij}},j=0,...,(m-1)\times l---\left(5\right)$

式中ω为松弛系数，i,j分别代表网格序号和采样序号，C_ij,D_ij,q_j,h_j为一个表 达式，表达式中设有A_(i+1)j,A_(i+w)j,A_(i-1)j,A_(i-w)j，A_(i+1)j,A_(i+w)j,A_(i-1)j,A_(i-w)j是第j条 光线通过第i个网格的相邻网格的长度

第二步，利用CS算法即结合全变差范数和最速下降法修正目标值。利用最 速下降法快速减小全变差范数，因而使迭代过程更快收敛或达到目标值。

为了防止梯度过大引起算法发散同时保证算法加速收敛，引入最速下降因 子ρ＝0.0000005，每迭代一次后调整因子为ρ_s＝0.997。

$\beta =\max \left(\right|n_{i,j}^{\left(k\right)}\left|\right)/\max \left(\right|f_{i,j}\left|\right)---\left(6\right)$

$n_{i,j\left(\mathrm{CS}\right)}^{\left(k\right)}=n_{i,j}^{\left(k\right)}-\beta \times \rho \times f_{i,j}---\left(7\right)$

ρ＝ρ×ρ_s   (8)

式(6)和式(7)中

$f_{i,j}=\frac{\partial n_{\mathrm{TV}}}{\partial n_{i,j}^{\left(k\right)}}=\frac{4n_{i,j}^{\left(k\right)}-n_{i+1,j}^{\left(k\right)}-n_{i-1,j}^{\left(k\right)}-n_{i,j-1}^{\left(k\right)}}{a_{i,j}}+\frac{n_{i,j}^{\left(k\right)}-n_{i+1,j}^{\left(k\right)}}{a_{i+1,j}}+\frac{n_{i,j}^{\left(k\right)}-n_{i-1,j}^{\left(k\right)}}{a_{i-1,j}}+\frac{n_{i,j}^{\left(k\right)}-n_{i,j+1}^{\left(k\right)}}{a_{i,j+1}}+\frac{n_{i,j}^{\left(k\right)}-n_{i,j-1}^{\left(k\right)}}{a_{i,j-1}}---\left(9\right)$

式(9)中

$a_{i,j}=\sqrt{\frac{(n_{i+1,j}^{\left(k\right)}-n_{i,j}^{\left(k\right)})+(n_{i,j}^{\left(k\right)}-n_{i-1,j}^{\left(k\right)})+(n_{i,j}^{\left(k\right)}-n_{i,j-1}^{\left(n\right)})+(n_{i,j+1}^{\left(k\right)}-n_{i,j}^{\left(k\right)})}{2{\Delta }^2}+{ϵ}^2}---\left(10\right)$

式中△是采样间隔，ε是很小的正数，本文取ε＝10^-8。

第三步，初期数值模拟实验显示，在极少数投影条件下，DARRT重建的流 场分布较平滑，但峰值误差较大，而CS重建的流场能很快达到峰值，但重建流 场毛刺较多。为此引入一个权重因子λ(0≤λ≤1)，通过调整λ的大小来发挥 DARRT和CS优势，同时抑制DARRT和CS各自的不足。迭代过程：

$n_{\mathrm{DARRT}+\mathrm{CS}}^{\left(k\right)}=(1-\lambda )n_{\mathrm{DARRT}}^{\left(k\right)}+\lambda n_{\mathrm{cs}}^{\left(k\right)}---\left(11\right)$

3)由DARRT+CS算法迭代得到二维温度分布，使用单阈值分割法对重建 后的二维温度分布进行阈值分割，用Sobel算子计算二维温度分布的梯度，选取 梯度域确定燃烧火焰边界。多幅阈值分割后的二维温度分布按位置顺序排列， 使用基于面绘制及体绘制方法的VTK(Visualization Tool Kit)技术，实现三维 火焰表面结构和参量数据场显示。其中，面绘制使用Marching Cubes算法(移 动立方体算法)，体绘制采用Ray Casting算法(光线投射法)，利用VTK函数 库结合C++语言进行编程实现。当完成周期内所有时刻的重建，即形成脉动火 焰的四维结构和参量显示序列，其原理如下。

对重建后的二维温度场施加sobel算子，获得温度场的梯度分布。

该算子包含两组3×3的矩阵，分别为横向及纵向，将之与温度场作平面卷 积，即可分别得出横向及纵向温度梯度。如果以A代表原始温度场，Gx及Gy 分别代表经横向及纵向边缘检测的温度分布，其公式如下：

温度场每一个点的横向及纵向梯度通过以下公式结合，来计算该点梯度的 大小G：

$G=\sqrt{G_x^2+G_y^2}$

在由sobel算子计算得出的梯度分布中选取最大的梯度域确定燃烧火焰的边 界，即对原始温度场进行阈值分割。把多幅阈值分割后的二维温度分布按位置 顺序排列获得用于三维重建的温度数据场，这种由多幅二维温度分布数据所组 成的体数据属离散的结构化数据(所谓结构化数据是指在逻辑上组成三维数组的 空间离散数据，每个元素都有确定的层号、行号、列号)。通常结构化数据又分2 种：均匀网格结构化数据和规则网格结构化数据。前者在X、Y、Z等3个方向上 网格之间的距离都相等，后者在X、Y、Z等3个方向的距离互不相等，但在同一 个方向上每层间隔的距离是相等的。光学CT所获得的数据属后者，即规则网格 结构化数据，在X方向的最小间隔为单层温度场离散点间的横向距离，Y方向的 最小间隔为单层温度场离散点间的纵向距离，Z方向的最小间隔为每层温度数据 场的距离。

在规则网格结构化数据中，每个网格称为1个体素(Voxe1)。每个小立方体代 表1个体素。对于每个体素的8个顶点，可根据具体情况定义其函数值，这个函 数值可以是标量值、灰度值、光强度值或透明度值等。相应的光学CT所获得的 三维空间数据场又称为光学CT获得的体数据，每个顶点的函数值也称为其体数 据的函数值，在本专利中函数值为温度。

Marching Cubes算法的基本原理:在三维数据场中构造等值面,找出经过该等 值面的体元(cubes)，求出该体元内的等值面并计算出相关参数,便于使用常用的 软件包或图形硬件绘制出等值面。

MarchingCubes算法的过程可以描述如下:

1)从三维体数据中每次读出2张切片，形成一层(layer)。

2)2张切片上下相对应的4个点构成一个立方体(cube)。

3)从左至右、从前到后的顺序处理一层中的立方体(抽取每个立方体中的等 值面)，然后从下到上顺序处理到n-1层，算法就结束，故名为Marching Cubes。 对于每一个立方体而言，它的8个顶点的灰度值可以直接从输入数据中得到,要抽 取的等值面的阈值由用户指定，然后根据体数据的信息就可以提取出等值面的 三角网格表达，进而对物体表面进行显示。

Ray Casting算法：一种基于图像空间直接体绘制的经典算法，对体数据的颜 色和透明度进行重采样后通过累加得到重建图像。应用射线穿过数据场时颜色 和透明度累加的Alpha混合算法，即沿射线的方向对采样物体进行颜色和透明度 的从前至后的累加算法。

光学CT三维温度场重建中颜色和透明度累加的Alpha混合算法：由于漫反 射、环境光和镜面反射光的作用，透明物体本身反射出光线，但同时物体也会 遮挡一部分光线，透明物体反射光和衰减光的能力可分别用颜色和透明度来表 示。透明度的本质是指光穿过物体的能力，光穿过物体波长会发生变化。光穿 过三维数据集时，间隔均匀采样过程中光线波长的变化是累加的，透明物体的 渲染本质上是将透明物体的颜色和其后物体的颜色进行混合。在使射线穿过数 据场的光学CT三维温度场重建中，透明度和颜色的累加应用Alpha混合算法实 现。假设间隔采样的一个采样点是A，沿光线方向与A相邻的一个采样点是B，那 么透过B去看A，其视觉效果可以用C来表示，C即是B和A用Alpha混合算法计算 的结果。Alpha混合算法的计算公式可近似表示为：

I(C)＝α(B)*I(B)+(1-α(B))*I(A)

α(C)＝α(B)+(1-α(B))*α(A)

其中，α(A)、α(B)和α(C)是采样点A和B及其视觉效果C的透明度，I(A)、 I(B)和I(C)分别是采样点A和B及其视觉效果C的光强度值。透明度的取值范围 是从0到1，0代表完全透明，1代表完全不透明。

颜色和透明度的累加有2种方法：一种是沿射线从前向后累加；另一种是沿 射线从后向前累加。从前向后累加的计算效率要高于从后向前累加，在从前向 后累加的过程中，一旦颜色达到最大值，透明度值累加到完全不透明就可以结 束采样，可提高采样的效率。因此，在使射线穿过数据场的光学三维温度场重 建方法中，选用从前向后累加Alpha混合算法进行颜色和透明度的计算。

VTK(Visualization ToolKit)是美国Kitware公司开发的一套C++类库，它吸收 了众多优秀的图像处理和图形生成算法，是一个源代码开放、资源共享的软件。 分别使用VTK中的类vtkMarchingCubes和vtkvolumeRaycastFunction实现面绘制 和体绘制。