CN101178814B - 一种融合解剖与功能成像信息数据场的半透明体绘制方法 - Google Patents

Abstract

本发明公开了一种融合解剖与功能成像信息数据场的半透明体绘制方法，该方法包括用于半透明体绘制的光学模型TVROM及其实现方法，以及基于小波变换，以区域标准差与区域能量相结合的融合规则进行解剖与功能成像信息融合的方法。用于半透明体绘制的光学模型TVROM在光线吸收与发射模型的基础上，考虑了阴影、直接散射与间接散射等因素。在医学三维可视化中，用本发明能方便、经济、高效率地实现融合解剖与功能成像信息的三维数据场的半透明体绘制，并能充分显示重建对象的内部隐含分界面及内部细节的详细信息，从而满足医学运用的精确性要求。

Description

一种融合解剖与功能成像信息数据场的半透明体绘制方法

技术领域

本发明涉及一种融合解剖与功能成像信息数据场的半透明体绘制方法，属于医学图像处理及应用技术领域。

背景技术

在医学图像可视化的运用中，医生常常需要透过组织表面看清内部的病变结构和细节信息，确定肿瘤的相对位置。因此，能透过对象表面，显示内部的隐含分界面(内部不同介质之间的界面)及其内部细节的半透明体绘制技术被广泛应用于医学可视化中。

直接体绘制技术直接应用体绘制方程将体数据映射为二维图像，其实质是一个三维离散数据场的重采样和图像合成的过程。光学模型描述三维数据是如何产生、反射、阻挡以及散射光线的。基于各类光学模型，可以得到相应的体绘制方程。体绘制方程在重采样的基础上，依据光学模型计算出全部采样点对屏幕像素的贡献，亦即通过采用不同的方法累积数据属性来获得不同的绘制效果。

结合Blinn-Phong面明暗化模型的光线吸收与发射模型是应用最广泛的一种体绘制光学模型。尽管该模型能基本描述许多明显的半透明组织表面，然而，其难以表述半透明组织或者散射占优的组织的详细信息。Nelson Max在论文《Optical models for direct volumerendering》(IEEE Transactions on Visualization and Computer Graphics，1995，1(2)：99-108)证明了多重散射对于体绘制的效果显著，并且给出了计算多散射的方法。然而，庞大的计算量与苛刻的内存要求使得精确的光学模型难以在实时的医学工程中得到应用。因此，在实际的体绘制应用中，一些相对容易实现的光学经验模型被普遍采用。如，蔡文立和彭延军分别在论文《基于输运方程的混合式体绘制模型》(计算机学报，1995，18(5)：330-338)与论文《体绘制中显示隐含分界面的一种方法及其实现》(软件学报，2002，13(9)：1887-1892)中提出了在基本的光线吸收与发射模型中加入直接散射部分，并用非真实感绘制技术加大隐含分界面的显示效果。此外，Kniss Joe在论文《A model for volumelighting and modeling》(IEEE Transactions on Visualization and Computer Graphics，2003，9(2)：150-162)提出了一个适用于高档硬件实现的体绘制的光学经验模型。LacrouteP.在论文《Fast Volume Rendering Using a Shear-Warp Factorization of the ViewingTransformation》(Ph.D.dissertation，USA.Stanford University，1995：29-43)提出的Shear-Warp算法是目前最快的基于CPU的体绘制方法。

目前已申请的有关体绘制方法的发明中，体绘制的方法(02800260.1)、一种用于虚拟手术的体绘制成像方法(200410018565.9)、具有封闭轮廓多平面重定格式的体绘制数据的可视化(200510065505.7)、采用多图形处理器加速方格片元过滤的实时体绘制方法(200410086020.1)、CT图像的快速渐进式直接体绘制三维重建方法(200510042734.7)、滑动纹理的体绘制(200510087829.0)、医学图像中基于GPU硬件加速的体绘制方法(200510110665.9)、基于分形的体绘制可视化方法(200610117049.0)、空腔性脏器内壁快速体绘制方法(200610118894.X)、基于分块的蒙特卡洛体绘制方法(200610117570.4)，这些发明都没有涉及体绘制的本质——光学模型的内容。然而，直接体绘制技术直接应用体绘制方程将体数据映射为二维图像，其实质是一个三维离散数据场的重采样和图像合成的过程。光学模型描述三维数据是如何产生、反射、阻挡以及散射光线的。基于各类光学模型，可以得到相应的体绘制方程。因此，光学模型决定着最终的三维重建结果的效果。而使用明暗处理的直接体绘制(200580029305.4)则是用梯度参数对常用的Phong光照模型作了简单的补充。在三维数据场的半透明体绘制中，为了充分显示重建对象的内部隐含分界面及内部细节的详细信息，除了需要考虑光线的吸收、发射，还要考虑阴影、直接散射与间接散射等因素。

此外，在医学图像三维可视化中，体绘制所绘制的体数据来自于各类医学影像设备采集的多模态医学图像。电子计算机断层扫描(Computerized Tomography CT)、核磁共振成像(Magnetic Resonance Imaging MRI)等解剖图像以较高的分辩率提供了人体的解剖形态信息，但无法反映脏器的功能代谢信息；正电子发射断层成像(Positron EmissionTomography PET)、单光子发射计算机断层成像(Single Photon Emission ComputerizedTomography SPECT)等功能图像反映了脏器的代谢水平和血流状况，对肿瘤病变呈现“热点”，但无法提供脏器的解剖细节。因此，为了满足医学运用的精确性要求，要求三维数据场融合有解剖与功能成像的信息。

目前已申请的涉及到解剖与功能成像信息融合相关内容的发明中，容积—容积融合的可视化(03805385.3)提出利用多个容积的方法来融合解剖与功能成像信息的方法；一种多模态医学体数据三维可视化方法则侧重于研究配准的方法。上述两个发明都没涉及解剖与功能成像信息融合的融合规则问题。然而，为了满足医学运用的精确性要求，要求三维数据场融合有解剖与功能成像的信息。而在小波域里，融合解剖与功能成像的信息，需要考虑融合规则问题。

发明内容

本发明的目的在于克服现有技术的不足，提供一种能在普通硬件条件下实现的融合解剖与功能成像信息数据场的半透明体绘制方法。

本发明的目的通过如下技术方案实现：

一种融合解剖与功能成像信息数据场的半透明体绘制方法，包括如下步骤：

(1)形成融合了解剖与功能成像信息的三维数据场：基于小波变换，以区域标准差与区域能量相结合的融合方法进行解剖与功能成像的信息融合；在小波域里，采用区域标准差和能量作为小波系数的活动性水平测度；

(2)应用半透明体绘制的光学模型TVROM对融合解剖与功能成像信息数据场进行半透明体绘制：

a、再建游程编码数据结构；

b、在错切空间里，利用光学模型TVROM，合并中间图像；

c、在变形空间里，用纹理映射绘制最终图像。

所述步骤(1)基于小波变换，以区域标准差与区域能量相结合的融合方法进行解剖与功能成像的信息融合包括如下步骤：

①对CT、PET图像进行配准；

②选取双正交滤波器组的Daubechies 9/7双正交小波分别对配准后的CT和PET进行三层双正交小波分解；

③分别计算CT和PET图像中窗口的小波系数平均值；

④在图像的小波变换域内，对图像进行融合；

⑤对CT图像和PET图像经小波变换之后的逼近系数进行求和后除以2的处理；

⑥利用小波系数中的逼近系数进行小波反变换可得融合图像。

所述快速再建游程编码数据结构是指，在交互式操作中，当不透明度传递函数改变时，在对图像进行重采样前，利用不透明度传递函数的相关性，和已存在的min-max Octree分类数据结构，对体数据进行加速分类；随后，利用再分类后的分类编码数据结构，快速再建游程编码数据结构。

所述在错切空间里，利用光学模型TVROM，合并中间图像包括如下步骤：

a)预设置传递函数，将其设为查询表；

b)预设置直接散射的权重值；

c)合成中间图像：由前向后序次，循环对每个切片进行处理，合成中间图像，在处理第k个切片时，依次计算直接散射部分、计算光照累积光亮度、计算间接散射部分；将下一个切片的间接光照缓存数据存放于当前切片的间接光照缓存，并将下一个切片的间接光照缓存都设置为值1；处理下一个切片。

所述用于半透明体绘制的光学模型TVROM为：

$I(x_1,\stackrel{\→}{\mathrm{\ω}})=T(0,l)I(x_0,\stackrel{\→}{\mathrm{\ω}})+$

${\∫}_0^{l}T(s,l)*(S^D\left(s\right)*L_g*\exp (-{\∫}_s^{l_g}\mathrm{\τ}\left(t\right)\mathrm{dt})+S^I\left(s\right)*L_g*\exp (-{\∫}_s^{l_g}{\mathrm{\τ}}_i\left(x\right)\mathrm{dx}))\mathrm{ds}$

式中，

指沿着视线方向

在空间坐标x处的光线强度，

指视线方向，x为体素的空间三维坐标(x，y，z)，s、l表示沿着视线距离，T(s，l)指沿着视线从x(s)到x(l)过程中的衰减， $T(s,l)=\exp (-{\∫}_s^l\mathrm{\τ}\left(t\right)\mathrm{dt}),$ τ指光线衰减系数；x(s)指沿着视线距离s处的三维坐标；S^D(s)为直接散射的值，S^I(s)为间接散射的值；L_g为光源的光亮，其下角g表示沿着ω′的光线方向，ω′为入射光方向；l_g表示沿着光线方向的距离，其下角g表示沿着ω′的光线方向；τ_i指光线间接衰减系数，T_i(s，l)指沿着光线从x(s)到x(l)过程中的间接衰减， $T_i(s,l)=\exp (-{\∫}_s^l{\mathrm{\τ}}_i\left(t\right)\mathrm{dt}).$

相对于现有技术本发明具有如下优点和有益效果：

1)用于半透明体绘制的光学模型TVROM在光线吸收与发射模型的基础上，考虑了阴影、直接散射与间接散射等因素，结合该模型与改进Shear-Warp方法的半透明体绘制技术，能高效率地重建三维医学图像，并能充分显示内部隐含分界面及内部细节的详细信息。

2)能在普通硬件条件下实现的用于半透明体绘制的光学模型TVROM。在医学三维可视化中，用本发明能方便、经济、高效率地实现融合解剖与功能成像信息的三维数据场的半透明体绘制。

3)采用基于小波变换的融合方法，以区域标准差与区域能量相结合的融合规则对CT、PET图像进行融合处理，能有效地保留了多模态图像的边缘和纹理特征，避免了融合图像出现模糊现象，为PET图像描述的病变区域的准确定位提供了解剖信息。

附图说明

图1为根据本发明——一种融合解剖与功能成像信息数据场的半透明体绘制方法的流程图。

图2a为用本发明提供的用于半透明体绘制的光学模型TVROM及实现方法对头部进行半透明体绘制的结果。

图2b为光学模型TVROM的直接散射部分结果。

图2c为光学模型TVROM的间接散射部分。

图2d为用本发明提供的用于半透明体绘制的光学模型TVROM及实现方法对脑颅进行半透明体绘制的结果。

图3a为从正面观察半透明三维显示的靶区与病人内部结构的结果。

图3b为从背面观察半透明三维显示的靶区与病人内部结构的结果。

图4a为体绘制结果。

图4b为图4a白线所指向的横断面(融合解剖与功能成像信息数据场的横断剖面)。

图4c为融合解剖与功能成像信息数据场的横断剖面所对应的CT图像。

图4d为融合解剖与功能成像信息数据场的横断剖面所对应的PET图像。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明，但本发明要求保护的范围并不局限于实施例表述的范围。

本发明实验平台为PentiumIV 3.0G，1G RAM，显卡为NVIDIA Quadro FX 1400。断层序列图像物理设备生成，经PACS系统路由到本实验平台。本实验平台用本发明方法对三维数据场数据进行处理，最终图像在显存里用纹理绘制方式绘制。

实施例1

用本发明提供的能在普通硬件条件下(本发明的实验平台为PentiumIV 3.0G，1G RAM，显卡为NVIDIA Quadro FX 1400)实现的用于半透明体绘制的光学模型TVROM及其实现方法，对医学序列CT图像(512*512*377(12位))进行体绘制。实施例1效果如图2a、图2b。

用于半透明体绘制的光学模型TVROM在光线吸收与发射模型的基础上，考虑了阴影、直接散射与间接散射等因素，结合该模型与改进Shear-Warp方法的半透明体绘制技术，能高效率地重建三维医学图像，并能充分显示内部隐含分界面及内部细节的详细信息。

用于半透明体绘制的光学模型TVROM如下：

$I(x_1,\stackrel{\→}{\mathrm{\ω}})=T(0,l)I(x_0,\stackrel{\→}{\mathrm{\ω}})+$

${\∫}_0^{l}T(s,l)*(S^D\left(s\right)*L_g*\exp (-{\∫}_s^{l_g}\mathrm{\τ}\left(t\right)\mathrm{dt})+S^I\left(s\right)*L_g*\exp (-{\∫}_s^{l_g}{\mathrm{\τ}}_i\left(x\right)\mathrm{dx}))\mathrm{ds}---\left(1\right)$

式中，

指沿着视线方向

在空间坐标x处的光线强度，指视线方向，x为体素的空间三维坐标(x，y，z)。s、l表示沿着视线距离，T(s，l)指沿着视线从x(s)到x(l)过程中的衰减， $T(s,l)=\exp (-{\∫}_s^l\mathrm{\τ}\left(t\right)\mathrm{dt}),$ τ指光线衰减系数；x(s)指沿着视线距离s处的三维坐标。S^D(s)为直接散射的值，S^I(s)为间接散射的值。L_g为光源的光亮，其下角g表示沿着ω′的光线方向，ω′为入射光方向。l_g表示沿着光线方向的距离，其下角g表示沿着ω′的光线方向。τ_i指光线间接衰减系数，T_i(s，l)指沿着光线从x(s)到x(l)过程中的间接衰减， $T_i(s,l)=\exp (-{\∫}_s^l{\mathrm{\τ}}_i\left(t\right)\mathrm{dt}).$

通常，在方向为空间x处的散射光强度

$S(x,\stackrel{\→}{\mathrm{\ω}})=r(x,\stackrel{\→}{\mathrm{\ω}},{\mathrm{\ω}}^{\′})i(x,{\mathrm{\ω}}^{\′})---\left(2\right)$

式中，i(x，ω′)为沿ω′方向入射，到达空间x处的光强度。

为双向反射分布函数(BRDF bi-directional reflection distribution function)。式(2)为通用的明暗化公式。

在一随机分布介质空间V，当光线穿过空间V时，由于光线与介质粒子的相互作用，在光的传播方向上，会发生吸收和散射的消光现象，因此，在x处的

应该由准直辐射亮度和散射辐射亮度组成。而散射辐射亮度由直接散射项与间接散射项两部分组成。设x处的一个粒子的双向反射分布函数

$r(x,\stackrel{\→}{\mathrm{\ω}},{\mathrm{\ω}}^{\′})=a\left(x\right)\mathrm{\τ}\left(x\right)p(\stackrel{\→}{\mathrm{\ω}},{\mathrm{\ω}}^{\′})---\left(3\right)$

式中，α(x)指该粒子的反射率，为相函数，表示散射的方向性。

本质上，结合Blinn-Phong面明暗化模型的光线吸收与发射模型考虑了直接散射的影响。若令

表示直接散射的相函数，则

$p^D(\stackrel{\→}{\mathrm{\ω}},{\mathrm{\ω}}^{\′})=|\stackrel{\→}{N}\·{\mathrm{\ω}}^{\′}|+{(\stackrel{\→}{N}\·\frac{\stackrel{\→}{\mathrm{\ω}}+{\mathrm{\ω}}^{\′}}{|\stackrel{\→}{\mathrm{\ω}}+{\mathrm{\ω}}^{\′}|})}^n---\left(4\right)$

式中，

为r点的梯度向量，也就是r点隐含分界面的法向量。考虑到隐含分界面的正反面方向，对取了绝对值。

若考虑梯度模值高的区域，则用类似Blinn-Phong面明暗化模型的方法，若在梯度模值较低的区域，则主要用阴影实现明暗化计算。另外，考虑用非真实感绘制技术加大隐含分界面的显示效果。例如，给直接散射和间接散射一个权重。

设直接散射为S^D(s)。直观地讲，愈接近隐含分界面，面散射亮度愈大，反之，则愈小。因此，这里的权函数本质上式一个边界检测函数，可以令该体素的权值w^D与梯度模值成正比。则有：

$S^D\left(s\right)=w^D\left[R\left(s\right)((1-w\left(s\right))+w\left(s\right)p^D(\stackrel{\→}{\mathrm{\ω}},{\mathrm{\ω}}^{\′}))\right]---\left(5\right)$

式中，w^D、w(s)为权值。R(s)为沿着视线距离s处面反射颜色值。

另一方面，考虑间接散射的影响。令间接散射

$S^I\left(s\right)=w^{I}R\left(s\right)p^I(\stackrel{\→}{\mathrm{\ω}},{\mathrm{\ω}}^{\′})---\left(6\right)$

式中，w^I为权值。

$p^I(\stackrel{\→}{\mathrm{\ω}},{\mathrm{\ω}}^{\′})\≈\frac{1}{\mathrm{\σ}\left(x\right)}\underset{p=1}{\overset{n}{\mathrm{\Σ}}}f(x,{\stackrel{\→}{\mathrm{\ω}}}_p^{\′},\stackrel{\→}{\mathrm{\ω}})\frac{\mathrm{\Δ}{\mathrm{\Φ}}_{p,i}(x,{\stackrel{\→}{\mathrm{\ω}}}_p^{\′})}{\frac{4}{3}{\mathrm{\πr}}^3}---\left(7\right)$

式中，ΔΦ_p，i为通过间接散射影响该体素的相关的n个光子的亮度和，

为包含这n个的光子的最小球形体积。

实施过程如下：

步骤一、预设置传递函数，将其设为查询表。

T(s，l)指沿着视线从x(s)到x(l)过程中的衰减，定义α为这段距离的不透明度，则α＝1-T(s，l)。同理，定义α_i为这段距离的间接不透明度，则α_i＝1-T_i(s，l)。 $\exp (-{\∫}_s^{l_g}{\mathrm{\τ}}_i\left(x\right)\mathrm{dx})$ 可以近似表示为

$\exp (-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\τ}(i\·\mathrm{\Δx})\·\mathrm{\Δx})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}\exp (-\mathrm{\τ}(i\·\mathrm{\Δx})\·\mathrm{\Δx})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}{t}_i---\left(8\right)$

式中，t_i＝exp(-τ(i·Δx)·Δx)可以看作是沿光线方向上第i个区段的间接透明度。

为快速实现光学模型算法，先预设置Δx＝1的间接不透明度值α_i的查询表。同理，预设置Δx＝1的不透明度值α的查询表。

步骤二、预设置直接散射的权重值

直观地讲，愈接近隐含分界面，直接散射亮度愈大，反之则愈小。因此，这里的权函数本质上是一个边界检测函数，因此，直接散射的权重w^D与梯度模值成正比。另外，若考虑到式(6)的影响，则可设置间接散射的权重w^I与体素值成正比。

步骤三、合成中间图像

本发明基于改进的Shear-Warp体绘制的方法实现。本发明提供的光学模型将在错切变换时合成中间图像的过程中应用。三维空间离散数据场的观察方向常常式由用户任意给定的，因而使得由物体空间到图像空间的变换也是任意的。在本发明的实现过程中，将三维离散数据场变换到一个中间坐标系，在这个中间坐标系中，观察方向与坐标系的一个轴并行。光源的方向是任意的，为叙述和执行方便，不妨将光源设于与视点处于同一侧的位置(是否在同一侧以夹角是否大于90°为依据)。若光源方向在另外一侧，只需要更改光线合成次序则可。

在本发明提供的光学模型TVROM的实现过程中，除了原来需要的缓存外，还需要新开辟三个缓存：当前切片的直接光照缓存、间接光照缓存以及下一个切片的间接光照缓存，缓存大小均为中间图像大小。该半透明体绘制方法也是按由前向后序次，循环对每个切片进行处理，合成中间图像。在处理第k个切片时，过程如下：

①计算直接散射部分。

$C_{\mathrm{post}}^D=C_{\mathrm{pre}}^D+{\mathrm{\α}}_{\mathrm{now}}^D*C_{\mathrm{now}}^D*(1-O_{\mathrm{pre}}^D)*I_{\mathrm{pre}}^D---\left(9\right)$

式中，C_pre ^D为计算直接散射部分之前的中间图像，在当前某一位置的RGB中的R分量的值；C_post ^D为计算直接散射部分之后的中间图像，在当前同一位置的RGB中的R分量的值；O_pre ^D为计算直接散射部分之前的中间图像当前位置的不透明度值；I_pre ^D为计算直接散射部分之前的中间图像当前位置的直接光照累积光亮度；α_now ^D为当前体素重采样不透明度；C_now ^D为当前体素重采样的R分量值。这些变量的上标D表示计算直接散射部分。

$O_{\mathrm{post}}^D=O_{\mathrm{pre}}^D+{\mathrm{\α}}_{\mathrm{now}}^D*(1-O_{\mathrm{pre}}^D)---\left(10\right)$

式中，O_post ^D为计算直接散射部分之后的中间图像当前位置pos的不透明度值。

由于光的合成顺序由后至前，有

$I_{\mathrm{post}}^D=(1-{\mathrm{\α}}_{\mathrm{now}}^D)*I_{\mathrm{pre}}^D---\left(11\right)$

式中，I_post ^D为计算直接散射部分之后的中间图像当前位置的直接光照累积光亮度。

中间图像的RGB中的G、B分量的直接散射部分的值可以用同样的方法求得。

②计算光照累积光亮度。

由光的合成顺序由后至前，有

$I_{\mathrm{post}}^I=I_{\mathrm{pre}}^I*(1-{\mathrm{\α}}_{\mathrm{now}}^I)---\left(12\right)$

式中，I_pre ^I为计算间接散射部分之前的中间图像当前位置的间接光照累积光亮度；I_post ^I为计算直接散射部分后的中间图像当前位置的间接光照累积光亮度；α_now ^I为当前体素重采样的间接不透明度值。这些变量的上标I表示计算间接散射部分。

③计算间接散射部分。

$C_{\mathrm{post}}^I=C_{\mathrm{pre}}^I+{\mathrm{\α}}_{\mathrm{now}}^I*C_{\mathrm{now}}^I*(1-O_{\mathrm{pre}}^D)*I_{\mathrm{post}}^I*I_a^I---\left(13\right)$

式中，C_pre ^I为添加间接散射部分前的中间图像当前位置的RGB中的R分量的值；C_post ^I为添加间接散射部分后的中间图像当前位置的RGB中的R分量的值；C_now ^I为当前体素在计算间接散射部分时重采样的R分量值。

式(13)中，I_a ^I取六个像素点(下一个切片的间接光照缓存在pos位置的像素点、当前切片的间接光照缓存在pos位置的像素点以及其相邻四个像素点)的平均值。这里，I_a ^I实质上是式(7)中 $\underset{p=1}{\overset{n}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{\mathrm{\Φ}}_{p,i}(x,{\stackrel{\→}{\mathrm{\ω}}}_p^{\′})}{\frac{4}{3}{\mathrm{\πr}}^3}$ 的近似值。中间图像的RGB中的G、B分量的间接散射部分的值同理求得。

④将下一个切片的间接光照缓存数据存放于当前切片的间接光照缓存，并将下一个切片的间接光照缓存都设置为值1。

⑤处理下一个切片。

实施效果：

处理结果如图2所示。三维数据场的半透明体绘制速度：一帧2.87秒。

图2a为用本发明提供的用于半透明体绘制的光学模型TVROM及实现方法对头部进行半透明体绘制的结果。图2b为光学模型TVROM的直接散射部分结果。图2c为光学模型TVROM的间接散射部分。图2d为用本发明提供的用于半透明体绘制的光学模型TVROM及实现方法对脑颅进行半透明体绘制的结果。由图2可见，本发明的方法能充分显示重建对象的内部隐含分界面及内部细节的详细信息。

实施例2

用本发明——一种融合解剖与功能成像信息数据场的半透明体绘制方法对源于怀疑肺癌的某女性病人的CT、PET序列断层图像进行处理。CT、PET的序列断层图像各266张，用GE公司的Discovery ST PET-CT对此怀疑肺癌的某女性病人进行扫描得到，CT原始分辨率为512×512，PET原始分辨率为128×128。

如图1为对照融合解剖与功能成像信息数据场的半透明体绘制方法的流程图，实施过程如下：

步骤一、形成融合解剖与功能成像信息的三维数据场。在此过程中，基于小波变换，以区域标准差与区域能量相结合的融合规则进行解剖与功能成像信息融合。

在图像的小波变换域内，对已经配准的解剖与功能成像进行信息融合。

在临床应用中，医生主要关心病变的位置特征。解剖图像以丰富的纹理清晰地描绘了人体组织的形态结构。图像融合的目的是为功能图像描绘的病变区域提供定位参考。因此，选用的活动性水平测度应能反映纹理模式。

纹理是一种视觉感知的图像局部特征的综合。纹理特征提取的目的是获得每一个像素点的特征向量。纹理特征描述一般使用能量、方差、加权平均、标准差和最大概率等。图像某一位置的纹理特征与这一位置周围的灰度变化规律密切相关。图像平滑区域包含的象素灰度彼此接近，而粗糙区域的象素灰度有较大变化。因而，可以采用区域标准差作为系数的活动性水平测度。此外，局部能量可以反映图像信号变化的绝对强度，而信号变化强度大的点反映了图像的显著特征，因此，依据局部能量可以给出图像特征的统一描述，同时反映图像的清晰度。基于上述原因，本发明采用区域标准差和能量作为系数的活动性水平测度。

a、基于区域标准差的活动性水平测度A_X为：

$A_X(i,j)=\sqrt{\underset{s\∈S,t\∈T}{\mathrm{\Σ}}\mathrm{\ω}(s,t){[D_X(i+s,j+t,k,l)-D_X(i,j)]}^2}---\left(14\right)$

式中，X＝CT，PET。D_CT(i，j)、D_PET(i，j)分别表示CT和PET图像中窗口的小波系数平均值，(i，j)表示窗口中心所处位置，k为分解层数(k＝1，2，3)，l为频带(水平、垂直或对角)，(s，t)表示以(i，j)为中心的窗中的位置，ω(s，t)为以(i，j)为中心的窗中(s，t)处系数的权重，离窗中心越远权重越小，且 $\underset{s\∈S,t\∈T}{\mathrm{\Σ}}\mathrm{\ω}(s,t)=1,$ S、T分别为窗的水平、垂直指标，一般取3×3或5×5。。

令δ_CT，δ_PET表示基于局域标准差的活动性水平测度分配给CT、PET的权重，

$\left\{\begin{array}{c}{\mathrm{\δ}}_{\mathrm{CT}}=\frac{{\left[A_{\mathrm{CT}}(i,j)\right]}^{\mathrm{\α}}}{{\left[A_{\mathrm{CT}}(i,j)\right]}^{\mathrm{\α}}+{\left[A_{\mathrm{PET}}(i,j)\right]}^{\mathrm{\α}}}\\ {\mathrm{\δ}}_{\mathrm{PET}}=\frac{{\left[A_{\mathrm{PET}}(x,y)\right]}^{\mathrm{\α}}}{{\left[A_{\mathrm{CT}}(i,j)\right]}^{\mathrm{\α}}+{\left[A_{\mathrm{PET}}(i,j)\right]}^{\mathrm{\α}}}\end{array}\right.---\left(15\right)$

式中，α为可调节参数。当α＞0时，活动性水平测度大的所占权重大；反之，活动性水平测度小的所占权重大。

b、基于区域能量的活动性水平测度B_X为：

$B_X(i,j)=\underset{s\∈S,t\∈T}{\mathrm{\Σ}}\mathrm{\ω}(s,t)D_X^2(i+s,j+t,k,l)---\left(16\right)$

令ε_CT，ε_PET为基于区域能量活动性水平测度分配给CT、PET的权重，

$\left\{\begin{array}{c}{\mathrm{\ϵ}}_{\mathrm{CT}}=\frac{B_{\mathrm{CT}}(i,j)}{B_{\mathrm{CT}}(i,j)+B_{\mathrm{PET}}(i,j)}\\ {\mathrm{\ϵ}}_{\mathrm{PET}}=\frac{B_{\mathrm{PET}}(i,j)}{B_{\mathrm{CT}}(i,j)+B_{\mathrm{PET}}(i,j)}\end{array}\right.---\left(17\right)$

c、结合局部标准差活动性水平测度和区域能量活动性水平测度，得到融合后图像的小波系数D_F为：

D_F(i，j)＝[δ_CTD_CT(i，j)+δ_PETD_PET(i，j)]×λ+[ε_CTD_CT(i，j)+ε_PETD_PET(i，j)]×μ(18)

式中，λ、μ为可调节参数，λ+μ＝1，随着μ增大图像亮度增强，随着λ增大图像边缘增强，适当调整λ、μ，可以尽量避免丧失边缘信息，同时消除图像模糊。

步骤一的实现过程有：

①对CT、PET图像进行配准。

②选取双正交滤波器组的Daubechies 9/7双正交小波分别对配准后的CT和PET进行三层双正交小波分解。

③分别计算CT和PET图像中窗口的小波系数平均值D_CT(i，j)和D_PET(i，j)。

$\left\{\begin{array}{c}{D}_X(i,j)=\underset{s\∈S,t\∈T}{\mathrm{\Σ}}\mathrm{\ω}(s,t)D_X(i+s,j+t,k,l)\\ X=\mathrm{CT},\mathrm{PET}\end{array}\right.---\left(19\right)$

④在图像的小波变换域内，对图像进行融合。本发明采用区域标准差和能量作为系数的活动性水平测度。基于区域标准差的活动性水平测度A_X如式(14)，基于局域标准差的活动性水平测度分配给CT、PET的权重δ_CT、δ_PET如式(15)，这里，α取值为1.8。基于区域能量的活动性水平测度B_X如式(16)，基于区域能量活动性水平测度分配给CT、PET的权重ε_CT，ε_PET如式(17)。结合局部标准差活动性水平测度和区域能量活动性水平测度，得到融合后图像的小波系数D_F如式(18)。

⑤对CT图像和PET图像经小波变换之后的逼近系数C_J ^CT和C_J ^PET进行处理。由于图像模糊表示其细节信息(或高频信息)丢失较多，相比之下，其整体信息(或低频信息)保持较好，因此两幅图像经小波分解后其逼近系数之间的差异要远小于小波系数之间的差异，故融合之后的逼近系数可

由如下式确定：

${\hat{C}}_J^F=(C_J^{\mathrm{CT}}+C_J^{\mathrm{PET}})/2---\left(20\right)$

⑥利用以上得到的全部小波系数D_F以及

中的逼近系数进行小波反变换可得融合图像F。

步骤二、用实现的用于半透明体绘制的光学模型TVROM及其实现方法对融合解剖与功能成像信息数据场进行半透明体绘制。光学模型TVROM的半透明体绘制方法如图1的虚线框图内部分所示。

对于在普通硬件条件下实现光学模型TVROM的实时实现问题的解决。本发明解决的方法主要利用不透明度传递函数的相关性，以及min-max Octree数据结构，实现对体数据的快速分类编码，并在此基础上，快速再建游程编码RLE数据结构，以满足在不透明度传递函数和观察角度频繁改变时，也能实现快速重建。

实现过程如下：

①快速再建游程编码(RLE)数据结构。

首先，在交互式操作中，当不透明度传递函数改变时，在对图像进行重采样前，可以利用不透明度传递函数的相关性，以及已存在的min-max Octree分类数据结构，对体数据进行加速分类。即，利用不透明度传递函数的相关性，对原有的min-max Octree进行再利用，不需要更改图中所有结点的数据属性，而仅仅更改某些结点的数据属性，从而提高效率。

接着，利用再分类后的分类编码数据结构，快速再建游程编码(RLE)数据结构。RLE数据结构是以扫描行为基本单位的，当重新建立RLE数据结构时，可以不再重新完全地遍历整个原始数据，也无需对原RLE数据结构完全进行更改，而只是对需要的部分扫描行(依据分类编码数据结构判断)进行遍历处理即可。

而后，利用新的RLE数据结构实现任意角度的三维重建。

②在错切空间里，利用光学模型TVROM，合并中间图像。该部分的实现过程与实施例1的实现过程相同。包括a)预设置传递函数，将其设为查询表；b)预设置直接散射的权重值；c)合成中间图像，在处理第k个切片时，主要包括计算直接散射部分、计算光照累积光亮度、计算间接散射部分几个部分。

③在变形空间里，用纹理映射绘制最终图像。

实施效果：

处理结果如图3a、图3b、图4a～图4d所示。融合解剖与功能成像信息数据场的半透明体绘制速度：一帧1.97秒。

图3a、图3b为半透明三维显示的靶区与病人内部结构的结果。图3a为从正面观察的结果，图3b为从背面观察的结果。由此可以确定靶区的空间位置。

图4a为体绘制结果，图4b为图4a白线所指向的横断面(融合解剖与功能成像信息数据场的横断剖面)，图4c为融合解剖与功能成像信息数据场的横断剖面所对应的CT图像，图4d为融合解剖与功能成像信息数据场的横断剖面所对应的PET图像。图4b、图4c、图4d十字的交点指向病人的同一位置。从图4c可以看出，右肺下叶背段有一结节影；图4d的左下部有一亮斑，显示有一显像药物的高吸收区，但病灶区的解剖细节不清晰；图4b则清楚地显示了CT图像中的节影区与PET图像中的癌渗透区的对应关系。由此可见，采用基于小波变换的融合方法，以区域标准差与区域能量相结合的融合规则对CT、PET图像进行融合处理后，有效地保留了多模态图像的边缘和纹理特征，避免了融合图像出现模糊现象，为PET图像描述的病变区域的准确定位提供了解剖信息。

由三维显示的病人三维图像(融合解剖与功能成像信息数据场)可见：右肺下叶背段见一结节状软组织密度影，大小约2.4×2.5cm，出现异常放射性浓聚。由此，得出结论：右肺下叶背段占位(2.4×2.5cm)，糖代谢增高——肺癌。

Claims (4)

1.一种融合解剖与功能成像信息数据场的半透明体绘制方法，其特征在于包括如下步骤：

(1)形成融合了解剖与功能成像信息的三维数据场：基于小波变换，以区域标准差与区域能量相结合的融合方法进行解剖与功能成像的信息融合；在小波域里，采用区域标准差和能量作为小波系数的活动性水平测度；

(2)应用半透明体绘制的光学模型TVROM对融合解剖与功能成像信息数据场进行半透明体绘制：

a、再建游程编码数据结构：在交互式操作中，当不透明度传递函数改变时，在对图像进行重采样前，利用不透明度传递函数的相关性，和已存在的min-max Octree分类数据结构，对体数据进行分类；随后，利用再分类后的分类编码数据结构，再建游程编码数据结构；

b、在错切空间里，利用光学模型TVROM，合并中间图像；

c、在变形空间里，用纹理映射绘制最终图像。

2.根据权利要求1所述的融合解剖与功能成像信息数据场的半透明体绘制方法，其特征在于，所述步骤(1)基于小波变换，以区域标准差与区域能量相结合的融合方法进行解剖与功能成像的信息融合包括如下步骤：

①对CT、PET图像进行配准；

②选取双正交滤波器组的Daubechies 9/7双正交小波分别对配准后的CT和PET进行三层双正交小波分解；

③分别计算CT和PET图像中窗口的小波系数平均值；

④在图像的小波变换域内，对图像进行融合；

⑤对CT图像和PET图像经小波变换之后的逼近系数进行求和后除以2的处理；

⑥利用小波系数中的逼近系数进行小波反变换可得融合图像。

3.根据权利要求1所述的融合解剖与功能成像信息数据场的半透明体绘制方法，其特征在于，所述在错切空间里，利用光学模型TVROM，合并中间图像包括如下步骤：

a)预设置传递函数，将其设为查询表；

b)预设置直接散射的权重值；

c)合成中间图像：由前向后序次，循环对每个切片进行处理，合成中间图像，在处理第k个切片时，依次计算直接散射部分、计算光照累积光亮度、计算间接散射部分；将下一个切片的间接光照缓存数据存放于当前切片的间接光照缓存，并将下一个切片的间接光照缓存都设置为值1；处理下一个切片。

4.根据权利要求1所述的融合解剖与功能成像信息数据场的半透明体绘制方法，其特征在于，所述用于半透明体绘制的光学模型TVROM为：

$I=(x_1,\stackrel{\→}{\mathrm{\ω}})=T(0,l)I(x_0,\stackrel{\→}{\mathrm{\ω}})+$

${\∫}_0^{l}T(s,l)*(S^D\left(s\right)*L_g*\exp \left({\∫}_s^{l_g}\mathrm{\τ}\left(t\right)\mathrm{dt}\right)+S^l\left(s\right)*L_g*\exp (-{\∫}_s^{l_g}{\mathrm{\τ}}_i\left(x\right)\mathrm{dx}))\mathrm{ds}$

式中，

指沿着视线方向

在空间坐标x处的光线强度，

指视线方向，x为体素的空间三维坐标(x，y，z)，s、l表示沿着视线距离，T(s，l)指沿着视线从x(s)到x(l)过程中的衰减，

τ指光线衰减系数：x(s)指沿着视线距离s处的三维坐标；S^D(s)为直接散射的值，S^l(s)为间接散射的值；L_g为光源的光亮，其下角g表示沿着ω′的光线方向，ω′为入射光方向；l_g表示沿着光线方向的距离，其下角g表示沿着ω′的光线方向；τ_i指光线间接衰减系数，T_i(s，l)指沿着光线从x(s)到x(l)过程中的间接衰减，

$T_i(s,l)=\exp (-{\∫}_s^l{\mathrm{\τ}}_i\left(t\right)\mathrm{dt}).$

Images