CN1916965A - 基于通用图形显示卡的被测体正投影与反投影方法 - Google Patents

Abstract

本发明属于投影探测与成像技术领域，其特点在于该方法含有以下步骤：把被测物体的重建体数据用三维纹理存储表示；预先计算射线出射点至平面探测器上个探测器之间的方向向量数组，并存储到通用图形显示卡中形成方向向量纹理；利用方向向量纹理与方向向量绕旋转轴旋转来实现任一投影角度的正投影过程的一次完成，无中间结果转存；在反投影过程中，利用三维纹理索引技术使虚拟切片组的放置位置不再受到限制。本发明提高了正投影过程的速度与精度，也增加了反投影过程的灵活性。

Description

基于通用图形显示卡的被测体正投影与反投影方法

技术领域

本发明属于被测体投影探测与成像技术领域。

背景技术

基于投影的重建是指在已知某被测体的一系列投影数据的情况下，通过一定的算法来恢复该数据体内部数据的过程。CT(计算机断层扫描)图像重建就是一种典型的基于投影的重建技术，它利用X射线的透射作用得到被测体的投影数据，然后再根据这些投影数据恢复被测体的内部结构。CT在医学影像、工业检测以及安全检查等领域有着十分广泛的应用。

大规模平板探测器(探测器呈面阵排列)是近年来出现的一种X射线探测装置，采用平板探测器的CT设备具有X射线利用率高，采样速度快，空间分辨率高等特点，因此这种探测器被广泛应用于新一代的CT设备上^[1]。平板探测器的使用使得直接采集二维投影数据成为可能，进而人们可以直接重建三维数据体。可以说平板探测器的出现为基于投影的直接三维重建方法奠定了硬件基础。采用平板探测器的CT设备一般使用锥形X射线束，因此以下将以锥束为例进行说明。如图1所示，被测体被放置在图中虚线正方体区域内，射线出射点与平板探测器围绕旋转轴y轴旋转若干角度，得到被测体的若干个投影角度的二维透视投影，然后利用这些投影数据根据一定的重建算法重建出被测体的内部结构，这就是典型的基于投影的CT图像三维重建过程。

通用图形显示卡(以下简称显卡)是PC机的重要硬件设备，随着图形处理芯片对图形、图像处理能力迅速提高，显卡已经突破传统的显示和娱乐功能，以其强大的流数据处理能力和高度的操作灵活性逐渐成为通用科学计算领域中的重要工具^[2]。三维重建在应用于临床或工程时，在重建速度、精度等方面都有严格要求，而显卡技术正是以追求更为实时、更加真实的3D景物模拟为目标的。三维重建中所使用的数据，规模庞大且操作方式简单，非常符合流数据处理对象的特点，这使得基于显卡实现算法加速成为一种适合三维重建特点的硬件加速方法。

在直接三维重建中，我们要重建的是三维数据体，三维数据体有两种存储形式，一种是三维数组，即T3D[长][宽][高]，另一种是一系列二维数组，如长×T2D[宽][高]，在图形显示卡中它们相应地被存储为一个三维纹理和一组二维纹理，不难看出，这组二维纹理中的每一个都相当于三维纹理的一个切片，如果以垂直于重建立方体表面方向建立直角坐标系，如图1，那么这些切片都垂直于坐标系主轴方向，且切片的各边分别平行于其余两个坐标系主轴，因而以下简称这些切片为重建体正切片。虽然存储为三维纹理会使重建更为简单方便，但由于图形显示技术的限制，在已有的技术中重建数据体被存储为后一种形式。

正投影过程与反投影过程是CT图像重建的两个基本组成部分，因而基于通用图形显示卡的CT图像重建加速将主要围绕这两个过程进行。

在基于投影的三维重建中，所谓的正投影是指模拟X射线投影的方式将三维的重建数据体投影到一个二维平面上形成拟透视投影数据的过程。在已有的技术中，正投影过程是利用重建体正切片逐片实现的，如图2所示，在子图a、b、c、d中实线表示重建体正切片，正投影过程按照a、b、c、d的顺序进行；子图e、f、g、h分别对应于a、b、c、d中的探测器平面，表示各重建体正切片在探测器平面上的投影；子图i表示系统的内存或图形处理器的累积缓存，每得到一片重建体正切片的投影都会在i上累加一次，同一投影角度下的所有重建体正切片的投影都累加过后，就得到该投影角度下的模拟投影数据。

在这种正投影方法中，由于重建体数据被存储为一组二维数组，也就是重建体正切片组，因而这些切片的方向在整个重建过程中都是固定的，即垂直于某一主轴(x或y或z)方向，不能随投影角度的改变而改变。又由于X射线与重建体正切片越接近平行，得到的投影结果越不准确，当完全平行时，几乎得不到投影结果，所以，当绕y轴旋转时，只能使用垂直于x轴或z轴的重建体正切片组存储重建体数据。另外，每一个重建体正切片的计算都需要转存一次投影结果数据，由于图形处理器与内存之间的数据转存需要耗费大量时间，因而转存到系统内存会严重影响重建速度，而转存到图形处理器的累积缓存虽然省时，但累积缓存的数据存储位数限制又会降低拟投影数据的精度，进而影响重建数据精度。

反投影过程是指将某一特定投影角度下的二维数据映射到三维的重建数据体的每个体素上的过程。这里的二维数据可以为任意数据，重建数据体与正投影过程定义相同。

投影纹理映射是比较常用的映射技术。所谓纹理映射简单地说就是通过计算某一点的纹理坐标值来索引该点的纹理取值，而这里的投影纹理映射主要是指通过计算在一点透视投影关系下某一体素在二维投影数据平面上的投影位置来取得该体素的更新值。如图3所示，连接某体素v与射线出射点A并延长交探测器平面于Pro(v)位置，Pro(v)位置的纹理取值即可用于更新体素v的值。在计算机具体实施时，投影纹理映射可分为投影纹理坐标计算和纹理映射两个步骤。

纹理坐标可以通过设置纹理矩阵并用初始纹理坐标与该纹理矩阵相乘得到。在投影纹理映射技术中，纹理矩阵可设置为

           TM＝Tans₁·Rot·Tans₂·Proj·Sc·Tans₃

TM表示纹理矩阵，Tans₁、Rot、Tans₂、Proj、Sc、Tans₃为OpenGL等图形开发语言提供的基本操作矩阵。其中，Tans₁代表平移D_AO(D_AO为旋转轨道半径)个单位距离，Rot代表旋转θ角(θ为投影角度)，Tans₂代表平移-D_AO个单位距离，Proj代表透视关系矩阵(由射线锥角决定)，Sc代表放大0.5倍，Tans₃代表平移0.5个单位距离。初始纹理坐标设置为体素坐标v(x，y，z)，则v所对应的投影纹理坐标T_pro(s_T，t_T，p_T，q_T)为：

          T_pro(s_T，t_T，p_T，q_T)＝v(x，y，z，1)·TM

计算得到的投影纹理坐标即可以用来索引探测器平面上的纹理。投影纹理坐标也是一种纹理坐标，如同用数组下标来索引数组元素一样，图形硬件用纹理坐标来索引纹理取值，不同的是，纹理坐标是连续的，当纹理坐标非整时，图形硬件会自动根据我们选择的插值方式由邻近的纹素值计算出该点的纹理值。如图3所示，当选择线形插值时Pro(v)的纹理值，即体素v的更新值d_v计算公式如下：

            d_v＝abI₁+a(1-b)I₂+(1-a)bI₃+(1-a)(1-b)I₄

其中I₁，I₂，I₃，I₄为四个相邻探测器D₁，D₂，D₃，D₄的取值；Pro(v)表示体素v在探测器平面上的映射位置；a，b定义如图3所示。

当然，无论是二维投影纹理索引还是三维纹理索引，OpenGL等图形开发语言都为我们提供了可以直接调用的函数，其中所需的插值运算将由图形硬件自动完成。

我们可以利用纹理索引函数来进行纹理映射操作。这里用到的纹理映射，简单地说，就是将一幅二维图像直接映射到一个二维平面上，所以反投影过程需要将三维的重建数据体分割成若干个二维切片来分别完成投影，我们以下将这些二维切片称为虚拟切片。

在已有的技术中，虚拟切片组的位置必须与重建体正切片组的位置相同，即虚拟切片必须与某一主轴方向(如图1中的x轴或z轴)垂直且切片的各边须分别与其余两个主轴平行，这样才能使反投影得到的数据可以直接用于更新重建体数据体。如图4所示，反投影过程也是利用虚拟切片逐片进行的，图中探测器平面数据分别向每一个虚拟切片(实线)投影，得到的反投影数据用来更新相应位置的重建体正切片。所有的重建体正切片更新完毕后，反投影过程结束^[4]。

在这种反投影方法中，由于使用一组重建体正切片来存储三维数据，造成在反投影过程中虚拟切片组只能在两个垂直于主轴的方向(图1中的x轴和z轴)放置，否则无法更新重建数据体，从而影响了重建的灵活性，给进一步的重建加速制造了障碍。

发明内容

本发明的目的在于充分利用目前通用图形处理器所提供图形绘制技术，提高正投影过程的速度，增加反投影过程实现的灵活性，进而实现对现有方法在重建速度与重建质量方面的改进，并为进一步的改进打下基础。

本发明的特征在于：所述方法分为正投影与反投影两个阶段，其中：

所述的正投影方法依次含有以下步骤：

步骤(1)，把所述被测体置于一个立方体中，并将此立方体在计算机中存储为三维数组，称为重建数据体，在以该立方体的几何中心为原点，以垂直于该重建立方体表面方向建立x、y、z三维直角坐标系；

步骤(2)，以y轴为旋转轴，以位于所述z轴延长线上一点A(0，0，z_A ⁰)为x射线出射点，以呈面阵排列的x射线平板探测器为接收平面，连接x射线出射点A与所述平板探测器上的一个探测器D(x_D ⁰，y_D ⁰，z_D ⁰)，得到该探测器在投影角度为0度时的方向向量(s_D，t_D，p_D)：

$\left[\begin{array}{c}{s}_D\\ t_D\\ p_D\end{array}\right]=\left[\begin{array}{c}{x}_D^0/M\\ y_D^0/M\\ (z_D^0-z_A^0)/M\end{array}\right]$

其中 $M={({\left(x_D^0\right)}^2+{\left(y_D^0\right)}^2+{(z_D^0-z_A^0)}^2)}^{\frac{1}{2}}$

对于所述平板探测器内的所有探测器，则得到一个与该探测器大小相同的方向向量数组；

步骤(3)，利用所述通用图形显示卡将步骤(2)得到的方向向量数组存储为一个二维纹理，其中每个方向向量的三个分量s_D，t_D，p_D分别存储在纹理的R、G、B三个颜色通道中，形成二维方向向量纹理；

步骤(4)，在所述x射线出射点与平板探测器相对位置保持不变的条件下，围绕y轴旋转θ度，该计算机计算各个探测器的方向向量(s，t，p)：

$\left[\begin{array}{c}s\\ t\\ p\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\θ}& 0& \sin \mathrm{\θ}\\ 0& 1& 0\\ -\sin \mathrm{\θ}& 0& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}{s}_D\\ t_D\\ p_D\end{array}\right]$

此时射线出射点A的坐标为(x_A ^θ，y_A ^θ，z_A ^θ)为

$\left[\begin{array}{c}{x}_A^{\mathrm{\θ}}\\ y_A^{\mathrm{\θ}}\\ z_A^{\mathrm{\θ}}\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\θ}& 0& \sin \mathrm{\θ}\\ 0& 1& 0\\ -\sin \mathrm{\θ}& 0& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}0\\ 0\\ -\mathrm{OE}\end{array}\right]$

其中OE为射线出射点A到旋转中心即坐标原点的距离；

步骤(5)，该计算机在X射线出射点A到探测器平面之间按设定的对被测重建数据体的采样密度划分采样切片，再计算X射线出射点到采样点的距离即采样深度为R₁时的纹理采样坐标(x，y，z)：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=R_1\left[\begin{array}{c}s\\ t\\ p\end{array}\right]+\left[\begin{array}{c}{x}_A^{\mathrm{\θ}}\\ y_A^{\mathrm{\θ}}\\ z_A^{\mathrm{\θ}}\end{array}\right]$

同理得到各个探测器在不同采样深度下的纹理采样坐标；

步骤(6)，将存储内容为该重建数据体的三维数组及步骤(5)计算得到的同一探测器不同采样深度的纹理采样坐标分别代入到OpenGL图形开发语言提供的三维纹理索引函数中，得到同一探测器对应的不同采样深度的纹理值，将这些纹理值相加就得到该探测器的投影值，

同理，该计算机可以计算得到在此投影角度θ下其他探测器的投影值；

所述的反投影方法依次含有以下步骤：

步骤(I)，构筑一个与正投影过程相同的立方体，在以该立方体的几何中心为原点，以垂直于该重建立方体表面方向建立x、y、z三维直角坐标系，把该立方体在计算机中存储为三维数组，称为重建数据体，并将此重建数据体分割为多个二维切片，以下把这些二维切片称为虚拟切片；

步骤(II)，当虚拟切片位置与原重建数据体正切片位置相同时，按以下步骤进行，所述重建数据体正切片是指所述立方体的一个二维切片，且该切片垂直于直角坐标系主轴即旋转轴y轴，而切片的各边又分别平行于其余两个坐标系主轴；

步骤(II.a)，对虚拟切片上的任一体素V(x，y，z)，计算其对应的投影纹理坐标T_pro(s_T，t_T，p_T，q_T)：

             T_pro(s_T，t_T，p_T，q_T)＝V(x，y，z，1)·TM

其中TM为纹理矩阵

             TM＝Tans₁·Rot·Tans₂·Proj·Sc·Tans₃

Tans₁、Rot、Tans₂、Proj、Sc、Tans₃为OpenGL图形开发语言提供的基本操作矩阵，其中Tans₁代表平移D_AO个单位距离，D_AO为旋转轨道半径，

Rot代表旋转θ角，θ为投影角度，

Tans₂代表平移-D_AO个单位距离，

Proj代表透视关系矩阵，由OpenGL图形开发语言根据射线锥角决定，

Sc代表放大0.5倍，

Tans₃代表平移0.5个单位距离；

步骤(II.b)，将二维数组和步骤(II.a)所述的投影纹理坐标代入到OpenGL图形开发语言提供的投影纹理索引函数，即输出所述体素V(x，y，z)的更新值d_v，同理，其他体素的更新值也可得到；

步骤(III)，当虚拟切片位置与原重建数据体正切片位置不一致时，按以下步骤进行：

步骤(III.a)，按照步骤(II.a)和步骤(II.b)的方法计算此时虚拟切片上各个体素的更新值；

步骤(III.b)，将步骤(III.a)得到的所有虚拟切片上的更新值依照虚拟切片的空间位置组合为一个虚拟更新数据体，若在原重建数据体坐标系下的坐标为V(x，y，z)，在所述虚拟更新数据体坐标系下的坐标为V′(x′，y′，z′)，则

              V′(x′，y′，z′)＝TN·V(x，y，z)

其中TN为纹理矩阵，根据原重建数据体坐标系和虚拟更新数据体坐标系的位置关系由OpenGL图形开发语言提供的基本操作矩阵依先后顺序相乘等到，所属先后顺序是指该院重建数据体坐标系通过平移或/和旋转坐标的方式重合到虚拟更新数据体坐标系的先后顺序；

步骤(III.c)，将所述的虚拟更新数据体和步骤(III.b)的得到的虚拟更新数据体坐标系下的坐标V′(x′，y′，z′)代入OpenGL图形开发语言提供的三维纹理索引函数，即可得到步骤(II.b)所述的体素V(x，y，z)的更新值d_v；

步骤(IV)，用步骤(II.b)或步骤(III.b)所述的更新值d_v更新原重建数据体的值。

本发明提高了正投影过程的速度与精度，增加了反投影过程的灵活性。

附图说明

图1CT图像直接三维重建系统结构示意图。

图2已有技术的正投影过程示意图。

图3投影纹理映射示意图。

图4已有技术的反投影过程示意图。

图5本发明的概述示意图。

图6硬件设备示意图。

图7方向向量纹理的计算与存储，a.方向向量的计算，b.方向向量的存储。

图8正投影过程。

图9虚拟切片位置与原重建体正切片位置相同(垂直于旋转轴y轴)的示意图，a.x-z方向剖面图，b.y-z方向剖面图。

图10虚拟切片位置与原重建体正切片位置不一致的示意图，a.x-z方向剖面图，b.组合虚拟更新数据体，c.计算原重建数据体坐标系下的坐标。

图11模型重建结果，a.1次迭代后的重建体中心切片，b.3次迭代后的重建体中心切片，c.模型中心切片。

具体实施方式

如图5所示，本发明提出了一种基于通用图形显示卡的重建数据体正投影与反投影方法，该方法的思路在于：将重建体数据用三维纹理存储表示；预先计算方向向量纹理并存储；利用方向向量纹理来实现某一投影角度的整个正投影过程的一次完成，无中间结果转存；在反投影过程中，利用三维纹理坐标插值技术使虚拟切片组的放置位置不再受到限制。

本发明所涉及内容为基于投影的三维重建方法的基本组成部分，在应用于各种具体重建算法时，视情况采取相应的技术内容，算法的其余部分，可采取与传统方法相同的处理办法。算法所使用的投影数据须为平面二维数据，可由平板探测器或其他二维输入设备获得。

本发明虽均以锥形射线束为例，但本发明并不拘泥于锥形射线束。对于其他射束类型，如平行射线束、扇形射线束等，本发明也适用。

如图6所示，本发明的数据来源为平面探测器(如CT的平板探测器)所采集的二维数据，将采集的数据送入安装有可编程通用图形显示卡的普通微型计算机内进行处理。

本发明的具体实现可利用各种三维图形开发语言(如OpenGL、DirectX等)与图形绘制语言(GLSL、Cg等)。使用的通用图形显示卡须支持三维纹理、可编程等技术。

1.正投影过程

在某一投影角度下(以0度为例)，连接射线出射点与每一个探测器都可以得到一个相应的方向向量，对于所有的探测器来说，就可以得到一个与探测器阵列大小相同的方向向量数组，把这个数组存储为一个二维纹理就是方向向量纹理。下面我们以探测器D的方向向量计算为例说明方向向量纹理的计算与存储过程。如图7(a)所示，设射线出射点A的坐标为(0，0，z_A ⁰)，探测器D的坐标为(x_D ⁰，y_D ⁰，z_D ⁰)，则D探测器在投影角度为0度时所对应的方向向量(s_D，t_D，p_D)为

$\left[\begin{array}{c}{s}_D\\ t_D\\ p_D\end{array}\right]=\left[\begin{array}{c}{x}_D^0/M\\ y_D^0/M\\ (z_D^0-z_A^0)/M\end{array}\right]$

其中 $M={({\left(x_D^0\right)}^2+{\left(y_D^0\right)}^2+{(z_D^0-z_A^0)}^2)}^{\frac{1}{2}}$

然后我们将(s_D，t_D，p_D)存储在方向向量纹理中D探测器对应的位置上，在如图7(b)所示，s_D、t_D、p_D三个分量分别存储在二维纹理的R、G、B三个颜色通道中。同理我们可以计算并存储其他探测器所对应的方向向量。

由于在正投影过程中射线出射点与探测器平面的相对位置关系始终保持不变，因而其他投影角度下各个探测器所对应的方向向量可以根据0度时的方向向量(也就是我们存储的方向向量纹理)经计算得到。

这里我们以投影角度为θ时探测器D的模拟投影值计算为例介绍整个正投影过程。首先我们调入前面计算好的方向向量纹理，然后利用纹理索引函数根据探测器D的位置索引出D的方向向量(s_D，t_D，p_D)，最后我们计算θ投影角度时D探测器的模拟投影值。如图8所示，仍假设旋转轴为y轴，则在投影角度为θ时方向向量(s，t，p)为

$\left[\begin{array}{c}s\\ t\\ p\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\θ}& 0& \sin \mathrm{\θ}\\ 0& 1& 0\\ -\sin \mathrm{\θ}& 0& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}{s}_D\\ t_D\\ p_D\end{array}\right]$

投影角度为θ时射线出射点A的坐标为(x_A ^θ，y_A ^θ，z_A ^θ)为

$\left[\begin{array}{c}{x}_A^{\mathrm{\θ}}\\ y_A^{\mathrm{\θ}}\\ z_A^{\mathrm{\θ}}\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\θ}& 0& \sin \mathrm{\θ}\\ 0& 1& 0\\ -\sin \mathrm{\θ}& 0& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}0\\ 0\\ -\mathrm{OE}\end{array}\right]$

其中OE为射线出射点到旋转中心的距离。

这样，我们就可以根据此角度下射线的方向向量计算出采样深度(即射线出射点到采样点的距离)为R₁时的纹理采样坐标

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=R_1\left[\begin{array}{c}s\\ t\\ p\end{array}\right]+\left[\begin{array}{c}{x}_A^{\mathrm{\θ}}\\ y_A^{\mathrm{\θ}}\\ z_A^{\mathrm{\θ}}\end{array}\right]$

当然，由于图形硬件对纹理坐标范围的限制，(x，y，z)还需要标准化之后才能用于纹理采样，具体的标准化方法可根据所采用的具体图形硬件的要求(如归一化)来制定。

同理我们也可以计算得到采样深度为R₂、R₃、R₄时的纹理采样坐标。将存储内容为重建数据体的三维数组及计算得到的纹理采样坐标分别代入到OpenGL等图形开发语言提供的纹理索引函数中，便可以得到同一探测器对应的不同采样深度的纹理值，将这些纹理值相加就得到探测器D的投影值。同理可以计算得到其他探测器在投影角度为θ时的投影值。正投影过程的GPU计算部分伪码如下：

FOR所有探测器

{

   (s_RV，t_RV，p_RV)＝LoadTexture(RV，m，n)；     //m，n为探测器位置编号

                                                     //RV为方向向量纹理

   (s，t，p)＝(s_D，t_D，p_D)×RotateMatrix(θ)；//θ为投影角度

   标准化处理                                      //根据具体设备采取相应办法

   Proj＝0；                                       //初始化投影值Proj

   FOR i＝D_min to D_max                          //D_min、D_max分别为采样深度

   {                                               //最小值与最大值

     (x，y，z)＝i×(s，t，p)+(0，0，-OE)×RotateMatrix(θ)；

     Proj+＝LoadTexture(V，x，y，z)；             // V为重建体纹理

   }

}

其中LoadTexture()为纹理索引函数，RotateMatrix()为旋转矩阵。

如图8所示，在整个正投影过程中，由于使用了类似球面坐标的计算方法，从而形成了沿主投影方向的均匀球面纹理采样，这种方法相对于以往的沿固定方向(主轴方向)的平面采样方法来说，可以更加真实地模拟透视投影的过程，避免了由于投影方向变化所造成的投影失真。而且整个正投影过程都是在一次绘制中进行的，没有中间结果的转存，节省了时间。

2.反投影过程

由于我们采用了三维纹理来存储重建体数据，所以在反投影过程中虚拟切片的位置不再局限于重建体正切片的位置，因此增大了反投影过程实现的灵活性。以下我们以两种常见的反投影形式为例进行说明。

a)虚拟切片位置与原重建体正切片位置相同

既然虚拟切片的位置可以任意放置，那么它也自然可以放置在原重建体正切片位置，这样放置的优点在于，可以用反投影得到的数据直接来更新原建体数据体。

如图9所示，这里我们采用的是垂直于y轴(旋转轴)的放置方式(现有方法不能垂直于y轴放置)。实线表示虚拟切片，它们的边分别与x、z轴平行，且所在平面垂直于y轴。反投影过程仍然按照虚拟切片逐片进行，以任一虚拟切片上的任一体素v(x，y，z)为例，按照传统的投影纹理映射方法设置纹理矩阵TM，计算v的投影纹理坐标

       T_pro(s_T，t_T，p_T，q_T)＝v(x，y，z，1)·TM

利用投影纹理坐标T_pro(s_T，t_T，p_T，q_T)索引二维数据纹理得到体素v的更新数据d_v，则d_v便可以直接用来更新重建数据体。

b)虚拟切片位置与原重建体正切片位置不一致

这种情况相对第一种情况来说更具有普遍性，这里以比较常用的虚拟切片与探测器平面平行放置为例进行说明。如图10(a)所示，图中实线表示虚拟切片，虚线表示原重建体正切片。这种放置方式最大的优点在于反投影映射面积大，更有利于探测器平面数据在反投影过程中的有效利用。

按照虚拟切片逐片进行反投影的计算方法和第一种情况一样，但这时得到的d_v并不能直接用于更新重建数据体，因为此时的虚拟切片与原重建体正切片方向不一致，因此我们将各虚拟切片得到的更新数据依照虚拟切片的空间位置组合为一个虚拟更新数据体，用一个三维纹理进行存储(如图10(b)所示)，则在虚拟切片与探测器平面平行放置的情况下，虚拟更新数据体与重建数据体绕旋转轴(y轴)成θ′角(反投影角度)。那么我们就可以利用三维纹理坐标插值技术来计算原重建体正切片在虚拟更新数据体中的采样位置。如图10(c)所示，这里仍以v点为例进行说明，假设其在原重建体坐标系下的坐标为(x，y，z)，在虚拟更新数据体坐标系下的坐标为(x′，y′，z′)，则

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\end{array}\right]=\left[\begin{array}{ccc}{\cos \mathrm{\θ}}^{\′}& 0& {\sin \mathrm{\θ}}^{\′}\\ 0& 1& 0\\ {-\sin \mathrm{\θ}}^{\′}& 0& {\cos \mathrm{\θ}}^{\′}\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$

跟正投影过程一样，(x′，y′，z′)也需要标准化之后才能用于虚拟更新数据体纹理索引，可根据具体所采用的图形硬件的要求采取相应的标准化方法(如归一化)。索引得到的虚拟更新数据体纹理值就是v点的更新值d_v′，这个值即可用来更新重建数据体对应位置的体素值。

实验仿真结果：

我们采用三维shepp-logan模型作为仿真试验模型，模型参数设计参见文献[5]，模型大小为128×128×128。数据体重建算法采用SART算法^[3]，其正反投影过程采用本技术中提出的基于通用图形处理器的方法来实现。投影角度从-90°至100°间隔2°采样，共有投影数据96个，投影数据大小为128×128。试验结果见下表：

                          表1试验仿真结果

迭代次数	所需时间	均方误差	相关系数	绝对误差
					1	2.3s	0.02758	0.98192	0.10241
3	6.9s	0.01363	0.99086	0.09685

试验用PC配置为CPU：AMD2800+ RAM：512M×2 GPU：Radon1800XT 512M 256bit

归一化均方误差计算公式为： ${\mathrm{\ϵ}}_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-x_i^{\′})}^2/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i^2$

相关系数计算公式为： ${\mathrm{\ϵ}}_2={\left[\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(x_i-\stackrel{\‾}{x})(x_i^{\′}-{\stackrel{\‾}{x}}^{\′})\right]/\left[\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i^{\′}-{\stackrel{\‾}{x}}^{\′})}^2\right]}^{\frac{1}{2}}$

归一化平均绝对误差计算公式为： ${\mathrm{\ϵ}}_3=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|x_i-x_i^{\′}|/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i$

其中x_i为模型体素值， x为模型体素均值；x_i′为重建体体素值， x′为重建体体素均值。重建结果见图11。

Claims (1)

1.基于通用图形显示卡的被测体正投影与反投影方法，其特征在于：所述方法分为正投影与反投影两个阶段，其中：

所述的正投影方法依次含有以下步骤：

步骤(1)，把所述被测体置于一个立方体中，并将此立方体在计算机中存储为三维数组，称为重建数据体，在以该立方体的几何中心为原点，以垂直于该重建立方体表面方向建立x、y、z三维直角坐标系；

步骤(2)，以y轴为旋转轴，以位于所述z轴延长线上一点A(0，0，z_A ⁰)为x射线出射点，以呈面阵排列的x射线平板探测器为接收平面，连接x射线出射点A与所述平板探测器上的一个探测器D(x_D ⁰，y_D ⁰，z_D ⁰)，得到该探测器在投影角度为0度时的方向向量(s_D，t_D，p_D)：

$\left[\begin{array}{c}{s}_D\\ t_D\\ p_D\end{array}\right]=\left[\begin{array}{c}{x}_D^0/M\\ y_D^0/M\\ (z_D^0-z_A^0)/M\end{array}\right]$

其中 $M={({\left(x_D^0\right)}^2+{\left(y_D^0\right)}^2+{(z_D^0-z_A^0)}^2)}^{\frac{1}{2}}$

对于所述平板探测器内的所有探测器，则得到一个与该探测器大小相同的方向向量数组；

步骤(3)，利用所述通用图形显示卡将步骤(2)得到的方向向量数组存储为一个二维纹理，其中每个方向向量的三个分量s_D，t_D，p_D分别存储在纹理的R、G、B三个颜色通道中，形成二维方向向量纹理；

步骤(4)，在所述x射线出射点与平板探测器相对位置保持不变的条件下，围绕y轴旋转θ度，该计算机计算各个探测器的方向向量(s，t，p)：

$\left[\begin{array}{c}s\\ t\\ p\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\θ}& 0& \sin \mathrm{\θ}\\ 0& 1& 0\\ -\sin \mathrm{\θ}& 0& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}{s}_D\\ t_D\\ p_D\end{array}\right]$

此时射线出射点A的坐标为(x_A ^θ，y_A ^θ，z_A ^θ)为

$\left[\begin{array}{c}{x}_A^{\mathrm{\θ}}\\ y_A^{\mathrm{\θ}}\\ z_A^{\mathrm{\θ}}\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\θ}& 0& \sin \mathrm{\θ}\\ 0& 1& 0\\ -\sin \mathrm{\θ}& 0& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}0\\ 0\\ -\mathrm{OE}\end{array}\right]$

其中OE为射线出射点A到旋转中心即坐标原点的距离；

步骤(5)，该计算机在X射线出射点A到探测器平面之间按设定的对被测重建数据体的采样密度划分采样切片，再计算X射线出射点到采样点的距离即采样深度为R₁时的纹理采样坐标(x，y，z)：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=R_1\left[\begin{array}{c}s\\ t\\ p\end{array}\right]+\left[\begin{array}{c}{x}_A^{\mathrm{\θ}}\\ y_A^{\mathrm{\θ}}\\ z_A^{\mathrm{\θ}}\end{array}\right]$

同理得到各个探测器在不同采样深度下的纹理采样坐标；

步骤(6)，将存储内容为该重建数据体的三维数组及步骤(5)计算得到的同一探测器不同采样深度的纹理采样坐标分别代入到OpenGL图形开发语言提供的三维纹理索引函数中，得到同一探测器对应的不同采样深度的纹理值，将这些纹理值相加就得到该探测器的投影值，

同理，该计算机可以计算得到在此投影角度θ下其他探测器的投影值；

所述的反投影方法依次含有以下步骤：

步骤(I)，构筑一个与正投影过程相同的立方体，在以该立方体的几何中心为原点，以垂直于该重建立方体表面方向建立x、y、z三维直角坐标系，把该立方体在计算机中存储为三维数组，称为重建数据体，并将此重建数据体分割为多个二维切片，以下把这些二维切片称为虚拟切片；

步骤(II)，当虚拟切片位置与原重建数据体正切片位置相同时，按以下步骤进行，所述重建数据体正切片是指所述立方体的一个二维切片，且该切片垂直于直角坐标系主轴即旋转轴y轴，而切片的各边又分别平行于其余两个坐标系主轴；

步骤(II.a)，对虚拟切片上的任一体素V(x，y，z)，计算其对应的投影纹理坐标T_pro(s_T，t_T，p_T，q_T)：

            T_pro(s_T，t_T，p_T，q_T)＝V(x，y，z，1)·TM

其中TM为纹理矩阵

            TM＝Tans₁·Rot·Tans₂·Proj·Sc·Tans₃

Tans₁、Rot、Tans₂、Proj、Sc、Tans₃为OpenGL图形开发语言提供的基本操作矩阵，其中

Tans₁代表平移D_AO个单位距离，D_AO为旋转轨道半径，

Rot代表旋转θ角，θ为投影角度，

Tans₂代表平移-D_AO个单位距离，

Proj代表透视关系矩阵，由OpenGL图形开发语言根据射线锥角决定，

Sc代表放大0.5倍，

Tans₃代表平移0.5个单位距离；

步骤(II.b)，将二维数组和步骤(II.a)所述的投影纹理坐标代入到OpenGL图形开发语言提供的投影纹理索引函数，即输出所述体素V(x，y，z)的更新值d_v，同理，其他体素的更新值也可得到；

步骤(III)，当虚拟切片位置与原重建数据体正切片位置不一致时，按以下步骤进行：

步骤(III.a)，按照步骤(II.a)和步骤(II.b)的方法计算此时虚拟切片上各个体素的更新值；

步骤(III.b)，将步骤(III.a)得到的所有虚拟切片上的更新值依照虚拟切片的空间位置组合为一个虚拟更新数据体，若在原重建数据体坐标系下的坐标为V(x，y，z)，在所述虚拟更新数据体坐标系下的坐标为V′(x′，y′，z′)，则

                   V′(x′，y′，z′)＝TN·V(x，y，z)

其中TN为纹理矩阵，根据原重建数据体坐标系和虚拟更新数据体坐标系的位置关系由OpenGL图形开发语言提供的基本操作矩阵依先后顺序相乘等到，所属先后顺序是指该院重建数据体坐标系通过平移或/和旋转坐标的方式重合到虚拟更新数据体坐标系的先后顺序；步骤(III.c)，将所述的虚拟更新数据体和步骤(III.b)的得到的虚拟更新数据体坐标系下的坐标V′(x′，y′，z′)代入OpenGL图形开发语言提供的三维纹理索引函数，即可得到步骤(II.b)所述的体素V(x，y，z)的更新值d_v；

步骤(IV)，用步骤(II.b)或步骤(III.b)所述的更新值d_v更新原重建数据体的值。

Images