CN101853518B - 基于各向异性体数据的错切变形体绘制方法 - Google Patents

Abstract

基于各向异性体数据的错切变形体绘制方法，属于图像处理技术领域。包括构造三维体数据场、错切变形分解、重采样、合成中间图像、变形得到最终图像等步骤；本发明主要对原始错切变形算法的坐标变换、重采样、不透明度校正等步骤做了改进，可以在图像数据构造成逻辑的三维体数据场之后，不需要通过断层插值操作使其各向同性，而是直接基于各向异性数据场做体绘制操作。不仅减少了预处理时间，也提高了用于绘制的图像数据精度。另外在实施过程中，采用了可变重采样的方式，鼠标交互时降低采样率提高绘制速度，绘制最终图像时提高采样率增强图像质量，同时满足了绘制速度和图像精度的要求。

Description

基于各向异性体数据的错切变形体绘制方法

技术领域

本发明属于图像处理技术领域，涉及可视化三维图像，尤其是针对各向异性体数据直接进行错切变形的体绘制方法。

背景技术

医学图像在医生诊断中的辅助作用越来越明显，将二维切片图像序列转变为具有直观立体效果的三维图像，可以大大提高医疗诊断和治疗规划的准确性和科学性。传统的三维可视化方法分为面绘制和体绘制两种，而体绘制又主要包括：

一、光线投射法：在设定的观察方向下，从屏幕的每一个像素点发出一条射线并穿过三维数据场，沿着射线做重采样和图像合成，最终形成该像素的图像值。光线投射方法生成的图像质量高，但由于需要在三维空间中做大量的重采样操作，绘制速度比较慢，难以满足实时绘制的要求。

二、错切变形法：将三维离散数据场的投影变换分解为三维空间的错切变换和二维空间的变形操作，重采样过程由三维降为二维，大大减少了计算量，使得三维数据场的体绘制可以在图形工作站上以接近实时的速度实现，是一个颇具应用前景的体绘制方法。

三、足迹表法：逐层、逐行、逐个的计算每一个数据点对屏幕像素的贡献并加以合成，形成最终的图像。该方法的一个缺点也是计算量大，难以满足实时绘制的要求。

四、体元投射法：将三维数据场看作多个大小和形状相同的体元，计算每一个体元在屏幕上的投影。由于每一个体元在图像屏幕上的平行投影在几何上一致，利用此空间相关性可以减少计算量，适用于体元比较大的数据场。

错切变形法因其接近实时的绘制速度得到了广泛的应用。传统的错切变形法是基于各向同性体数据场，即二维切片图像序列中相邻两切片图像之间的距离等于二维切片图像内相邻像素点之间的距离。对于二维切片图像序列中相邻两切片图像之间的距离大于二维切片图像内相邻像素点之间的距离的体数据场，传统的错切变形法在进行体绘制之前，需要通过切片间插值技术实现各向同性后再进行绘制。而图像插值是一个具有很大任意性的问题，不合适的插值算法直接影响插值后的图像精度，并影响最后绘制的图像质量。随着技术的不断更新，新一代CT和MRI设备分辨率不断提高，切片间距也不断变小，甚至已接近像素间距。在此背景下，直接基于各向异性体数据进行体绘制已经成为一种技术需求。

发明内容

本发明提出一种基于各向异性体数据的错切变形体绘制方法，具有图形绘制运算量小、速度快和图形绘制质量高的特点。该方法采用错切变形法直接针对各项异性体数据进行体绘制，无需对图像序列作插值运算；同时，该方法重采样间距可以任意设定，在交互过程中可采取低采样策略以提高绘制速度，交互完毕后采取高采样策略以提高绘制图像质量，可同时满足速度和质量的要求。

本发明的技术方案如下：

基于各向异性体数据的错切变形体绘制方法，如图1所示，包括以下步骤：

步骤1：读入图像数据构造三维体数据场。

步骤1-1：读入二维图像序列构造三维体数据场；设所述二维图像序列为尺寸为M×N的L张二维切片图像。

二维图像(又叫二维切片)序列可以是包含图像附加信息的DICOM文件或BMP文件等，也可以是只包含用于绘制的图像数据的二进制文件。为方便后序处理，对于噪声较大的二维图像应进行预处理，如滤波等操作，减少图像噪声，然后再构造三维体数据场。

步骤1-2：建立三维体数据场的世界坐标；

对于M×N×L的体数据场，其中第k层、第i行、第j列的体数据的三维下标为(i，j，k)，1≤i≤M，1≤j≤N，1≤k≤L。设二维切片图像内横向和纵向的像素间距为cx和cy，相邻二维切片图像的间距为cz，则该三维体数据场的世界坐标C_w表示为：

$C_w=\left[\begin{array}{cccc}\mathrm{cx}& 0& 0& 0\\ 0& \mathrm{cy}& 0& 0\\ 0& 0& \mathrm{cz}& 0\\ 0& 0& 0& 1\end{array}\right]\·\left[\begin{array}{c}i\\ j\\ k\\ 1\end{array}\right]$

步骤2：错切变形分解。

步骤2-1：设置初始视矩阵M_view(通常初始视矩阵M_view可设置为单位矩阵)；

步骤2-2：根据鼠标移动，设置视矩阵M_view；

将二维平面鼠标的移动映射为三维球体上向量的转动，根据三维球体上向量的转动角度和方向，利用OpenGL的glRotatef()函数和glGetFloatv()函数得到旋转矩阵，进而得到视矩阵M_view。

步骤2-3：寻找主轴，得到世界坐标到标准世界坐标的转换矩阵P；

计算视矩阵M_view所对应的观察点在世界坐标系中的坐标值，其中绝对值最大的坐标值对应的轴即为主轴。由于错切变形体绘制的重采样是在二维平面中进行，所以需要在三维数据场中确定重采样平面，与重采样平面垂直的轴称为主轴。主轴定义为世界坐标系中与观察视线夹角最小的轴。根据三维空间几何原理，夹角最小等价于夹角的余弦值最大，也就等价于视线上的观察点在该坐标轴上的投影值绝对值最大。所以只需要计算出观察点在世界坐标系中的坐标，绝对值最大的坐标值对应的轴即为主轴。标准世界坐标就是以主轴为第三坐标轴的右手坐标系，通过计算得到的主轴，即可以得到世界坐标到标准世界坐标的转换矩阵P。

步骤2-4：计算标准视矩阵M′_view；

标准视矩阵M′_view的计算公式为M′_view＝M_view·P^-1。视矩阵M_view为世界坐标到观察坐标的转换矩阵，标准视矩阵M′_view为标准世界坐标到观察坐标的转换矩阵，则有关系M′_view＝M_view*P^-1。标准视矩阵M′_view即用于错切变形分解的矩阵。

步骤2-5：采用平行投影方式分解标准视矩阵M′_view；

将标准视矩阵分解为错切矩阵M_shear和变形矩阵M_warp，分解公式为M′_view＝M_warp·M_shear。其中错切矩阵M_shear用于将标准世界坐标空间转换为错切空间，错切空间中观察视线与主轴平行。错切矩阵的形式如下：

$M_{\mathrm{shear}}=\left[\begin{array}{cccc}1& 0& \mathrm{sx}& {\mathrm{ti}}^{\′}\\ 0& 1& \mathrm{sy}& {\mathrm{tj}}^{\′}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$

其中sx，sy为错切因子，用于将视线方向错切到与主轴平行的方向；ti′、tj′为平移因子，用于将错切空间的原点平移至中间图像的最左上角以方便计算。由于错切因子表示的是单位错切距离，而在本发明中的标准世界坐标系中，切片间的距离并不是单位距离，所以平移因子ti′、tj′需要根据切片间的距离做比例运算进行校正。校正公式为：

${\mathrm{ti}}^{\′}=\▿*\mathrm{ti},$ ${\mathrm{tj}}^{\′}=\▿*\mathrm{tj}$

其中，ti、tj分别表示传统错切变形体绘制方法中的平移因子，

表示主轴方向切片(即重采样平面)间的距离。当Z轴为主轴时，

X轴为主轴时，

Y轴为主轴时， $\▿=\mathrm{cy}.$

步骤3：重采样。

采用平行投影方式，在错切空间中，确定视线经过二维切片图像的重采样点；并利用线性加权方式计算每一个重采样点的数据值。数据场各向异性时，二维切片图像内原始数据点成矩形阵列分布，每一个重采样点的采样权重不一样。但是在一行内，垂直方向权重一致；在一列内，水平方向权重一致。

步骤4：合成中间图像。

步骤4-1：不透明度校正；

对沿视线内的每一个重采样点，根据其数据值和不透明度传递函数获取每一个重采样点的不透明度值α，并对每一个重采样点的不透明度值α进行校正，校正公式为：

${\mathrm{\α}}_{\mathrm{corrected}}=1-e^{-\mathrm{\φ\ΔS}}$

$=1-{\left[e^{-\mathrm{\φ\ΔS}0}\right]}^{\frac{\mathrm{\Δs}}{\mathrm{\Δs}0}}$

$=1-{[1-\mathrm{\α}]}^{\frac{1}{\cos \mathrm{\θ}}*\frac{\mathrm{\Δc}}{\mathrm{\ΔS}0}}$

其中，α_corrected为校正后的不透明度值，θ为视线与主轴的夹角，Δc为重采样平面间的距离，Δs₀为不透明度值α所对应的基准距离。

步骤4-2：合成中间图像；

对沿视线内的每一个重采样点的颜色值c和校正后的不透明度值α_corrected，利用前向合成算子合成中间图像的像素值，合成公式为：

$\left\{\begin{array}{c}{\mathrm{\α}}_{\mathrm{out}}={\mathrm{\α}}_{\mathrm{in}}+(1-{\mathrm{\α}}_{\mathrm{in}}){\mathrm{\α}}_i\\ c_{\mathrm{out}}=c_{\mathrm{in}}+(1-{\mathrm{\α}}_{\mathrm{in}})c_i\end{array}\right.$

其中，α_in、c_in分别表示累积的不透明度值和颜色值；α_i、c_i分别表示当前重样点的不透明度值和颜色值；α_out、c_out分别表示经过当前重样点后的不透明度值和颜色值。当完成整个数据场区域的合成之后，就形成了中间图像。

步骤5：对中间图像做变形操作，形成最终图像。

对形成的中间图像，通过变形矩阵M_warp做变形操作，最后生成的就是最终图像。

步骤6：将最终的三维效果图像显示到屏幕。

将最终的图像数据利用OpenGL的纹理贴图功能绘制到屏幕上显示。

通过以上步骤，就能够对各向异性体数据直接进行错切变形体绘制，在屏幕上显示三维立体的最终图像，并可以通过鼠标交互旋转，从不同角度观看图像效果。

需要说明的是：

对于步骤3的重采样过程，重采样平面内的重采样点距离可以任意设定。在实际绘制中，当鼠标交互时，由于只是需要做观察方向的定位，所以采用较低的重采样率(对应较高的重采样距离)，以降低计算量，使图像可以以实时速度做旋转操作；在绘制最终图像时，采用较高的重采样率(对应较低的重采样距离)，以提高图像质量。

本发明的创新之处在于：

1、本发明可以直接基于各向异性体数据进行绘制，避免了各向同性处理的插值过程，不但降低了图像的预处理时间，也提高了绘制图像精度；

2、采用可变重采样率的方式，使绘制速度和绘制精度同时得到满足。

附图说明

图1为错切变形体绘制算法流程图；

图2为错切变形体绘制算法基本原理图；

图3为本发明中重采样示意图；

图4为本发明中不透明度校正示意图。

具体实施方式

本发明采样以上技术方案，利用VC开发环境和OpenGL图形库，实现了整个算法。数据源采用了DICOM格式的MRI图像以及二进制的单纯头骨数据做了测试，均绘制出了效果比较理想的三维图像，在鼠标交互时也能达到实时绘制速度。

Claims (2)

1.基于各向异性体数据的错切变形体绘制方法，包括以下步骤：

步骤1：读入图像数据构造三维体数据场；

步骤1-1：读入二维图像序列构造三维体数据场；设所述二维图像序列为尺寸为M×N的L张二维切片图像；

步骤1-2：建立三维体数据场的世界坐标；

对于M×N×L的体数据场，其中第k层、第i行、第j列的体数据的三维下标为(i，j，k)，1≤i≤M，1≤j≤N，1≤k≤L；设二维切片图像内横向和纵向的像素间距为cx和cy，相邻二维切片图像的间距为cz，则该三维体数据场的世界坐标C_w表示为：

$C_w=\left[\begin{array}{cccc}\mathrm{cx}& 0& 0& 0\\ 0& \mathrm{cy}& 0& 0\\ 0& 0& \mathrm{cz}& 0\\ 0& 0& 0& 1\end{array}\right]\·\left[\begin{array}{c}i\\ j\\ k\\ 1\end{array}\right]$

步骤2：错切变形分解；

步骤2-1：设置初始视矩阵M_view；

步骤2-2：根据鼠标移动，设置视矩阵M_view；

将二维平面鼠标的移动映射为三维球体上向量的转动，根据三维球体上向量的转动角度和方向，利用OpenGL的glRotatef()函数和glGetFloatv()函数得到旋转矩阵，进而得到视矩阵M_view；

步骤2-3：寻找主轴，得到世界坐标到标准世界坐标的转换矩阵P；

计算视矩阵M_view所对应的观察点在世界坐标系中的坐标值，其中绝对值最大的坐标值对应的轴即为主轴；

步骤2-4：计算标准视矩阵M′_view；

标准视矩阵M′_view的计算公式为M′_view＝M_view·P^-1；

步骤2-5：采用平行投影方式分解标准视矩阵M′_view；

将标准视矩阵分解为错切矩阵M_shear和变形矩阵M_warp，分解公式为M′_view＝M_warp·M_shear；其中错切矩阵M_shear用于将标准世界坐标空间转换为错切空间，形式如下：

$M_{\mathrm{shear}}=\left[\begin{array}{cccc}1& 0& \mathrm{sx}& {\mathrm{ti}}^{\′}\\ 0& 1& \mathrm{sy}& {\mathrm{tj}}^{\′}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$

其中sx，sy为错切因子，用于将视线方向错切到与主轴平行的方向；ti′、tj′为平移因子，用于将错切空间的原点平移至中间图像的最左上角以方便计算，且平移因子ti′、tj′需要根据切片间的距离做比例运算进行校正，校正公式为：

${\mathrm{ti}}^{\′}=\▿*\mathrm{ti},$ ${\mathrm{tj}}^{\′}=\▿*\mathrm{tj}$

其中，ti、tj分别表示传统错切变形体绘制方法中的平移因子，

表示主轴方向切片，即重采样平面间的距离；当Z轴为主轴时，X轴为主轴时，Y轴为主轴时，

步骤3：重采样；

采用平行投影方式，在错切空间中，确定视线经过二维切片图像的重采样点；并利用线性加权方式计算每一个重采样点的数据值；数据场各向异性时，二维切片图像内原始数据点成矩形阵列分布，每一个重采样点的采样权重不一样；但是在一行内，垂直方向权重一致；在一列内，水平方向权重一致；

步骤4：合成中间图像；

步骤4-1：不透明度校正；

对沿视线内的每一个重采样点，根据其数据值和不透明度传递函数获取每一个重采样点的不透明度值α，并对每一个重采样点的不透明度值α进行校正，校正公式为：

${\mathrm{\α}}_{\mathrm{corrected}}\text{=1-}{[1-\mathrm{\α}]}^{\frac{1}{\cos \mathrm{\θ}}*\frac{\mathrm{\Δc}}{\mathrm{\Δs}0}}$

其中，α_corrected为校正后的不透明度值，θ为视线与主轴的夹角，Δc为重采样平面间的距离，Δs₀为不透明度值α所对应的基准距离；

步骤4-2：合成中间图像；

对沿视线内的每一个重采样点的颜色值c和校正后的不透明度值α_corrected，利用前向合成算子合成中间图像的像素值，合成公式为：

$\left\{\begin{array}{c}{\mathrm{\α}}_{\mathrm{out}}={\mathrm{\α}}_{\mathrm{in}}+(1-{\mathrm{\α}}_{\mathrm{in}}){\mathrm{\α}}_i\\ c_{\mathrm{out}}=c_{\mathrm{in}}+(1-{\mathrm{\α}}_{\mathrm{in}})c_i\end{array}\right.$

其中，α_in、c_in分别表示累积的不透明度值和颜色值；α_i、c_i分别表示当前重样点的不透明度值和颜色值；α_out、c_out分别表示经过当前重样点后的不透明度值和颜色值；当完成整个数据场区域的合成之后，就形成了中间图像；

步骤5：对中间图像做变形操作，形成最终图像；

对形成的中间图像，通过变形矩阵M_warp做变形操作，最后生成的就是最终图像；

步骤6：将最终的三维效果图像显示到屏幕；

将最终的图像数据利用OpenGL的纹理贴图功能绘制到屏幕上显示；

通过以上步骤，就能够对各向异性体数据直接进行错切变形体绘制，在屏幕上显示三维立体的最终图像，并可以通过鼠标交互旋转，从不同角度观看图像效果。

2.根据权利要求1所述的基于各向异性体数据的错切变形体绘制方法，其特征在于，所述步骤3的重采样过程，重采样平面内的重采样点距离可以任意设定；在实际绘制中，当鼠标交互时采用较低的重采样率；在绘制最终图像时，采用较高的重采样率。

Images