用于 CT 图像断层重建的滤波函数建立方法及断层重建方 法 技术领域 本发明涉及图像处理技术领域， 尤其涉及一种用于 CT(ComputedTomography， 计算 机断层扫描 ) 图像断层重建过程的滤波函数建立方法及利用该方法进行 CT 图像断层重建 的方法。
背景技术 FBP(Filter Back Projection， 滤波反投影 ) 是一种用于 CT 图像重建过程的非常 经典的方法， 由于其具备计算量小、 图像重建速度快等特点， 已经广泛应用于工业 CT 图像 重建的各个领域。作为 FBP 方法实施过程中非常重要的阶段， 滤波过程是保证 CT 图像重建 质量的关键， 因此选择一个好的滤波函数就显得非常重要。目前应用比较广泛的滤波函数 有 1973 年 Ramachandran.G.N， Lakshminarayanan.A.V 所提出的 R-L 滤波函数和 1974 年 Shepp L.A 和 Logan B.F 所提出的 S-L 滤波器。
由于工业 CT 图像重建过程中存在各种各样的要求， 比如有的要求局部成清晰图 像， 有的要求消除 Gibbs( 吉布斯 ) 效应 ( 将具有不连续点的周期函数 ( 如矩形脉冲 ) 进行 傅立叶级数展开后， 选取有限项进行合成。 当选取的项数越多， 在所合成的波形中出现的峰 起越靠近原信号的不连续点。 当选取的项数很大时， 该峰起值趋于一个常数， 大约等于总跳 变值的 9％。这种现象称为吉布斯效应 )， 有的要求提高图像的空间分辨率以及有的要求 提高图像的密度分辨率等等。针对这些不同的实际要求， 需要设计出相适应的滤波函数来 满足不同的特定情况。目前， 研究人员已经在新型滤波函数的研发上做了很大的努力， 存 在有针对射线源强度不均匀性、 射线束硬化以及散射等影响而设计的滤波函数， 也存在有 基于 SR-CT(Synchrotron Radiation ComputedTomography， 同步辐射计算机断层扫描 ) 反 投影算法设计的滤波函数， 还存在有专门针对 CT 图像局部重建的滤波函数， 更存在有针对 Gibbs 现象设计的滤波函数。
尽管如前所述， 目前已存在有众多具备较强针对性的新型滤波函数， 但是它们基 本上均没有给出一个面对 CT 图像重建过程整体情况的设计思路和设计方向， 而其在通常 情况下， 一般仅仅只是给出一个函数， 抑或是给出一个基于该滤波函数的频域分布图和空 域分布图。因此目前所存在的滤波函数， 基本上均为直接借用已有的 R-L 滤波函数和 S-L 滤波函数来改造， 其不能满足某些特殊重建效果的需要 ； 此外， 目前所存在的特殊滤波函数 均为针对不同实际需求而设计， 其一般是根据个人经验进行的尝试， 不具备普适性。
发明内容 ( 一 ) 要解决的技术问题
本发明要解决的技术问题是 ： 提高 CT 图像重建过程中所应用的滤波函数的普适 性， 提高图像的重建质量。
( 二 ) 技术方案
为解决上述技术问题， 本发明提供了一种用于 CT 图像断层重建过程的滤波函数 建立方法，
S1 ： 确定作为滤波函数建立基础的基函数 ：
hbase(x) ＝ Ramp(x)×wbase(x) ；
其中， Ramp(x) 为斜坡函数， wbase(x) 为空域窗函数 ； 斜坡函数用于消除 CT 图像断 层重建过程中因反投影过程引入的星形伪影， 窗函数决定要建立的滤波函数的振幅大小和 截止频率 ；
S2 ： 根据精度需要， 选择基函数的个数， 并对所确定的基函数在空域 ( 时域 ) 进行 平移， 左移或右移若干个采样频率后， 选取一组基函数构建成基函数组 ： hbase(x-kd) ；
其中， k ＝ 0， ±1， ±2， ......±n， k 为采样频率个数， d 为空间采样频率， n 是基 函数的个数， hbase(x-kd) 表示索引号为 k 的基函数 ；
S3 ： 对所述基函数组进行加权求和， 得到在空域上建立的滤波函数 ：
其中， wk 表示索引号为 k 的基函数在所要建立的滤波函数中所占有的比重。 S4 ： 将所述滤波函数从空域向频域转换， 得到在频域上的滤波函数 ：其中，为空间采样间隔， i 是复数单位， ρ 是空间频率， Wbase(ρ) 为频域窗函数。 其中， 比重 wk 的选取依据为 ： 若要提高空间分辨率， 则提高索引号为偶数的基函数 的比重 wk ； 若要提高密度分辨率， 则提高索引号为奇数的基函数的比重 wk。
其中， 所述基函数 hbase 也为滤波函数。
其中， 所述基函数 hbase 为 R-L 滤波函数或 S-L 滤波函数。
其中， 所述斜坡函数为绝对值函数。
其中， 所述窗函数为 R-L 滤波函数。
本发明还提供了一种 CT 图像断层重建方法， 包括以下步骤 ：
S1’ 、 利用斜坡函数消除 CT 图像断层重建过程中因反投影过程引入的星形伪影 ；
S2’ 、 利用所述斜坡函数以及所选取的窗函数按照上述的方法构建滤波函数 ；
S3’ 、 利用所构建的滤波函数进行 CT 图像断层重建。
在步骤 S3’ 中， 通过滤波反投影算法进行 CT 图像断层重建。所述步骤 S3’ 具体 为： 利用所构建的滤波函数对 CT 图像进行滤波处理 ； 然后通过逆投影变换公式进行对处理 后的图像进行反投影变换。
( 三 ) 有益效果
本发明技术方案对比现有技术， 利用基滤波函数作为构建新滤波函数模型的基函 数， 采用平移若干个空间采样频率的方法来构建基函数族， 并基于平均加权的思想来设计 滤波函数， 通过本发明技术方案所设计出来的滤波函数， 在频域上反应为一组傅里叶余弦 级数， 证明了它的普适性。
附图说明 图 1 为本发明实施例的滤波函数建立方法的流程图 ；
图 2 为本发明具体实施方式中所涉及的图像采集装置 ；
图 3 为本发明具体实施方式中所涉及的加权平均方式示意图 ；
图 4 为本发明实施例中所涉及的滤波函数在频域上分布示意图 ；
图 5 为本发明实施例中所涉及的滤波函数在频域上分布示意图 ；
图 6 为本发明实施例中所涉及的滤波函数在空域远旁瓣上的分布示意图 ；
图 7 为本发明实施例中所涉及的滤波函数与已有的滤波函数在重建效果上的比 较示意图 ；
图 8 为本发明实施例中所涉及的滤波函数与已有的滤波函数在重建图像的中间 行灰度值分布情况示意图 ；
图 9 为本发明实施例中所涉及的滤波函数与已有的滤波函数在重建图像的海螺 断层的重建效果示意图 ；
图 10 为本发明实施例中所涉及的滤波函数与已有的滤波函数在重建图像的中间 行像素元的比较情况示意图。
具体实施方式 为使本发明的目的、 内容、 和优点更加清楚， 下面结合附图和实施例， 对本发明的 具体实施方式作进一步详细描述。
如图 1 所示， 本发明技术方案所提供的用于 CT 图像断层重建过程的滤波函数建立 方法， 包括如下步骤 ：
S1 ： 确定作为滤波函数建立基础的基函数 ：
滤波函数一般可以分为平滑型， 复原型和增强型， 平滑型滤波函数一般用于增强 的空间滤波， 其目的是消除噪声 ； 复原型滤波函数有 Wienner(WN) 滤波函数和 Metz(MZ) 滤 波函数， 这两种滤波函数有增强信号的功能， 也称为增强型滤波函数。由 FBP 重建算法的推 导过程可知， 在 FBP 算法重建过程中， 一般都采用平滑型滤波函数， 其主要原因是因为平滑 型滤波函数可以很好的消除因 FBP 重建过程中引入的星型伪影。FBP 算法重建过程中需要 有一个斜坡函数 (Ramp 函数 ) 来消除因反投影引起的星形伪影， 这也是引入滤波函数的最 重要的目的。因此一个滤波函数中， Ramp 函数部分是必须要求存在的。一般而言， 一个平 滑型滤波函数包括 Ramp 函数和窗函数。 Ramp 函数作用是消除反投影过程中的星形伪影， 而 窗函数则决定振幅大小和截止频率。该滤波模型设计需要选择一基函数作为设计基础 ：
hbase(x) ＝ Ramp(x)×wbase(x) ；
其中， Ramp(x) 为斜坡函数， wbase(x) 为窗函数 ；
作为基函数， hbase(x) 自己本身也应该是一个滤波函数， 可以选用常见的 R-L 滤波 函数， S-L 滤波函数或者其它常用的滤波函数等。
S2 ： 如图 3 所示， 对确定选用的基函数在空域进行平移， 根据精度需要， 选择基函 数的个数， 左移和右移若干次后， 选取其中一部分平移操作所对应得到的平移后的基函数， 构建为一组基函数组 ：
hbase(x-kd) ；
其中， k ＝ 0， ±1， ±2， ......±n， k 根据实际需要进行选择， d 为空间采样频率， n 为基函数的个数 ； 这样构成 2n+1 个函数， 作为基函数组。
S3 ： 对所述加权基函数组进行求和， 得到在空域上建立的滤波函数 ：
其中， wk 表示索引号为 k 的基函数在所要建立的滤波函数中所占有的比重。
从信号与系统的角度看， 判断一个滤波函数的好坏主要有一下三个因素 ： 主瓣， 近 旁瓣， 远旁瓣。 从空域的特性图上看， 主瓣越高越窄， 空间分辨率就越好。 远旁瓣的幅度和幅 值越小时， 吉布斯效应就越小， 有利于提高图像的密度分辨率。 因此， 根据实际需要， 通过加 大所述空域上建立的滤波函数偶数位后的函数元的权重来提高 CT 图像空间分辨率 ； 通过 加大所述空域上建立的滤波函数奇数位后的函数元的权重来提高 CT 图像的密度分辨率。
S4 ： 从理论上分析， 将滤波函数从空域向频域进行平移转换， 反映在频域上就是乘 i2πkρ/(2N-1) 以一个波动因子 e 。因此当上述滤波函数在空间邻域 2n+1 点进行加权平均时， 反 应到频域上所建立的滤波函数如下 ：
其中，为空间采样间隔， i 为复数单位， ρ 表示空间频率， 是空间长度的倒 ＝ 1 是归一化的表达式。这样构成了一个新的滤波函数设计模型。针对不数。此外，同的实际情况而言， 可以通过修改 wk 来得到满足需要的滤波函数。
当 wk ＝ w-k 时， 上式则可以化为 ：
令 t ＝ 2 πρ /(2N-1) ， 则 cos(2 π k ρ /(2N-1)) ＝ cos(kt) 。 cos(t) ， cos(2t)， ......， cos(kt)， 是一组傅里叶余弦级数， 任何偶函数均可以用这组余弦级数来 表示。在 CT 图像重建的过程中， 窗函数应该都是偶函数， 因此本发明技术方案针对 CT 图像 重建过程所设计提供的滤波函数， 具有适用于各种情况的普适性的。
上述方法所用到的投影数据可通过如图 2 所示的装置进行采集以及实时处理， 来 得到重建图像效果。该装置包括 ： 量程为 1200mm 的电控平移台 1， 量程为 600mm 的电控平 移台 2， 量程为 1200mm 的电控平移台 3， 平板探测器源 4， 连接处 5， 电控旋转台 6， 样品工作 平台 7， 系统平台 8， X 射线源 9。
实施例 1 下面结合附图来说明本发明技术方案在 CT 图像重建过程中具体的一种实施情 (1) 首先选用基滤波函数， 然后将该基函数， 通过平移基滤波函数构建一基函数况。
组。 这里选择 R-L 滤波函数作为窗函数， 因为用于 CT 图像重建的滤波函数的设计 主要是设计窗函数， 当然也可以使用现有的窗函数， 如 Hanning 窗函数， Hamming 窗函数， Blackman 窗函数， Butterworth 窗函数， Gaussian 窗函数。频域中的 R-L 滤波函数 ( 也就 是下文的窗函数 WR-L(ρ)) 表示为 ：
HR-L(ρ) ＝ |ρ|×Rect(ρ)，
其中 ：是矩形窗函数， 其中 N 是指空间采样间隔 ； 该矩形窗函数在空域中表示为 ：
hR-L(x) ＝ 2N2sinc(2xN)-N2sinc2(xN)，
hR-L(x-kd) ＝ {2N2sinc(2N(x-kd))-N2sinc2(N(x-kd))， 其中， k ＝ 0， ±1， ±2， ......±n， K 根据实际需要选取， 这样就构成 2n+1 个函数， 作为基函数组 ； 其中的函数 sinc(X) 的定义为 ：
(2) 构建新滤波函数模型
当 w0 ＝ 0.5， w1 ＝ 0.25 时，是 Hanning 窗函数；
当 w0 ＝ 0.54， w1 ＝ 0.23 时，是 Hamming窗函数 ；
当w0＝0 . 4 2 ，w1＝0 . 2 5 ，w2＝0 . 0 4时 ，是 Blackman 窗函数。
(3) 根据实际需要， 选择合适的 w0， w1， w2 等。下面选用一个包含 5 个基函数的基函数组来构建新函数模型。针对不同情况在频 域和空域中的分布规律选择我们合适的滤波函数， 图 4 是新建滤波函数模型不同权重取值 时在频域中的分布情况， 图 5 是滤波函数在空域主瓣的分布， 图 6 是滤波函数在空域远旁瓣 的分布。根据模型， 使 w0 ＝ 0.7， w1 ＝ 0.10， w2 ＝ 0.05。设计一个新的滤波函数 ：
H(ρ) ＝ |ρ|[0.7-0.05cos(2πρ/(2N-1))+0.025cos(4πρ/(2N-1))]，
运用上述步骤设计的滤波函数， 与已有的 S-L 滤波函数和 R-L 滤波函数对模拟的 Shepp-Logan 的模拟的投影数据进行滤波反投影重建。 图 7 是三种滤波器重建效果的比较， 其中左图为 R-L 滤波函数重建的效果， 中图为本发明技术方案所提供的滤波函数重建的效 果， 右图为 S-L 滤波函数重建的效果 ； 图 8 是重建后的图像的中间行灰度图的分布情况， 其中左图为 R-L 滤波函数重建的情况， 中图为本发明技术方案所提供的滤波函数重建的情 况， 右图为 S-L 滤波函数重建的情况 ； 然后对由图 1 系统得到的关于海螺投影的数据进行重 建， 三种滤波函数重建效果的比较如图 9 所示， 其中左图为 R-L 滤波函数重建的情况， 中图 为本发明技术方案所提供的滤波函数重建的情况， 右图为 S-L 滤波函数重建的情况 ； 图 10 是 S-L 滤波函数进行重建的结果和新设计的滤波器的重建效果在中间行灰度值的细节上 的比较， 其中， 左图为本发明技术方案所提供的滤波函数重建的情况， 右图为 S-L 滤波函数 重建的情况。 综合来看， 新设计的滤波函数的总体性能比已有的两个滤波函数的性能要好。
本发明技术方案从信号与系统的角度出发， 分析滤波函数在频域分辨的特点和其 滤波效果的关系 ： 从信号与系统的角度来看， 判断一个滤波函数的好坏主要有一下三个因 素： 主瓣， 近旁瓣， 远旁瓣。从空域的特性图上看， 主瓣越高越窄， 空间分辨率就越好。远旁 瓣的幅度和幅值越小时， 吉布斯效应就越小， 有利于提高图像的密度分辨率。 借用平均加权 的思想， 建立了相关的函数设计模型， 并发现这一模型具有普适性。 为滤波函数的设计提供 了一种设计思路。
本发明还提供了一种 CT 图像断层重建方法， 包括以下步骤 ：
S1’ 、 利用斜坡函数消除 CT 图像断层重建过程中因反投影过程引入的星形伪影 ；
S2’ 、 利用所述斜坡函数以及所选取的窗函数按照上述函数建立方法构建滤波函 数；
S3’ 、 利用所构建的滤波函数进行 CT 图像断层重建。
在步骤 S3’ 中， 通过滤波反投影算法 (FBP 算法 ) 进行 CT 图像断层重建。所述步 骤 S3’ 具体为 ： 利用所构建的滤波函数对 CT 图像进行滤波处理 ； 然后通过逆投影变换公式 对处理后的图像进行反投影变换， 通过逆投影变换公式对处理后的图像进行反投影变换是 现有技术。
图 4 表示几种滤波函数在频率上的比较， 线表示的 R-L 滤波函数在频域上的分布 情况， 空心圆分布表示设计 w0 ＝ 0.5， w1 ＝ 0.25， w2 ＝ 0 时的滤波函数在频域上的分布， 叉 号分布表示设计 w0 ＝ 0.1， w1 ＝ 0.45， w2 ＝ 0 时的滤波函数在频域上的分布情况， 加号分布 表示设计 w0 ＝ 0.15， w1 ＝ 0.25， w2 ＝ 0.175 时的滤波函数在频域上的分布， 三角形分布表 示设计 w0 ＝ 0.42， w1 ＝ 0.25， w2 ＝ 0.04 时的滤波函数在频域上的分布，
图 5 与图 6 是图 4 中几种滤波函数在空域中的分布情况， 空心圆分布表示设计 w0 ＝ 0.5， w1 ＝ 0.25， w2 ＝ 0 时的滤波函数在频域上的分布， 叉号分布表示设计 w0 ＝ 0.1， w1 ＝ 0.45， w2 ＝ 0 时的滤波函数在频域上的分布情况， 加号分布表示设计 w0 ＝ 0.15， w1 ＝ 0.25，w2 ＝ 0.175 时的滤波函数在频域上的分布， 三角形分布表示设计 w0 ＝ 0.42， w1 ＝ 0.25， w2 ＝ 0.04 时的滤波函数在频域上的分布， 其中图 5 是滤波中心点附近的分布情况， 即主瓣的 分布情况， 图 6 是远离中心点的分布情况， 即远旁瓣的分布情况。
图 7 是针对模拟的 Shepp-Logan 重建的一个效果比较， 其中 a 图是采用 R-L 滤波 函数重建的效果， 参照图 8 中的 (a) 部分 ( 表示中间行的像素灰度值分布 )， 其灰度值变化 幅度比较大， 在空间分辨率中具有严重的噪声， 其空间分辨率很不好， 同时其密度分辨率上 面出现了严重的过曝现象， 其分辨率也很不好。图 7 中 (c) 部分是采用 S-L 滤波函数的成 像效果比较光滑， 其空间分辨率比较好， 但是其不同密度区域的对比度比较差一点。图 7 中 (b) 是采用我们这种设计思想， 在空间分辨率与密度分辨上面的一个折中的情况下设计的 一个新滤波函数的重建效果图。它同时兼顾了空间分辨率与密度分辨率。
图 9 和图 10 是针对实际系统 ( 图 2) 扫描得到的投影图进行重建的效果比较。
以上所述仅是本发明的优选实施方式， 应当指出， 对于本技术领域的普通技术人 员来说， 在不脱离本发明技术原理的前提下， 还可以做出若干改进和变形， 这些改进和变形 也应视为本发明的保护范围。