CN100498839C - 一种多模态医学体数据三维可视化方法 - Google Patents

Abstract

本发明涉及一种多模态医学体数据的三维可视化方法。现有技术中的配准多基于二维-二维层面或二维-三维层面，同时不能完整的实现医学体数据配准和融合显示过程。本发明主要包括基于多分辨率的标准化互信息法配准和单点多模态直接体绘制法融合显示。其中基于多分辨率的标准化互信息法配准包括坐标变换、定位判据、多分辨率优化三个步骤；单点多模态直接体绘制法融合显示包括传递函数定义、光照模型的计算、图像合成、多模态数据的单点多模表示、融合显示五个步骤。采用本发明方法可以完整的实现医学体数据配准和融合显示过程，是真正三维-三维意义上的多模态可视化。

Description

一种多模态医学体数据三维可视化方法

技术领域

本发明属于生物医学工程领域，涉及一种多模态医学体数据的三维可视化方法。

背景技术

医学图像数据可分为两大类：一类是X线计算机断层(CT)和核磁共振(MR)等解剖图像，适合于显示人体的解剖细节，具有较高的空间分辨率；另一类是正电子发射断层(PET)、单光子发射断层(SPECT)等功能图像，适合于反映人体组织的功能改变和代谢异常等信息，其空间分辨率差。由于成像原理不同所造成的图像信息局限性，使得单独使用一种模态数据的效果并不理想。通过多模态可视化可以把多个扫描源或不同时间获取的同一病人或多个病人间的数据进行适当综合，将功能信息、新陈代谢及解剖信息等融合在一起为医生提供与解剖位置相关的功能变异，实现解剖图像与功能图像的信息互补，有利于发现新的信息，弥补单模态图像信息不完整、部分信息不准确或不确定引起的缺陷，使临床诊断和治疗、放疗的定位和计划设计、外科手术和疗效评估等更加准确完善。多模态可视化的关键在于多模态数据间的配准和融合显示。

配准是为了多模态体数据间空间上的对齐，主要有基于图像特征匹配和基于图像灰度信息相似性两类方法，目前已申请的配准方法如：刚体变换下基于轮廓多源图像配准方法(200310108731.X)主要针对多源二维图像采用轮廓法实现图像配准，属于基于图像特征的方法；一种基于互信息敏感区域的多模态医学图像配准方法(200310122434.0)提出一种基于互信息敏感区域的图像配准方法，该法在互信息敏感区域的选取上需要经验知识，因此并不是全自动的方法；图像配准方法(200410015656.7)提出基于快速相关法和最大化互信息法来进行图像配准，在该法中坐标变换是基于刚体变换的，即没有考虑弹性缩放问题，而互信法存在重叠区域易变性问题。

目前已申请的有关医学体数据显示方法的发明有：分布式心脏图像多维重建与交互可视化方法(200510020318.7)提出一种采用分布式结构实现旋转扫描二维图像快速重建三维数据场的交互显示的方法；具有封闭轮廓多平面重定格式的体绘制数据的可视化(200510065505.7)提出一种用于可视化与查看矢量相关的体绘制图像的方法；基于点重建的超大规模医学影像三维可视化方法(200310121173.0)提出一种基于点重建的超大规模医学影像三维可视化方法；脑外科手术导航系统中解剖图谱的三维可视化应用方法(03116689.X)提出一种脑外科手术导航系统中解剖图谱的三维可视化应用方法。以上这些发明都是针对单模态医学体数据的显示方法，没有解决多种模态的医学体数据的融合显示问题。容积-容积融合的可视化(03805385.3)提出利用多个容积的方法来融合显示图像的方法，该发明在三维体数据的融合显示采用提取二维层面，其实质是基于二维-二维的融合显示方式，没有解决真正基于三维立体的多模态数据融合显示问题。

以上申请的发明专利表明，在医学多模态体数据可视化领域，配准与融合显示的研究工作都是各自进行的，即专门从事配准研究或对已配准的多模态数据进行可视化，且配准多基于二维-二维层面或二维-三维层面。

发明内容

本发明针对现有技术的不足，从信息论的角度出发，提供一种基于多分辨率的标准化互信息的方法来对三维-三维间的医学体数据进行高效、全自动的弹性配准，并提供一种新的基于单点多模态的直接体绘制法来进行多模态医学体数据的融合显示的方法，从而完整的实现医学体数据配准和融合显示过程，是真正三维-三维意义上的多模态可视化。

本发明主要包括两大步骤：基于多分辨率的标准化互信息法配准和单点多模态直接体绘制法融合显示。

基于多分辨率的标准化互信息法配准包括坐标变换、定位判据、多分辨率优化三个步骤：

(1)坐标变换方法

设要进行多模态可视化的两个体数据分别称为参考体和浮动体，对两个待配准体数据分别定义一个物空间坐标系，其中X轴为沿行扫描方向，Y轴为沿列扫描方向，Z轴则沿从颅顶到颅底的方向。传统的刚体变换包括平移、旋转矩阵，分别以T、R表示，本发明通过增加比例因子S，解决弹性配准问题，假设A是参考体到浮动体的坐标变换函数，则A可表示成：

A＝T(t_x，t_y，t_z)*R(φ_x，φ_y，φ_z)*S(s_x，s_y，s_z)

t_x、t_y、t_z分别表示x、y、z轴方向的平移，φ_x、φ_y、φ_z分别表示绕x、y、z轴的旋转角度，s_x、s_y、s_z分别为x、y、z轴方向的缩放比例，利用坐标变换函数A将从浮动体的网格点提取的样本进行坐标变换，如果坐标变换后的坐标未映射到参考体的网格点上，这时通过三次线性插值法得到整数坐标。

(2)定位判据方法

在坐标变换过程中寻找参数：

${\stackrel{\→}{\mathrm{\α}}}^{*}=(t_x,t_y,t_y,{\mathrm{\φ}}_x,{\mathrm{\φ}}_y,{\mathrm{\φ}}_z,s_x,s_y,s_z),$

使之满足下式：

${\stackrel{\→}{\mathrm{\α}}}^{*}=\arg \underset{\stackrel{\→}{\mathrm{\α}}}{\max }{\mathrm{NI}}_{\mathrm{FR}}\left(\stackrel{\→}{\mathrm{\α}}\right)---\left(1\right)$

此时，两个体数据对应点间的标准化互信息NI_FR最大。通过计算两个体数据重叠部分体素的联合灰度直方图来估计浮动体的熵H_F、参考体的熵H_R及两者的联合熵H_FR，从而计算当前配准位置的标准化互信息NI_FR。

标准化互信息NI_FR计算公式为：

NI_FR＝(H_F+H_R)/H_FR          (2)

$H_F=-\underset{f=1}{\overset{n_F-1}{\mathrm{\Σ}}}p\left(f\right)\·{\log }_{2}p\left(f\right)---\left(3\right)$

$H_R=-\underset{r=1}{\overset{n_R-1}{\mathrm{\Σ}}}p\left(r\right)\·{\log }_{2}p\left(r\right)---\left(4\right)$

$H_{\mathrm{FR}}=-\underset{f=1}{\overset{(n_F-1)}{\mathrm{\Σ}}}\underset{r=1}{\overset{(n_R-1)}{\mathrm{\Σ}}}p(f,r)\·{\log }_{2}p(f,r)---\left(5\right)$

p(f)、p(r)和p(f，r)分别为浮动体和参考体的概率密度函数及它们的联合概率密度函数，n_F、n_R是参考体和浮动体数据灰度级别的个数，一般都标准值为255。

(3)多分辨率优化方法

由于计算配准判据的时间与采样点数成比例，通过欠采样浮动图像在低分辨率下优化可极大的提高速度。将图像配准过程分解为由粗到细：先将原始图像进行子采样，降低其分辨率，按照(1)和(2)的坐标变换和定位判据的方法在低分辨率下进行配准，得到一个最优的变换数据；然后以此变换参数为初始点，按照(1)和(2)的坐标变换和定位判据的方法进行更高分辨率的配准；由于前次配准的结果已经接近于最终结果，因此，后一次配准可以减小变换参数的搜索范围，从而达到减少迭代次数提高配准速度的效果。另外引入多分辨率搜索方法可以避免陷入局部极值。

配准后的多模态数据采用单点多模态直接体绘制法进行融合显示，包括传递函数定义、光照模型的计算、图像合成、多模态数据的单点多模表示、融合显示五个步骤。

(1)传递函数定义的方法

采集到的医学体数据多为灰度数据，直接体绘制中对每个体素赋予颜色值和不透明值，不同的人体组织器官具有不同的灰度分布，本发明中的传递函数则根据灰度直方图中不同的灰度分布赋予颜色值和不透明度值。

(2)光照模型的计算方法

在直接体绘制中，环境光和周围的点光源都会对最终的显示效果产生影响，在医学体数据的可视化中，光照的逼真问题并不是最重要的因素，为了提高绘制速度，本发明将Phone模型加以简化，采用单个点光源，简化公式为：

I＝k_aI_a+fI_l[k_d(NL)+k_s(NH)]       (6)

I为局部光亮度；k_a、k_d、k_s分别为环境、漫反射及镜面反射系数；I_a为环境光亮度分量；I₁、f为点光源光强及衰减因子；N为单位法向量；L为点光源的单位入射向量；H为点光源的单位入射向量与单位视线向量的平均值。

采用基于灰度梯度的中心差分法计算中间体素的单位法向量N，对于边界点处的法向量计算则采用向前或向后差分法。中心差分法可表示成：

Grad_x＝[f(x_i+1，y_j，z_k)-f(x_i-1，y_j，z_k)]/2Δx

Grad_y＝[f(x_i，y_j+1，z_k)-f(x_i，y_j-1，z_k)]/2Δy　  (7)

Grad_z＝[f(x_i，y_j，z_k+1)-f(x_i，y_j，z_k-1)]/2Δz

Grad_x表示x轴方向的灰度梯度，Grad_y表示y轴方向的灰度梯度，Grad_z表示z轴方向的灰度梯度。

(3)图像合成的方法

本发明采用由前向后的图像合成方法，在合成过程中，不透明度β必然逐步增大，当β值趋近于1时，说明该像素点的图像已接近于完全不透明，后面的体元不会再对该像素点的图像有所贡献，因而可以不必再计算了，从而可以省去无效的计算，提高计算速度。对于多个(2)中光照模型的计算结果由前向后进行图像合成，具体公式可表示成：

C_outβ_out＝C_inβ_in+C_nowβ_now(1-β_in)    (9)

β_out＝β_in+β_now(1-β_in)                           (10)

C_out、β_out为经过第i个体元后的颜色值和不透明度值，C_now、β_now为第i个体元的颜色值和不透明度值，C_in、β_in为进入第i个体元的颜色值和不透明度值。

(4)多模态数据的单点多模表示的方法

为了在一个像素点上同时显示多模态体数据信息，我们采用视觉颜色HLS模型表示(3)中图像合成的结果。人眼的颜色知觉主要包括三个要素：即色调(Hue)、饱和度(Saturation)和亮度(Luminance)，色调指光的颜色，不同波长的光呈现不同的颜色，具有不同的色调。饱和度指颜色的深浅或浓淡程度，其深浅与颜色中加入白色的比例有关。亮度就是人眼能感觉到的光的明暗程度，光波的能量越大，亮度就越大。

视觉颜色模型以人眼的视觉特征为基础，用色调、饱和度和亮度来表示各个多模态数据的颜色属性。如：HLS模型中的亮度表示参考体属性，而色调和饱和度则表示浮动体的颜色属性。通常以高空间分辨率的模态数据来表示HLS模型中的亮度信息。

(5)融合显示的方法

将HLS模型转化为RGB模型，并加以显示。

采用本发明方法可以完整的实现医学体数据配准和融合显示过程，是真正三维-三维意义上的多模态可视化，同时本发明具有以下特点：

1、用标准化互信息作为多模态医学体数据的配准测度，不需要对两种模态数据强度间关系的性质作任何假设，也不需要对图像作预分割或任何预处理，因此用这种方法来进行多模态体数据的配准，能提高自动化程度和配准精度；同时在判据中考虑了坐标变换中的弹性缩放，适合任何不同模态体数据的配准；且标准化互信息法可解决互信息法难以克服的重叠区域易变性问题。

2、通过多分辨率优化配准过程可极大的提高配准的速度，同时避免陷入局部极值。

3、单点单模态的表示仅仅提供多个单模态可视化结果的数学叠加，表示效果受限，单点多模态直接体绘制法可使每个像素点同时表示多种模态数据属性，能更有效的表现多模态融合信息。

具体实施方式

我们对CT和PET两种模态数据的三维可视化进行实例说明。令CT体数据为参考体，PET体数据为浮动体，各自的坐标为物空间坐标。首先我们设置多分辨率的粗配准为8个像素长度单位，对CT和PET数据进行坐标变换，并计算坐标变换后的标准化互信息，找出标准化互信息最大时的坐标参数，将该最优点作为起始点，缩小坐标变换的范围，同时将坐标变换单位长度减半，即提高分辨率，继续前面的坐标变换及标准化互信息判据过程，直至最优，在新的最优基础上，重复前面的步骤，由粗配准到精细配准步步逼近，直至单位长度为1，得到此时的坐标变换参数即为配准参数。

得到配准参数后，从像空间平面发出射线，确定采样点在射线上的分布间隔，采样点像空间坐标变换至CT物空间坐标，读取CT的相应灰度值，同时利用配准参数得到对应PET物空间中点的灰度值；CT和PET各自对采样点灰度样本分类，赋颜色值和不透明度值，计算光照效应，利用公式(9)(10)由前向后合成图像，把空间分辨率较高的CT颜色值映射到HSL模型中的亮度值，而PET颜色值映射到HSL中的色调和饱和度，最后将HSL模型转换为RGB模型得到融合显示结果。

Claims (1)

1、一种多模态医学体数据三维可视化方法，包括配准和融合显示，其特征在于配准采用基于多分辨率的标准化互信息法，具体包括坐标变换、定位判据、多分辨率优化，其中

(1)坐标变换方法

设要进行多模态可视化的两个体数据分别称为参考体和浮动体，对两个待配准体数据分别定义一个物空间坐标系，其中X轴为沿行扫描方向，Y轴为沿列扫描方向，Z轴则沿从颅顶到颅底的方向，传统的刚体变换包括平移、旋转矩阵，分别以T、R表示，通过增加比例因子S，解决弹性配准问题，假设A是参考体到浮动体的坐标变换函数，则A可表示成：

A＝T(t_x，t_y，t_z)*R(φ_x，φ_y，φ_z)*S(s_x，s_y，s_z)

t_x、t_y、t_z分别表示x、y、z轴方向的平移，

分别表示绕x、y、z轴的旋转角度，s_x、s_y、s_z分别为x、y、z轴方向的缩放比例，利用坐标变换函数A将从浮动体的网格点提取的样本进行坐标变换，如果坐标变换后的坐标未映射到参考体的网格点上，通过三次线性插值法得到整数坐标；

(2)定位判据方法

在坐标变换过程中寻找参数：

${\stackrel{\→}{\mathrm{\α}}}^{*}=(t_x,t_y,t_y,{\mathrm{\φ}}_x,{\mathrm{\φ}}_y,{\mathrm{\φ}}_z,s_x,s_y,s_z),$

使之满足下式：

${\stackrel{\→}{\mathrm{\α}}}^{*}=\arg \underset{\stackrel{\→}{\mathrm{\α}}}{\max }{\mathrm{NI}}_{\mathrm{FR}}\left(\stackrel{\→}{\mathrm{\α}}\right)$

此时，两个体数据对应点间的标准化互信息NI_FR最大，通过计算两个体数据重叠部分体素的联合灰度直方图来估计浮动体的熵H_F、参考体的熵H_R及两者的联合熵H_FR，从而计算当前配准位置的标准化互信息NI_FR，计算公式为：

NI_FR＝(H_F+H_R)/H_FR

$H_F=-\underset{f=1}{\overset{n_F-1}{\mathrm{\Σ}}}p\left(f\right)\·{\log }_{2}p\left(f\right)$

$H_R=-\underset{r=1}{\overset{n_R-1}{\mathrm{\Σ}}}p\left(r\right)\·{\log }_{2}p\left(r\right)$

$H_{\mathrm{FR}}=-\underset{f=1}{\overset{(n_F-1)}{\mathrm{\Σ}}}\underset{r=1}{\overset{(n_R-1)}{\mathrm{\Σ}}}p(f,r)\·{\log }_{2}p(f,r)$

p(f)、p(r)和p(f，r)分别为浮动体和参考体的概率密度函数及它们的联合概率密度函数，n_F、n_R是参考体和浮动体数据灰度级别的个数，采用标准值255；

(3)多分辨率优化方法

将图像配准过程分解为由粗到细：先将原始图像进行子采样，降低其分辨率，在低分辨率下按照(1)和(2)的坐标变换和定位判据的方法进行配准，得到一个最优的变换数据；然后以此变换参数为初始点，按照(1)和(2)的坐标变换和定位判据的方法进行更高分辨率的配准；在前次配准的基础上减小变换参数的搜索范围；

配准后的多模态数据采用单点多模态直接体绘制法进行融合显示，具体包括传递函数定义、光照模型的计算、图像合成、多模态数据的单点多模表示、融合显示，其中

(1)传递函数定义的方法

将采集到的医学体数据根据人体不同组织和器官的分类对每个体素赋予颜色值和不透明值，其中分类采用灰度直方图的方法；

(2)光照模型的计算方法

对Phone模型加以简化，采用单个点光源，公式为：

I＝k_aI_a+fI_l[k_d(NL)+k_s(NH)]

I为局部光亮度；k_a、k_d、k_s分别为环境、漫反射及镜面反射系数；I_a为环境光亮度分量；I₁、f为点光源光强及衰减因子；N为单位法向量；L为点光源的单位入射向量；H为点光源的单位入射向量与单位视线向量的平均值，

采用基于灰度梯度的中心差分法计算中间体素的单位法向量N，表示成：

Grad_x表示x轴方向的灰度梯度，Grad_y表示y轴方向的灰度梯度，Grad_z表示z轴方向的灰度梯度；

(3)图像合成的方法

对于多个(2)中光照模型的计算结果由前向后进行图像合成，具体公式表示成：

C_outβ_out＝C_inβ_in+C_nowβ_now(1-β_in)

β_out＝β_in+β_now(1-β_in)

C_out、β_out为经过第i个体元后的颜色值和不透明度值，C_now、β_now为第i个体元的颜色值和不透明度值，C_in、β_in为进入第i个体元的颜色值和不透明度值；

(4)多模态数据的单点多模表示的方法

采用视觉颜色HLS模型表示(3)中图像合成的结果，该模型根据人眼的颜色知觉三个要素，即色调、饱和度和亮度来表示各个多模态数据的颜色属性；

(5)融合显示的方法

将HLS模型转化为RGB模型，并加以显示。