CN102499705B - 一种用于消除层析成像中环状伪影的方法及系统 - Google Patents

Abstract

本发明公开了一种用于消除层析成像中环状伪影的方法，包括：A、对层析成像中360度投影图像分别进行负对数运算；B、将负对数运算后的图像进行叠加；C、根据计算的旋转轴在投影图上的投影位置PRA，计算叠加图像的对称性；D、获得不正常像素APs的位置；E、判断是否仍存在不正常像素，如果有，则执行下一步；如果没有不正常像素，则完成环状伪影的消除；F、判断该不正常像素是否为新出现的不正常像素，如果为新出现的不正常像素，则使用插值或低通滤波方法在所有的投影图像中进行修正；如果为上次迭代中就存在的不正常像素，则按比例从每张投影图像中进行修正；G、执行步骤C。本发明公开了一种用于消除层析成像中环状伪影的系统。

Description

一种用于消除层析成像中环状伪影的方法及系统

技术领域

本发明涉及成像技术领域，尤其涉及一种用于消除层析成像中环状伪影的方法及系统。

背景技术

随着X射线探测器制造工艺日趋成熟，X射线断层成像(computedtomography，CT)技术越来越多地应用于生物医学研究中，如小动物成像、乳房X射线CT、牙科CT，介入治疗检测等。除此之外，CT技术在工业的无损检测方面也有广泛的应用。

然而，X射线探测器由于有大量的像素点(百万量级)，很难保证所有的像素对X射线有相同的响应。由于常用的CT重建算法是滤波反投影算法，因此这种响应的不一致会在重建图像中带来严重的环状伪影(MedicalPhysics，Vol.28，No.5，May2001)。

现有的去除环状伪影的方法主要可以分为两类，分别是在重建图像中进行校准和在投影图中进行校准。

在重建图像中进行校准的方法又可分为两类：一类是使用标准体模，即在成像前使用均匀的体模(比如水或者空气)进行成像，然后以体模的重建结果做为模板直接对样品成像重建后的图像进行校准(美国专利：5774519，中国发明专利申请号200510103923.0，中国发明专利申请号：200980130829.0)，这种方法的缺点在于对于时变的响应不正常点是无效的；另一类是直接对样品成像重建后的图像进行校准，这种方法比较直观，它或者直接对圆环进行匹配滤波(中国发明专利申请号：201010524976.0，中国发明专利申请号：201010524976.0)，或者是将重建图像转化到极坐标系内对直线进行滤波(Physics in Medicine and Biology54(2009)3881-3895，中国发明专利专利号：ZL200410069802.4，中国发明专利专利号：ZL200710031403.2)，这类方法的缺点比较明显，即对强度较高的图像伪影无效，并且会导致图像细节的丢失。

在投影图中进行校准主要是使用正弦图校准，其将投影图在不同角度下的投影结果放到同一图像中，由于响应不均匀的点在不同投影角度下有类似的表现，因此在正弦图中，不正常的点将表现为一条直线，通过对直线的检测或滤波即可将伪影去除(Physics in Medicine and Biology55(2010)6911-6930)。这种方法的缺点在于对强度较弱的不正常像素点检测灵敏度有限，并且有时难于将响应不正常点与图像中原有的信号区分，从而导致错误的判断。

发明内容

有鉴于此，本发明的目的在于提供一种用于消除层析成像中环状伪影的方法及系统，用于高效彻底的消除层析成像中的环状伪影。

本发明提供了一种用于消除层析成像中环状伪影的方法，包括：

A、对层析成像中360度投影图像分别进行负对数运算；

B、将负对数运算后的图像进行叠加；

C、根据计算的旋转轴在投影图上的投影位置PRA，计算叠加图像的对称性；

D、获得不正常像素APs的位置；

E、判断是否仍存在不正常像素，如果有，则执行下一步；如果没有不正常像素，则完成环状伪影的消除；

F、判断该不正常像素是否为新出现的不正常像素，如果为新出现的不正常像素，则使用插值或低通滤波方法在所有的投影图像中进行修正；如果为上次迭代中就存在的不正常像素，则按比例从每张投影图像中进行修正；

G、执行步骤C。

本发明还提供了一种用于消除层析成像中环状伪影的系统，其特征在于，包括：

图像重建装置，用于对层析成像中360度投影图像分别进行负对数运算，并将负对数运算后的图像进行叠加；

全角度投影叠加对称性检测装置，用于根据计算的PRA，计算所述叠加图像的对称性；

响应不正常像素点修正装置，用于获取不正常像素的位置，并对不正常像素进行修正。

本发明是基于层析成像中全角度投影图像的叠加对称性，即通过检测投影叠加图像沿旋转轴投影位置的对称性确定响应不正常的像素点位置，从而进行各个投影图像的校准。本发明能准确找出投影图中响应不正常的像素点位置，在校准时针对响应不正常像素点进行处理，无需对正常像素点进行处理，因此在消除成像伪影的时候能保持成像细节信息，并且不会引入新的误差。

附图说明

图1为本发明实施例提供的一种用于消除层析成像中环状伪影的方法流程图；

图2为本发明实施例提供的获取不正常像素位置方法的原理图；

图3为本发明实施例提供的获取不正常像素点的位置的方法流程图；

图4为本发明实施例提供的获取不正常像素直线的位置的方法流程图；

图5为本发明实施例提供的一种用于消除层析成像中环状伪影的系统结构图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明作进一步的详细描述。

本发明实施例基于层析成像中全角度投影图像的叠加对称性，即通过检测投影叠加图像沿旋转轴投影位置的对称性确定响应不正常的像素点位置，从而进行各个投影图像的校准。本发明实施例能准确找出投影图中响应不正常的像素点位置，在校准时针对响应不正常像素点进行处理，无需对正常像素点进行处理，因此在消除成像伪影的时候能保持成像细节信息，并且不会引入新的误差。

图1为本发明实施例提供的一种用于消除层析成像中环状伪影的方法流程图，包括：

步骤101、对层析成像中360度投影图像分别进行负对数运算。

$p(x,y,\mathrm{\θ})=-\ln \left(\frac{I}{I_0(x,y,\mathrm{\θ})}\right)$

其中p(x，y，θ)是在投影角度为θ时的负对数运算结果，I是入射X射线强度，I₀(x，y，θ)是投影角度为θ的投影图上各点测量的出射X射线强度。

步骤102、将负对数运算后的图像进行叠加。将每个角度的p(x，y，θ)进行求和叠加得到SOP(i，j)：

$\mathrm{SOP}(i,j)=\underset{\mathrm{\θ}=0}{\overset{2\mathrm{\π}}{\mathrm{\Σ}}}p(i,j,\mathrm{\θ})$

其中，SOP(i，j)是投影叠加图像，i，j分别表示投影图或投影叠加图像的行坐标和列坐标。

步骤103、根据计算的旋转轴在投影图上的投影位置(projection ofrotation axis，PRA)，计算叠加图像的对称性。

本实施例定义投影叠加图像SOP的对称性SY(i，j)由下式给出：

SY(i，j)＝SOP(i，j)-SOP(i，n-j+1)i∈[1，m]，

其中，SOP是投影叠加图像，(i，j)表示SOP中像素点的坐标位置，m和n分别表示SOP的行数和列数，假设n是偶数并且旋转轴在投影图中的投影位置位于直线

SY(i，j)为叠加图像的对称性，i，j分别表示投影图或投影叠加图像的行坐标和列坐标。

步骤104、获得不正常像素的位置。图2是获取不正常像素位置方法的原理图，由于不正常像素点在对称性图像中表现出局部的奇异性，因此使用奇异性检测的方法获得不正常像素点的位置；同时由于部分不正常像素点可能位于关于旋转轴对称的直线上，因此通过在叠加图像中检测线奇异性可获得不正常像素直线。

其中，获取不正常像素点的位置的方法如图3所示，包括：

步骤10410、首先使用二维二进小波变换对SY进行分解。

步骤10411、由于不正常像素点(abnormal pixels，APs)表现为SY中局部的奇异性，因此使用一个阈值的方法可以初步分离出APs。在SY中手动选取一块区域，这块区域要求是没有任何成像物体的并且没有任何可见的APs，阈值th_WT即选定为本区域内第一层小波系数的最大值。由于靠近PRA的APs会带来更为严重的伪影，因此本实施例对阈值乘了一个权重函数，使得越靠近PRA的位置使用的阈值越低：

$\mathrm{thr}\_\mathrm{WT}=(1-\frac{1}{x/20+4})\×\mathrm{th}\_\mathrm{WT}$

这里x是像素点与PRA的以像素个数计的距离。

步骤10412、凡是第一层小波系数大于thr_WT的位置被初步标记为APs。位置标记图中如果有空洞则对其进行填充。

步骤10413、由于使用的小波有一阶消失矩，因此有些APs旁边的像素点也会被错误地认为APs。本实施例再采用另一个阈值的方法来排除这些点。首先定义另一个参数pa_SY如下式所描述：

$\mathrm{pa}\_\mathrm{SY}(i,j)=|\mathrm{SY}(i,j)-\frac{1}{8}\mathrm{SN}(i,j)|$

这里SN(i，j)是SY(i，j)的8连通域求和。上述初步确定的APs不会参与8连通域求和中。

步骤10414、与th_WT类似，阈值thr_SY是pa_SY在上述手动选定区域内的最大值。它也乘以相同的权重函数：

$\mathrm{thr}\_\mathrm{SY}=(1-\frac{1}{x/20+4})\×\mathrm{th}\_\mathrm{SY}$

使用这个阈值，将会进一步确定APs的位置。

步骤10415、上面的步骤能够确定APs对的位置，然而无法确认一对AP中哪个位置的像素是异常点，因此需要进一步的处理，首先对标记APs的图像MA进行镜像扩展：

$\mathrm{nMA}(i,j)=\left\{\begin{array}{cc}\mathrm{MA}(i,j)& \mathrm{ifj}\≤\frac{n}{2}\\ \mathrm{MA}(i,n-j+1)& \mathrm{ifj}>\frac{n}{2}\end{array}\right.i\∈[1,m],j\∈[1,n]$

其中(i，j)代表MA中像素点的坐标位置，m和n分别表示SOP的行数和列数。

步骤10416、然后使用拉普拉斯解(LS方法)，以nMA矩阵为模板对SOP进行处理。

步骤10417、再将修正后的SOP和未修正的SOP相减，那些变化的绝对值小于

的将被认为是正常点。

步骤10418、最后对nMA中存在的空洞进行填充。

此外，考虑到在平板中会有线状的APs，而有的线状的APs是关于PRA对称，这时使用上述的方法就会无效。因此本实施例还进一步获取不正常像素直线，其方法如图4所示，包括：

步骤10419、使用二维的二进小波变换对SOP进行分解。

步骤10420、使用基于小波变换和分型的方法检测SOP中每一点的奇异性。李氏指数为负，且小波系数高于噪声水平的点会被初步认为APs，并使用矩阵SAP来标记APs，SAP为非正常响应像素点的位置标记图像。

步骤10421、然后沿着线状APs的方向对SAP进行叠加：

$\mathrm{LA}\left(i\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{SAP}(i,j)$

步骤10422、由于通常线状的APs会比其临近的直线有更多的APs，因此通过曲线LA使用阈值的方法容易地检测出线状APs的位置。而在图像重建中由于使用FDK算法，图像的滤波正是沿着这些线状APs的方向，因此这些线状的APs不会带来显著的成像伪影。因此，本实施例使用了一个高阈值，这个阈值被选为是LA中位数的5倍。

针对异常点的位置图和异常线的位置图进行或运算，从而得到最终APs位置矩阵。

步骤105、判断步骤104的执行结果中是否仍存在不正常像素，如果有，则判断是否为新出现的不正常像素点；如果没有不正常像素，则完成环状伪影的消除。如果为新出现的不正常像素点，则使用插值或低通滤波方法(如LS方法)在所有的投影图像中对其进行修正；如果为上次迭代中就存在的不正常像素点，则按比例从每张投影图中进行修正。

判断是否有新的APs，如果是新的APs，则使用LS方法在每张投影图中对相应位置进行修正；如果是原有的APs则按比例从每张投影图中对相应位置进行修正。计算每张投影图片的修正因子的方法为：

PV_re(i，j)＝PV_c(i，j)+SY/N

这里PV_c(i，j)是投影图经LS方法修正后的结果，而PV_re(i，j)是按比例修正的结果，SY是点(i，j)在对称性结果(SY(i，j))中进行低通滤波或插值的结果。修正后的结果将重新计算投影叠加，并进行新一轮的异常点检测。上述判断将迭代进行直到没有APs存在，这时图像的APs已经全部被补偿完毕。

执行步骤103。

为证明以上方法的有效性，使用基于平板探测器的锥形束微型CT进行成像，并使用本发明的方法进行伪影去除。图2是使用微型CT对一只小鼠成像的结果，从图2(a)和(c)中能看到环状伪影，而在使用本实施例的方法进行处理后如图2(b)和(d)中已经完全消除环状伪影。图3是使用微型CT对鸡大腿骨进行成像的结果。从图3(a)中可以看到环状伪影，而使用本实施例方法处理后从图3(b)中已经将伪影彻底消除。可见本实施例的方法对在体小动物、离体的骨组织都有好的伪影消除效果，证明该方法有较好的适应性。除了环状伪影，由于本方法是直接对不正常像素点进行处理，因此可以同时消除不正常像素点带来的放射线状的伪影。此外本方法也可同样适用于工业CT中，包括kV和MV的工业CT。

本发明实施例还提供了一种用于消除层析成像中环状伪影的系统，如图4所示，包括：

图像重建装置41，用于对层析成像中360度投影图像分别进行负对数运算，并将负对数运算后的图像进行叠加；

全角度投影叠加对称性检测装置42，用于根据计算的PRA，计算所述叠加图像的对称性；

响应不正常像素点修正装置43，用于获取不正常像素的位置，并对不正常像素进行修正。

该系统能准确找出投影图中响应不正常的像素点位置，在校准时针对响应不正常像素点进行处理，无需对正常像素点进行处理，因此在消除成像伪影的时候能保持成像细节信息，并且不会引入新的误差。

总之，以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (10)

1.一种用于消除层析成像中环状伪影的方法，其特征在于，包括：

A、对层析成像中360度投影图像分别进行负对数运算；

B、将负对数运算后的图像进行叠加；

C、根据计算的旋转轴在投影图上的投影位置PRA，计算叠加图像的对称性；

D、获得不正常像素APs的位置；

E、判断是否仍存在不正常像素，如果有，则执行下一步；如果没有不正常像素，则完成环状伪影的消除；

F、判断该不正常像素是否为新出现的不正常像素，如果为新出现的不正常像素，则使用插值或低通滤波方法在所有的投影图像中进行修正；如果为上次迭代中就存在的不正常像素，则按比例从每张投影图像中进行修正；

G、执行步骤C。

2.根据权利要求1所述的用于消除层析成像中环状伪影的方法，其特征在于，所述对层析成像中360度投影图像分别进行负对数运算的方法具体包括：

$p(x,y,\mathrm{\θ})=-\ln \left(\frac{I}{I_0(x,y,\mathrm{\θ})}\right)$

其中p(x，y，θ)是在投影角度为θ时的负对数运算结果，I是入射X射线强度，I₀(x，y，θ)是投影角度为θ的投影图上各点测量的出射X射线强度。

3.根据权利要求2所述的用于消除层析成像中环状伪影的方法，其特征在于，所述将负对数运算后的图像进行叠加的方法具体包括：

将每个角度的p(x，y，θ)进行求和叠加得到SOP(i，j)：

其中，SOP(i，j)是投影叠加图像，i，j分别表示投影图或投影叠加图像的行坐标和列坐标。

4.根据权利要求3所述的用于消除层析成像中环状伪影的方法，其特征在于，所述计算叠加图像的对称性的方法具体包括：

SY(i，j)＝SOP(i，j)-SOP(i，n-j+1)i∈[1，m]，

其中，SOP是投影叠加图像，(i，j)表示SOP中像素点的坐标位置，m和n分别表示SOP的行数和列数，假设n是偶数并且PRA位于直线

SY(i，j)为叠加图像的对称性，i，j分别表示投影图或投影叠加图像的行坐标和列坐标。

5.根据权利要求4所述的用于消除层析成像中环状伪影的方法，其特征在于，所述获得不正常像素的位置的方法包括使用奇异性检测的方法获得不正常像素点的位置以及通过在叠加图像中检测线奇异性的方法获得不正常像素直线的位置。

6.根据权利要求5所述的用于消除层析成像中环状伪影的方法，其特征在于，所述获得不正常像素点的位置的方法具体包括：

使用二维二进小波变换对SY进行分解；

在SY中手动选取一块区域，该区域没有任何成像物体且没有任何可见的APs，阈值th_WT选定为该区域内第一层小波系数的最大值；

对阈值乘一个权重函数，使得越靠近PRA的位置使用的阈值越低：

其中，x是像素点与PRA的以像素个数计的距离；

将第一层小波系数大于thr_WT的位置初步标记为APs，位置标记图中如果有空洞则对其进行填充；

定义参数pa_SY：

$\mathrm{pa}\_\mathrm{SY}(i,j)=|\mathrm{SY}(i,j)-\frac{1}{8}\mathrm{SN}(i,j)|,$ SN(i，j)是SY(i，j)的8连通域求和；

阈值thr_SY是pa_SY在所述手动选定区域内的最大值，并乘以相同的权重函数：

$\mathrm{thr}\_\mathrm{SY}=(1-\frac{1}{x/20+4})\×\mathrm{th}\_\mathrm{SY}$

对标记APs的图像MA进行镜像扩展：

$\mathrm{nMA}(i,j)=\left\{\begin{array}{cc}\mathrm{MA}(i,j)& \mathrm{ifj}\≤\frac{n}{2}\\ \mathrm{MA}(i,n-j+1)& \mathrm{ifj}>\frac{n}{2}\end{array}\right.i\∈[1,m],j\∈[1,n],$ 其中(i，j)表示MA中像素点的坐标位置，m和n分别表示SOP的行数和列数；

使用拉普拉斯解LS方法，以nMA矩阵为模板对SOP进行处理；

将修正后的SOP和未修正的SOP相减，变化的绝对值小于

的像素点为是正常点；

对nMA中存在的空洞进行填充。

7.根据权利要求5或6所述的用于消除层析成像中环状伪影的方法，其特征在于，所述获得不正常像素直线的位置的方法具体包括：

使用二维的二进小波变换对SOP进行分解；

使用基于小波变换和分型的方法检测SOP中每一点的奇异性，初步设定李氏指数为负且小波系数高于噪声水平的点为APs，并使用矩阵SAP来标记APs，SAP为非正常响应像素点的位置标记图像；

沿着线状APs的方向对SAP进行叠加：

$\mathrm{LA}\left(i\right)=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{SAP}(i,j)$

通过曲线LA使用阈值的方法检测出线状APs的位置，该阈值为是LA中位数的5倍。

8.根据权利要求7所述的用于消除层析成像中环状伪影的方法，其特征在于，获取了不正常像素点和不正常像素直线后还进一步包括：

针对不正常像素点的位置和不正常像素直线的位置进行或运算，得到最终APs位置矩阵。

9.根据权利要求8所述的用于消除层析成像中环状伪影的方法，其特征在于，所述在投影图像中进行修正的方法中，获取每张投影图片的修正因子的方法为：

PV_re(i，j)＝PV_c(i，j)+SY/N

其中，PV_c(i，j)是投影图经LS方法修正后的结果，而PV_re(i，j)是按比例修正的结果，SY是点(i，j)在对称性结果(SY(i，j))中进行低通滤波或插值的结果。

10.一种用于消除层析成像中环状伪影的系统，其特征在于，包括：

图像重建装置，用于对层析成像中360度投影图像分别进行负对数运算，并将负对数运算后的图像进行叠加；

全角度投影叠加对称性检测装置，用于根据计算的PRA，计算所述叠加图像的对称性；

响应不正常像素点修正装置，用于获取不正常像素的位置，并对不正常像素进行修正。

Images