CN101958006B - 一种基于x光影像的物体三维成像方法 - Google Patents

Abstract

本发明公开了一种基于X光影像的物体三维成像方法，包括以下步骤：步骤一，X光数据采集，由X光影像采集设备环绕着待重建物体每隔一定的固定角度间隔拍摄一幅X光平面二维影像从而构成二维影像序列{P_t}，其中t＝1，2，3...n，每幅X光平面二维影像中包括X光影像采集设备的探测器中的像素下标(i，j)，即像素对应的列位置和行位置；步骤二，三维物体重建，将采集的X光平面二维影像进行重建，获得待重建物体的三维体素表示；步骤三，三维物体显示，对物体的三维体素表示进行渲染从而获得待重建物体的三维影像。本发明利用若干幅X光影像重建出相应物体的三维表示的技术，从而提供了一种更加全面、生动、立体的观察物体内部结构的方法。

Description

一种基于X光影像的物体三维成像方法

技术领域

本发明属于计算机应用技术领域，特别是涉及一种基于X光影像的物体三维成像方法。

背景技术

由于三维模型可以提供比二维图像更多的空间信息和观察维度，因此，三维模型近年来越来越多地应用到工业、医疗、教育等领域。X光影像是获取物体内部结构信息的一种重要的手段，但其影像所呈现的结果为物体的二维投影，因此，无法从更多的观察角度去了解物体的情况，而且X光片所形成的影像中物体内部在射线方向的投影重叠在一起，无法直观地了解物体内部各处的情况。CT是另一种得到广泛应用的医疗诊断器材，它可以用于了解物体内部的情况，但是具有辐射量大、无法应用手术过程中的诊断等缺点。

现有的X光影像采集设备包括发射源和探测器，如图1所示，只能完成单个平面的X光影像采集，如何利用现有的采集设备，在成本最低，对现有设备改动最小的情况下，实现三维成像是目前领域研究的一大难点。

发明内容

发明目的：本发明所要解决的技术问题是针对现有技术的不足，提供一种基于X光影像的物体三维成像方法，从而利用一组的二维平面X光影像，通过本发明处理，获得相应物体在计算机内的三维表示，从而实现在计算机上对物体进行三维观察。

为了解决上述技术问题，本发明公开了一种基于X光影像的物体三维成像方法，包括以下步骤，分别是：

步骤一，X光数据采集，X光影像采集设备(包括发射源和探测器)环绕着待重建物体每隔一定的固定角度间隔拍摄一幅X光平面二维影像P_t，并存储为DICOM文件格式，从而构成二维影像序列{P_t}，其中t＝1，2，3...n，每幅X光平面二维影像中包括X光影像采集设备的探测器中的像素下标(i，j)，即像素对应的列行位置；

步骤二，三维物体重建，将采集的X光平面二维影像进行重建，获得待重建物体的三维体素表示；

步骤三，三维物体的显示，对物体的三维体素表示进行渲染从而获得待重建物体的三维影像。

本发明步骤一中，X光数据采集包括：

(1)X光影像采集设备环绕对物体进行环绕拍摄。

(2)在数据采集的过程中，使用X光影像采集设备对物体进行环绕拍摄，每间隔θ度(通常为5度或2.5度)拍摄一幅X光平面二维影像，则可以得到一组连续的二维影像序列序列{P_t}，t＝1，2，3...n。

(3)环绕拍摄过程中X光影像采集设备具有如下的参数：

a)X光射线发射源到探测器的距离为D(单位为厘米)；

b)X光射线发射源到被拍摄物体的中心的距离为d(单位为厘米)；

c)探测器的分辨率为M×N个像素，每个像素所对应的正方形探测单元边长大小为Δ(单位厘米)；

由于三维重建过程是在世界坐标系下进行的，因此，需要将探测器上的下标转换为直角坐标系中的坐标，换言之，将探测器中的像素下标(i，j)转换为直角坐标系中的(x，y)。以像素的中心位置作为该像素的坐标位置，依据上述的参数设定，可以得到下面的计算公式可得(x，y)：

$\left\{\begin{array}{cc}x=(i-\frac{N}{2})\·\mathrm{\Δ},& i=\mathrm{0,1},...,N-1\\ y=(\frac{M}{2}-j)\·\mathrm{\Δ},& j=\mathrm{0,1},...,M-1\end{array}\right.---\left(1\right)$

其中，Δ为每个像素所对应的正方形探测单元边长，M×N为探测器的分辨率。

步骤二，三维物体重建由如下步骤组成：

(1)投影数据的权重修正：对于X光平面二维影像P_t上的每个像素值(i，j)，其修正权重的计算公式为：

$w(i,j)=\frac{D^2}{\sqrt{D^2+x^2+y^2}}---\left(2\right)$

公式(2)中的(x，y)由公式(1)得到。利用公式(2)得到的权重对每幅X光二维影像P_t上的每个像素的值进行修正，即：

${\hat{P}}_t(i,j)=P_t(i,j)\·w(i,j)=\frac{P_t(i,j)D^2}{\sqrt{D^2+{(\mathrm{i\Δ}-\mathrm{N\Δ}/2)}^2+{(\mathrm{M\Δ}/2-\mathrm{j\Δ})}^2}}---\left(3\right)$

由于权重矩阵对于每一幅X光平面二维影像都是相同的，因此，可以预先计算出w(i，j)存储在一个二维数组中，从而加速运算过程。

(2)投影数据的滤波。采用一维滤波器(例如ram-lak、shepp-logan、cosine等滤波器)对每幅X光平面二维影像的每一行进行滤波操作。假设滤波器为g，则对二维影像序列{P_t}(为方便公式化描述可采用

代替)，每一行数据的滤波操作实质上是指该行数据与g的卷积操作，，即如下的公式所示：

${{\hat{P}}^{\′}}_t(i,\·)={\hat{P}}_t(i,\·)*g---\left(4\right)$

其中(i，·)表示第i行所有数据，符号*表示卷积操作。为了加速上述的计算过程，通常可以采用快速傅里叶变换，将空间域里的卷积操作变换到频域后用乘法运算代替，即

${{\hat{P}}^{\′}}_t(i,\·)=F^{-1}\left\{F\right[{\hat{P}}_t(i,\·)]\·F[g\left]\right\}---\left(5\right)$

公式(5)中的F(·)和F^-1(·)分别表示傅里叶变换和反变换。

(3)计算重建物体的大小(即包含待重建物体的长方体区域大小)。为了便于计算机处理，重建结果通常使用三维数组表示，数组的每一个元素表示三维物体在空间中一个体素的值。换言之，重建结果应该为一个长方体(底面为正方形)。为了保证重建物体都置于X光射束之内，需要对长方体的三条边的大小进行计算。重建物体xy平面上的边长a的计算公式为：

$\frac{0.5a}{0.5\mathrm{N\Δ}}=\frac{d-0.5a}{D}\⇒a=\frac{\mathrm{dN\Δ}}{D+0.5\mathrm{N\Δ}}---\left(6\right)$

重建物体在z轴方向上的大小c为：

$c=\frac{\mathrm{dM\Δ}}{D+0.5\mathrm{M\Δ}}---\left(7\right)$

(4)体素模型的数组表示，将重建结果存放在数组V[XDIM][YDIM][ZDIM]中(其中XDIM，YDIM，ZDIM为设置参数，一般情况下只需保证探测器分辨率M×N中M、N的值分别都比XDIM(或YDIM)、ZDIM大，在此基础上可固定一组参数，比如当探测器分辨率为500×400时可给出一个参考取值为300、300、500，具体取值可在前述原则上通过实验进行适当调整)，则有：

1)重建物体的几何中心位于三维数组元素V[XDIM/2][YDIM/2][ZDIM/2]，与X光影像采集设备旋转中心重合；

2)每个体素在x，y，z方向上的大小分别为：

${\mathrm{\Δ}}_x^r=\frac{a}{\mathrm{XDIM}},$ ${\mathrm{\Δ}}_y^r=\frac{a}{\mathrm{YDIM}},$ ${\mathrm{\Δ}}_z^r=\frac{c}{\mathrm{ZDIM}}---\left(8\right)$

3)以体素的中心位置作为该体素的坐标位置，三维数组V中的元素V[i][j][k]所对应的体素空间坐标为：

$\left\{\begin{array}{c}x=(\mathrm{XDIM}/2-i)\·{\mathrm{\Δ}}_x^r\\ y=(j-\mathrm{YDIM}/2)\·{\mathrm{\Δ}}_y^r\\ z=(k-\mathrm{ZDIM}/2)\·{\mathrm{\Δ}}_z^r\end{array}\right.---\left(9\right)$

其中V[0][0][0]为包含待重建物体的长方体区域最下层的左下角处的体素。

(5)计算反投影坐标，该过程所涉及的问题是，如何将第t幅X光平面二维影像的数据反投影到体素模型中去。依据体素位置，计算其在二维影像P_t上投影的位置，将该位置所对应的像素值加到体素中去。体素的空间坐标为(x，y，z)，当前要进行反投影的X光平面二维影像P_t的拍摄角度为φ，则该体素在X光平面二维影像Pt上对应的坐标为：

$(x_p,y_p)=(\frac{D(y\cos \mathrm{\φ}-x\sin \mathrm{\φ}}{d+x\cos \mathrm{\φ}+y\sin \mathrm{\φ}},\frac{\mathrm{Dz}}{d+x\cos \mathrm{\φ}+y\sin \mathrm{\φ}})---\left(10\right)$

由公式(10)得到的空间坐标，可以进一步采用插值的方式得到该坐标所对应的像素投影值，本发明采用最近邻的方法插值，即：

其中i、j分别表示X光平面二维影像P_t上行坐标、列坐标位置令V[i][j][k]＝V[i][j][k]+P_t(i，j)，将每幅X光平面的二维影像反投影到三维体素模型中，获得待重建物体的三维体素表示。

(6)重建结果的去噪声，为了提高重建结果的鲁棒性，使用中值滤波对重建结果进行去噪声，窗口设置为3，即对体素V[i][j][k]，找到它的3-邻居集合，对集合中的值进行排序(升降序都可)，然后将V[i][j][k]设置为排序结果中排在中间位置的那个值。

步骤三，三维物体显示，显示的设备为通常的液晶显示屏，可以采用常用的光线跟踪算法来对重建结果进行体绘制。所述光线跟踪算法可以参见文献：[1]Appel A.Oncalculating the illusion of reality.Proceedings of IFIP Congress 1968，Volume 2，Edinburgh，1968：945-950或者文献：[2]Donald Hearn，M.Pauline Baker著.蔡士杰，宋继强，蔡敏等译.计算机图形学(第三版).北京：电子工业出版社.第482～495页。

有益效果：实验表明，采用本发明可以重建出物体的三维视图，这是X光平面二维影像所无法观察得到的。而且，由于多了一个可控制的观察维度，因此，观察者可以从更多的角度去了解被测物体的情况。本发明利用若干幅X光影像重建出相应物体的三维表示的技术，从而提供了一种更加全面、生动、立体的观察物体内部结构的方法，可以应用到医疗辅助诊断和工业产品检测等领域中。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1本发明一个实施例中X光平面二维影像采集设备结构示意图。

图2本发明实施例中二维影像采集设备探测器的像素下标与空间坐标变换示意图。

图3本发明计算重建物体大小的示意图。

图4本发明的方法流程图。

图5本发明反投影过程的流程图。

具体实施方式：

图1为本发明所述X光影像采集设备的大体构造，包括旋转臂4上的X光发射源1和探测器2，该设备的X光发射源发射锥型X光束，穿透被拍摄物体3后到达设备的探测器，由探测器检测到残余的X射线强度，并借助计算机处理最后得到拍摄物体的X光平面二维影像。X光影像采集设备的旋转臂4可以顺着滑槽沿图中双向箭头的方向移动，从而形成对被拍摄物体的环绕拍摄。通常总是拍摄的角度间隔是相同的。在本发明中，X光发射源与被拍摄物体的中心的距离用d表示，X光发射源与探测器的距离用D表示。

图2为本发明探测器中像素的下标与其空间坐标变换的示意图。探测器下标的坐标系为：左上角为原点，水平方向向右为正方向，垂直方向向下为正方向。而在空间坐标中，坐标原点为探测器的中心位置，垂直方向上的正方向为向上。若探测器的分辨率为M×N个像素，每个像素所对应的正方形探测单元边长大小为Δ(单位厘米)，则可以通过公式(1)来将下标转换为空间坐标，以像素的中心位置作为该像素的坐标位置，依据上述的参数设定，可以得到下面的计算公式可得(x，y)：

$\left\{\begin{array}{cc}x=(i-\frac{N}{2})\·\mathrm{\Δ},& i=\mathrm{0,1},...,N-1\\ y=(\frac{M}{2}-j)\·\mathrm{\Δ},& j=\mathrm{0,1},...,M-1\end{array}\right.---\left(1\right)$

其中，Δ为每个像素所对应的正方形探测单元边长，M、N与探测器的分辨率，M×N中值分别对应。

图3a和图3b为本发明计算重建物体在xy切面上大小的示意图。本发明只对置于X光辐射范围内的部分进行三维重建，这是因为当重建物体的大小超出X光辐射范围时会有信息缺失或冗余，进而影响重建精度，因此，需将重建物体的大小限制在X光射束内部。图3a中假设重建结果为一个底面为正方形的长方体。

如图3b所示，点O为底面为正方形的长方体的底面中心，重建物体xy平面上的边长a的计算公式为，E点位X光发射源位置：

$\frac{0.5a}{0.5\mathrm{N\Δ}}=\frac{d-0.5a}{D}\⇒a=\frac{\mathrm{dN\Δ}}{D+0.5\mathrm{N\Δ}},$

重建物体在z轴方向上的大小c为： $c=\frac{\mathrm{dM\Δ}}{D+0.5\mathrm{M\Δ}}.$

图4为本发明流程图，包括六大步骤：步骤10采集X光平面二维影像序列，将待检测的部位暴露在X光发射源和探测器之间，每间隔一定的角度拍摄一幅X光平面二维影像；步骤11使用公式(3)对X光平面二维影像的上每个像素值进行权重修正；步骤12采用一维滤波器对每幅X光平面二维影像的每一行进行滤波；步骤13将经过权重修正和滤波后的X光平面二维影像反投影到三维体素模型中；步骤14使用中值滤波算法对重建后的三维模型进行去除噪声；最后步骤15将最终重建三维物体重建结果显示出来，可以采用常用的光线跟踪算法来对重建结果进行体绘制。

因现在的图像开发库(如VTK)可以提供光线跟踪算法的方法类，在具体应用时可直接在重建体数据基础上根据具体开发方法调用。对屏幕上每一像素，执行下述3步操作：

1)：从视点出发通过该像素中心向场景发出一条光线，并求出该条光线与场景中物体的全部交点，

2)：将各交点沿光线方向排序，获得离视点最近交点，

3)：依据局部光照明模型计算该交点处的光亮度，并将所得光亮度值赋给该像素。

当所有屏幕像素都处理完毕时，即得到一幅真实感图形。

图5a、图5b为本发明中反投影算法的过程：

步骤20使用公式(6)和(7)计算重建结果在x轴的长度和在y轴上的长度均为a，在z轴方向上的长度为c，即

$a=\frac{\mathrm{dN\Δ}}{D+0.5\mathrm{N\Δ}}---\left(6\right)$

用M、c分别替代(6)式中的N、a，可以得到重建物体在z轴方向上的大小c为

$c=\frac{\mathrm{dM\Δ}}{D+0.5\mathrm{M\Δ}}---\left(7\right)$

步骤21使用公式(8)算出重建结果的体素大小

即

${\mathrm{\Δ}}_x^r=\frac{a}{\mathrm{XDIM}},$ ${\mathrm{\Δ}}_y^r=\frac{a}{\mathrm{YDIM}},$ ${\mathrm{\Δ}}_z^r=\frac{c}{\mathrm{ZDIM}}---\left(8\right)$

步骤22初始化三维数组V[XDIM][YDIM][ZDIM]为0，

步骤23设置控制X光平面二维影像循环的计数器t，初值赋为0，

步骤24检测计数器t是否小于X光平面二维影像数，若小于则进入步骤25，否则结束反投影过程，

步骤25初始化重建结果第一维下标计数器i为0，

步骤26检测计数器i是否小于XDIM，若小于则进入步骤28，否则转步骤27，

步骤27令计数器t加1，

步骤28初始化重建结果第二维下标计数器j为0，

步骤29检测计数器j是否小于YDIM，若小于则进入步骤31，否则转步骤30，

步骤30令计数器i加1，

步骤31初始化重建结果第三维下标计数器k为0，

步骤32检测计数器k是否小于ZDIM，若小于则进入步骤34，否则转步骤33，

步骤33令计数器j加1，

本发明提供了一种基于X光影像的物体三维成像方法的思路及方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确说明的各组成部分均可用现有技术加以实现。

Claims (2)

1.一种基于X光影像的物体三维成像方法，其特征在于，包括以下步骤：

步骤一，X光数据采集，由X光影像采集设备环绕着待重建物体每隔一定的固定角度间隔拍摄一幅X光平面二维影像从而构成二维影像序列{P_t}，其中t＝1，2，3...n，每幅X光平面二维影像中包括X光影像采集设备的探测器中的像素下标(i，j)，即像素对应的列位置和行位置；

步骤二，三维物体重建，将采集的X光平面二维影像进行重建，获得待重建物体的三维体素表示；

步骤三，三维物体显示，根据待重建物体的三维体素表示获得待重建物体的三维影像；

步骤一中，将探测器置于直角坐标系中，从而将探测器中的像素下标(i，j)转换为直角坐标系中的坐标(x，y)；

以像素的中心位置作为该像素的坐标位置，探测器中的像素下标(i，j)的直角坐标(x，y)为：

$\left\{\begin{array}{c}x=(i-\frac{N}{2})\·\mathrm{\Δ},i=\mathrm{0,1},...,N-1\\ y=(\frac{M}{2}-j)\·\mathrm{\Δ},j=\mathrm{0,1},...,M-1\end{array}\right.;$

其中，Δ为每个像素所对应的正方形探测单元边长，M×N为探测器的分辨率；

步骤二包括如下步骤：

对X光平面二维影像上的每个像素值进行权重修正；

采用一维滤波器对每幅X光平面二维影像的每一行进行滤波；

将上步骤中得到的二维影像反投影到三维体素模型中；

去除噪声，对三维体素V[i][j][k]，找到其3-邻居集合，对3-邻居集合中的值进行排序，然后将三维体素V[i][j][k]设置为排序结果中排在中间位置的值；

将二维影像反投影到三维体素模型中包括如下步骤：

计算待重建物体的大小，得到待重建物体在x轴的长度和在y轴上的长度均为a，在z轴方向上的长度为c，其中：

$a=\frac{\mathrm{dN\Δ}}{D+0.5\mathrm{N\Δ}},$ $c=\frac{\mathrm{dM\Δ}}{D+0.5\mathrm{M\Δ}},$

D为X光影像采集设备的X光射线发射源到探测器的距离，d为X光影像采集设备的X光射线发射源到待重建物体的几何中心的距离；

体素模型的数组表示，初始化三维数组V为0，并将重建结果存放在三维数组V中，待重建物体的中心对应于三维数组元素V[XDIM/2][YDIM/2][ZDIM/2]，其中XDIM、YDIM，ZDIM分别对应三维体素模型中用来度量x，y，z方向长度的体素个数，每个体素在x，y，z方向上的大小分别为

其中：

${\mathrm{\Δ}}_x^r=\frac{a}{\mathrm{XDIM}},$ ${\mathrm{\Δ}}_y^r=\frac{a}{\mathrm{YDIM}},$ ${\mathrm{\Δ}}_z^r=\frac{c}{\mathrm{ZDIM}},$

以体素的中心位置作为该体素的坐标位置，三维数组V中的元素V[i][j][k]所对应的体素空间坐标为：

$\left\{\begin{array}{c}x=(\mathrm{XDIM}/2-i)\·{\mathrm{\Δ}}_x^r\\ y=(j-\mathrm{YDIM}/2)\·{\mathrm{\Δ}}_y^r\\ z=(k-\mathrm{ZDIM}/2)\·{\mathrm{\Δ}}_z^r\end{array}\right.,$

其中V[0][0][0]为包含待重建物体的长方体区域最下层的左下角处的体素；

计算反投影坐标，体素的空间坐标为(x，y，z)，当前要进行反投影的X光平面二维影像的拍摄角度为φ，则该体素在X光平面二维影像上对应的坐标为：

$(x_p,y_p)=(\frac{D(y\cos \mathrm{\φ}-x\sin \mathrm{\φ})}{d+x\cos \mathrm{\φ}+y\sin \mathrm{\φ}},\frac{\mathrm{Dz}}{d+x\cos \mathrm{\φ}+y\sin \mathrm{\φ}}),$

采用插值法得到该坐标所对应的像素值，即：

其中i、j分别表示X光平面二维影像上的行坐标和列坐标；

令V[i][j][k]＝V[i][j][k]+P_t(i，j)，将每幅X光平面的二维影像反投影到三维体素模型中，获得待重建物体的三维体素表示。

2.据权利要求1所述的一种基于X光影像的物体三维成像方法，其特征在于，对X光平面二维影像上的每个像素值进行权重修正采用如下公式：

${\hat{P}}_t(i,j)=P_t(i,j)\·\mathrm{\ω}(i,j)=\frac{P_t(i,j)D^2}{\sqrt{D^2+{(\mathrm{i\Δ}-\mathrm{N\Δ}/2)}^2+{(\mathrm{M\Δ}/2-\mathrm{j\Δ})}^2}},$

其中，表示加权后X光平面二维影像上第i行第j列的像素值，ω(i，j)表示对应的权重。

Images