基于小波变换的多尺度锥束CT图像快速三维重建方法

摘要

一种基于小波变换的多尺度锥束CT图像三维重建快速算法，包括以下步骤：(1)对试件进行圆周锥束CT扫描，采集一组用于重建的投影图像；(2)根据对重建数据分辨率要求的不同，对采集到的投影图像进行相应尺度的小波变换，得到各尺度小波分解系数；(3)选择相应尺度的小波系数进行FDK重建；如果需要对小波一次分解后的所有分量分别进行FDK重建，则继续；(4)将得到的重建数据分别沿着y轴方向取断层图像，进行相应的小波逆变换；(5)将小波逆变换后得到的三维数据沿着z轴方向取断层图像，进而进行线性插值，得到最终重建图像。本发明能够在一定的精度内快速实现基于FDK算法的三维重建，在此过程中能够针对不同的要求来呈现不同分辨率的三维数据。

说明书

技术领域

本发明涉及一种用于计算机图像处理的基于小波变换的多尺度锥束CT图像的三维重建方法，属于计算机图像处理技术领域。

背景技术

X光计算机断层成像技术(X-ray Computer Tomography,XCT)，是利用X射线对物体进行不同角度的扫描而获得物体断层信息的成像技术，该技术涵盖了物理学、数学、计算机视觉、计算机图形学、数字信号处理等多方面的学科知识，并且它被广泛应用在很多科学领域，这使得人们观察物体内部结构的能力极大地增强。尤其是在医学临床诊断中，它被用来当做一种获取人体内部结构信息的最佳手段。

目前主流的商用CT机是扇束CT扫描机，也就是说这种CT机的X射线源每扫描一周只能获取被检测对象某一个断层的数据，这样导致X射线的使用效率低、扫描速度慢及重建速度慢。与扇束CT机相比，锥束CT由于采用了锥束X射线源和面阵探测器来采集被检测对象不同角度的投影数据，大大提高了X射线利用率、扫描速度和CT重建图像的分辨率，在医学临床检测和工业无损检测上的应用越来越广泛。一般地，锥束XCT重建算法分为两类，一类是迭代法，另一类是解析法(如FDK算法)。迭代法的计算量很大、效率低，难以适应现代商业实时性的要求。与迭代法相比，解析算法的重建效率高，需要的数据存储空间小。一般地解析法分为近似重建法和精确重建法，其中FDK算法是从扇束滤波反投影算法推导出来的锥束XCT近似重建算法，在X射线源锥角较小的情况下重建出来的断层图像质量非常好。由于其简单而且易实现的特点，使FDK算法的研究具有实际应用价值。

锥束CT完成一次三维扫描的数据量达几百GB甚至TB，重建算法复杂，运算量大，耗时较多，难以满足实际应用中实时性的要求，特别是在医学应用领域。其次，针对不同的应用和显示载体的限制，对投影数据重建后的三维数据的分辨率的要求是不同的，对一些重建后分辨率要求不高的情况，不需要耗时太多来重建高分辨率图像。

发明内容

针对锥束CT三维重建算法存在的耗时多以及不同的应用对重建后的三维数据的分辨率要求不同等问题，本发明提出在保证精度前提下快速实现三维重建的基于小波变换的多尺度锥束CT图像三维重建快速算法。

本发明的基于小波变换的多尺度锥束CT图像三维重建快速算法，包括以下步骤：

(1)对试件进行圆周锥束CT扫描，采集一组用于重建的投影图像；

(2)根据对重建数据分辨率的要求，对采集到的投影图像进行相应尺度的小波变换，得到各尺度小波分解系数；

(3)选择相应尺度的小波系数进行FDK重建；

已知得到分辨率为2N×2N的投影图像，对投影图像进行小波一次分解，进而对分解后的低频分量进行FDK重建，最终得到分辨率为(N/2)×(N/2)×(N/2)的重建数据；或者对小波一次分解后的所有分量分别进行FDK重建，得到四组重建数据，继续执行步骤(4)和(5)；

以一次小波分解为例，由步骤(1)得到一组分辨率为2N×2N的二维投影图像，对投影图像进行小波一次分解，得到四组分辨率为N×N的小波分解系数，若对重建数据分辨率要求较低，可仅对低频分量进行FDK重建，得到分辨率为(N/2)×(N/2)×(N/2)的三维重建数据；若对重建数据分辨率要求较高，可对小波分解后的所有分量分别进行FDK重建，得到四组重建数据，继续执行步骤(4)和(5)。

(4)将得到的重建数据分别沿着y轴方向取断层图像，进行对应的小波逆变换；

(5)将小波逆变换后得到的三维数据沿着z轴方向取断层图像，进行线性插值，得到最终重建图像。

本发明的方法能够在一定的精度内快速实现基于FDK算法的三维重建，在此过程中能够针对不同的要求来呈现不同分辨率的三维数据。

附图说明

图1是一幅图像在二维频域被分解为四个子带的框图。

具体实施方式

以下详细说明本发明基于小波变换的多尺度锥束CT图像三维重建快速算法的具体过程。

(1)对试件进行圆周锥束CT扫描，采集一组用于重建的投影图像；

锥束扫描及投影数据的参数如下表所示：

参数名称参数射线源到试件中心距离(mm)DSO射线源到探测器距离(mm)DSD探测器长度(mm)2L探测器分辨率2N×2N重建物体分辨率N×N×N重建物体边长长度(mm)L采集投影间隔(°)M采集投影个数360/m

按照上述参数进行锥束扫描，会得到360/m张，分辨率为2N×2N的二维投影图像P_n(u,v)，n＝1,2,…,360/m，u,v分别标记投影图像像素的行标和列标，它们的零点都以平面探测器的中心像素为基准。

(2)根据对重建数据分辨率的要求，对采集到的投影图像进行相应尺度的二维离散小波变换，得到各尺度小波分解系数；

二维离散小波变换：

二维小波变换由一维小波变换扩展而来，由一维尺度函数和相应的小波函数ψ，可得到一个二维尺度函数和三个二维小波ψ^H(x,y)，ψ^V(x,y),ψ^D(x,y)：

ψ^D(x,y)＝ψ(x)ψ(y)，

例如，Haar尺度函数定义为:

小波函数定义为:

那么，尺寸为M×N的函数f(x,y)的离散小波变换是：

$>W_{\psi }^i(j,m,n)=\frac{1}{\sqrt{MN}}\underset{x=0}{\overset{M-1}{\Sigma }}\underset{y=0}{\overset{N-1}{\Sigma }}f(x,y){\psi }_{j,m,n}^i(x,y),i=\{H,V,D\}$>

其中，j₀是任意的开始尺度，系数定义了在尺度j₀的f(x,y)的近似，系数对于j≥j₀附加了水平、垂直和对角方向的细节，通常令j₀＝0并且选择N＝M＝2^J,j＝0,1,2,…，J-1,m,n＝0,1,2,…，2^j-1，。

二维离散小波逆变换：

通过上式给出的和f(x,y)可通过离散反小波变换得到，即

假设初始尺度为j₀的一幅图像f(x,y),在二维频域可被分解为四个子带分别对应于图像f(x,y)的低频成分、垂直方向上的高频成分、水平方向上的高频成分和对角方向上的高频成分。从图1可以看出，分解后的子图像形成一个塔形结构，每层有3个子图，随着阶数的增加，图像的空间分辨率依次在行和列上减半。对低频子带还可以进一步分解，进行一次变换可得到4个子图像，进行k次变换就可得到3k+1个子图像，这就是基于小波变换的多尺度分解。

以一次小波分解为例，对第n幅投影图像P_n(u,v)进行一次小波分解，可分解为四个小波系数，由于投影图像分辨率为2N×2N，经过一次小波变换会分解为四个分辨率为N×N的图像，分别记为P_n1(u,v)，P_n2(u,v)，P_n3(u,v)，P_n4(u,v)对应于低频成分、垂直方向上的高频成分、水平方向上的高频成分和对角方向上的高频成分。所有投影图像同理，因此，得到了四组小波系数。

(3)选择相应尺度的小波系数进行FDK重建；

由于投影数据是离散的，为了便于计算，必须把FDK算法离散化，在这里引入几个标记：P_β(u,v)和分别表示在旋转角β下的投影数据和滤波处理后的投影数据，β是中心射线与Y轴的夹角,β和投影张数n及投影间隔m的关系是：β＝n×m；h(u)表示对第u行投影图像数据进行卷积操作的滤波函数，理论上要求滤波器为斜坡滤波器。具体FDK重建算法的离散化形式，分为以下三个步骤：

第一步：对投影数据P_β(u,v)使用加权因子进行修正，即：其中加权因子是锥束射线入射角的余弦函数，DSO是射线源到试件中心的距离，u,v分别标记投影图像像素的行标和列标，它们的零点都以平面探测器的中心像素为基准。

第二步：对修正后的投影数据逐行地进行一维滤波，即其中*是表示卷积,常用的滤波函数有R-L滤波函数和S-L滤波函数。例如，S-L滤波函数的采样序列为：n＝0,±1,±2,...。

第三步：对滤波后的投影数据做锥形束加权反投影，即：

式中，