CN103177425A - 切伦科夫荧光成像中γ射线的去除方法 - Google Patents

Abstract

本发明公开一种切伦科夫荧光成像中γ射线的去除方法。主要解决切伦科夫成像中γ射线干扰图像质量的问题，其实现步骤是：(1)获取图像数据；(2)图像预处理；(3)根据如下公式计算图像掩膜的宽度；(4)获取局部方差图像；(5)确定被干扰像素；(6)构造信号判定矩阵；(7)修复被干扰像素，获取初步修复图像；(8)判断初步修复图像是否满足处理要求；(9)对初步修复图像抑噪。本发明能够基本完全去除切伦科夫荧光图像中的γ射线，同时图像中切伦科夫荧光信号能够尽可能的保留下来，提高切伦科夫荧光成像的图像质量。

Description

切伦科夫荧光成像中γ射线的去除方法

技术领域

本发明涉及图像处理技术领域，更进一步涉及医学影像处理领域中的切伦科夫荧光成像去除干扰的方法。本发明可用于去除切伦科夫荧光成像过程中存在的γ射线干扰，提高切伦科夫荧光成像的图像质量。

背景技术

在切伦科夫荧光成像中，希望成像图像中荧光信号能够准确对应放射性核素活度，以反映核素标记的靶细胞的增殖和凋亡情况。但是由于产生切伦科夫光的同时，伴随着正电子湮灭或者γ衰变，面临着信号被γ射线干扰的问题，希望在去除γ射线干扰的同时，尽可能保留信号。

由于切伦科夫荧光信号的产生机制的特殊性，直接对原始图像选取阈值以去除干扰是很困难的。常用方法对γ射线干扰的像素使用平均化进行处理，不仅γ射线干扰不能完全去除，而且对周围的信号像素影响也较大。这些方法处理低放射性活度下采集的切伦科夫荧光图像的结果对原始信号影响之大以至于无法进行后续分析。

中国科学院自动化研究所在其申请的专利“基于切伦科夫效应的断层成像方法和系统”(专利申请号201080059913.0，公开号CN102741707A)中提出了一种切伦科夫三维断层成像系统。该专利申请采用高阶简化球谐近似模型和迭代半阈值算子相耦合的方法进行全域CLT重建，该系统仅使用切伦科夫荧光探测装置采集切伦科夫光学成像。该专利申请不足之处是，在采集平面光学图像时并没有考虑如何处理γ射线对切伦科夫成像的干扰。

Gregory S.Mitchell等人在“In vivo Cerenkov luminescence imaging：a new tool formolecular imaging.PhilTrans R Soc A 2011；369：4605-4619.”中提出了在切伦科夫图像采集的过程中通过CCD相机中的宇宙射线去除模块对切伦科夫荧光成像中γ射线进行处理。该方法对采集图像的相机要求较高，在一定程度上能够去除切伦科夫荧光成像过程中的γ射线。这种方法的不足之处有两个，其一是对于不同相机处理的结果不同，处理结果依赖于相机处理模块的性能。其二是不能完全去除γ射线的干扰，尤其是低放射性活度(等效表面10uCi以下)条件下图像处理的质量反而会下降。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种切伦科夫荧光成像中γ射线的去除方法，用于去除切伦科夫荧光成像中γ射线的干扰。

实现本发明的具体思路是，采用迭代处理的方法，逐步增加图像掩膜的宽度，计算图像的局部方差，根据最大类间方差法得到的阈值确定被干扰的像素，利用被干扰像素在掩膜范围内未被干扰像素的平均灰度值替换γ射线的像素值，根据处理的结果判断迭代是否结束，最后利用非截断小波抑噪算法滤除图像的噪声。

为实现上述目的，本发明的具体步骤如下：

(1)获取图像数据

1a)在避光条件下用CCD相机捕获切伦科夫荧光图像；

1b)在避光条件下用CCD相机捕获背景噪声图像；

(2)图像预处理

将采集的荧光图像减去背景噪声图像，得到被干扰图像；

(3)根据如下公式计算图像掩膜的宽度：

w＝2t+1

其中，w表示图像掩膜的宽度大小，t表示迭代的次数；

(4)获取局部方差图像

4a)使用平滑线性滤波器，对被干扰图像进行平滑线性滤波，获得图像中每个像素的局部均值；

4b)利用图像局部方差公式，计算被干扰图像中每个像素的局部方差，获取局部方差图像；

(5)确定被干扰像素

5a)采用最大类间方差法处理局部方差图像，获取判定阈值；

5b)将局部方差图像中方差大于判定阈值的像素标记为被干扰像素；

(6)构造信号判定矩阵

利用下面公式构造信号判定矩阵：

$S(i,j)=\left\{\begin{array}{cc}1& D(i,j)\∉C\\ 0& D(i,j)\∈C\end{array}\right.$

其中，S(i，j)表示构造的信号判定矩阵，(i，j)表示信号判定矩阵第i行第j列对应的坐标值，D(i，j)表示被干扰图像在第i行第j列对应坐标的像素值，C表示被干扰像素的集合；

(7)修复被干扰像素，获取初步修复图像

利用如下公式修复被干扰图像中被γ射线干扰的像素，获取初步修复图像：

$p(i,j)=\frac{\underset{(m,n)\∈A}{\mathrm{\Σ}}S(m,n)D(m,n)}{\underset{(m,n)\∈A}{\mathrm{\Σ}}S(m,n)}$

其中，p(i，j)表示修复后的被干扰像素，(i，j)表示修复后的被干扰像素对应的坐标值，(m，n)表示以坐标(i，j)为中心的图像掩膜范围内的图像像素对应的坐标值，A表示以坐标(i，j)为中心的图像掩膜范围内的图像像素的集合，∑·表示求和运算符，S(m，n)表示信号判定矩阵第m行第n列对应的值，D(m，n)表示被干扰图像第m行第n列对应坐标的像素值；

(8)判断初步修复图像是否满足处理要求

人工观察判断初步修复图像中是否存在被干扰的像素，如果是，将初步修复图像作为被干扰图像，返回步骤(3)；否则，执行下一步骤；

(9)对初步修复图像抑噪

采用非截断小波抑噪公式，对初步修复图像进行抑噪，得到去除干扰图像，去除γ射线过程结束。

本发明与现有技术相比具有如下优点：

第一，本发明采用修复切伦科夫图像中被γ射线的像素点的方法，克服了现有技术中处理结果依赖于相机处理模块性能的缺点，使得普通CCD相机采集到的切伦科夫图像，采用本发明的方法处理后，在获得较好的图像质量的同时降低了切伦科夫成像成本。

第二，本发明在去除γ射线时，采用计算γ射线周围未被γ射线干扰的像素点灰度平均值的方法，去除了γ射线周围可能存在的其他γ射线的影响，克服了现有技术中处理低放射性活度条件下的切伦科夫荧光图像时图像质量下降的缺点，使得通过本发明处理后，图像信号能够尽可能的保留下来，能够获得更高的成像质量，即使对低放射性活度下采集到的切伦科夫荧光图像有较好的处理效果。

附图说明

图1为本发明的流程图；

图2为本发明的仿真示意图。

具体实施方式

下面结合图1对本发明的具体实施步骤做进一步的详细描述。

步骤1.获取图像数据

1a)在避光条件下用CCD相机捕获切伦科夫荧光图像；

1b)在避光条件下用CCD相机捕获背景噪声图像；

在本发明实施例中，调整CCD相机与光源的成像平面相距0.3m，CCD相机的曝光时间是180s，将带孔的塑料板四个圆孔中加入放射性活度分别是8uCi的Na-I131，其他圆孔不做处理。

步骤2.图像数据预处理

将采集的荧光图像减去背景噪声图像，得到被干扰图像，完成图像数据预处理。该处理能够消除背景噪声对以后处理结果的影响。

步骤3.根据如下公式计算图像掩膜的宽度

w＝2t+1

其中，w表示图像掩膜的宽度大小，t表示迭代的次数。使用这样的方法设置图像掩膜的宽度，使得图像每次迭代之后图像的掩膜就会增大，这样处理的结果就尽可能的避免了图像中产生局部平均斑块，使得处理后图像较平滑。

步骤4.获取局部方差图像

使用平滑线性滤波器，对被干扰图像进行平滑线性滤波，获得图像中每个像素的局部均值。

利用图像局部方差公式，计算被干扰图像中每个像素的局部方差，获取局部方差图像。计算图像的局部方差能够提高图像的对比度，图像中的孤立像素点的局部方差较大，而图像中平滑的区域的局部方差就较小，从而便于区分图像中被γ射线干扰的像素。图像局部方差按如下公式计算得来：

$V(i,j)=\underset{(k,l)\∈A}{\mathrm{\Σ}}{[D(k,l)-M(k,l)]}^2$

其中，V(i，j)表示以坐标(i，j)为中心图像掩膜范围内图像的局部方差，(k，l)表示以坐标(i，j)为中心的图像掩膜范围内的图像像素的坐标值，A表示以坐标(i，j)为中心的图像掩膜范围内的图像像素的集合，∑·表示求和运算符，D(k，l)表示被干扰图像第k行第l列对应坐标的像素值，M(k，l)表示以像素坐标(k，l)为中心的图像掩膜范围内图像的局部均值。

步骤5.确定被干扰的像素

采用最大类间方差法处理局部方差图像，获取判定阈值。判定阈值可以将局部方差图像的像素点分为两类，一类是被γ射线干扰的像素，一类是正常的切伦科夫图像的像素。最大类间方差法按如下步骤进行：

第一步，任意选定局部方差图像中一个像素点，将该像素点的灰度值作为初始阈值；

第二步，利用初始阈值将局部方差图像分为背景和前景两部分；

第三步，按照下式计算局部方差图像的类间方差：

g＝w₀w₁(u₀-u₁)

其中，g表示局部方差图像的类间方差，w₀表示局部方差图像前景部分像素点总数占整幅图像的比值，w₁表示局部方差图像背景部分像素点总数占整幅图像的比值，u₀表示局部方差图像前景部分所有像素的灰度平均值，u₁表示局部方差图像背景部分所有像素的灰度平均值；

第四步，遍历局部方差图像中所有的灰度值，得到使局部方差图像的类间方差最大的阈值。

将局部方差图像中方差大于判定阈值的像素标记为被干扰像素。在判定结束后，得到被干扰像素的集合，在以下的步骤中只对被干扰像素进行替换处理。

步骤6.构造信号判定矩阵

利用下面公式构造信号判定矩阵：

$S(i,j)=\left\{\begin{array}{cc}1& D(i,j)\∉C\\ 0& D(i,j)\∈C\end{array}\right.$

其中，S(i，j)表示构造的信号判定矩阵，(i，j)表示信号判定矩阵第i行第j列对应的坐标值，D(i，j)表示被干扰图像在第i行第j列对应坐标的像素值，C表示被干扰像素的集合。信号判定矩阵中的坐标与被干扰图像中像素的坐标是对应的，当信号判定矩阵中某个坐标的值是1，被干扰图像中对应坐标的像素即是被γ射线干扰的像素；否则就是未被干扰的像素。

步骤7.修复被干扰像素，获取初步修复图像

利用如下公式修复被干扰图像中被γ射线干扰的像素，获取初步修复图像：

$p(i,j)=\frac{\underset{(m,n)\∈A}{\mathrm{\Σ}}S(m,n)D(m,n)}{\underset{(m,n)\∈A}{\mathrm{\Σ}}S(m,n)}$

其中，p(i，j)表示修复后的被干扰像素，(i，j)表示修复后的被干扰像素对应的坐标值，(m，n)表示以坐标(i，j)为中心的图像掩膜范围内的图像像素对应的坐标值，A表示以坐标(i，j)为中心的图像掩膜范围内的图像像素的集合，∑·表示求和运算符，S(m，n)表示信号判定矩阵第m行第n列对应的值，D(m，n)表示被干扰图像第m行第n列对应坐标的像素值。

利用该公式处理被干扰图像，由于使用了信号判定矩阵，只计算被干扰像素邻域周围未被干扰像素的灰度平均值，利用该灰度平均值替换被干扰像素的灰度值。这样基本完全排除了γ射线对初步修复图像灰度造成的影响，虽然不能够精确的恢复该像素点的原始值，但是根据图像的平滑性，用邻域的均值替换是合理的。

步骤8.判断初步修复图像是否满足处理要求

人工观察判断初步修复图像中是否存在被干扰的像素，如果是，将初步修复图像作为被干扰图像，返回步骤(3)；否则，执行下一步骤。

采用人工观察方法判断每次迭代后得到的初步修复图像，因为仅用一次处理不能保证有很好的处理效果，采用迭代的方法对图像多次处理，不仅去除γ射线的效果较好，而且处理后的图像更平滑。在现有阶段还没有对γ射线干扰定量的评价标准，不能够直接评价去除γ射线干扰的效果，本发明采用人工观察的方式对处理效果的评价。

步骤9.对初步修复图像抑噪

采用非截断小波抑噪公式，对初步修复图像进行抑噪，得到修复图像，修复结束。

由于初步修复图像在采集时和处理后都有可能产生噪声，利用非截断小波对图像去噪能够进一步的提高修复图像的质量。非截断小波抑噪公式如下：

$F=\mathrm{R\θ}(\mathrm{Dp},\stackrel{\‾}{D}p)$

其中，F表示抑噪后的图像，R表示非截断小波变换重建矩阵，θ表示逐点门限函数组，〔·〕表示矩阵合并运算符，D表示非截断小波变换分解矩阵，p表示初步修复图像，表示非截断小波变换分解矩阵的共轭矩阵。

下面结合附图2对本发明的去除γ射线结果做进一步的描述。

1.仿真条件

在CPU为core 2 2.3GHZ、内存2G、Window7操作系统上使用Matlab 2009a进行仿真。

2.仿真内容

对于切伦科夫荧光图像，现在还没有定量评价标准，本发明的效果只能从视觉效果上进行比较，即从γ射线亮斑的多少，切伦科夫荧光信号边界的清晰度进行比较。图2所示即为本发明仿真比较示意图，其中图2(a)为用现有技术中CCD相机自身配备的去除宇宙射线模块处理的切伦科夫荧光图像；图2(b)为使用本方法处理的切伦科夫荧光图像。在图2(b)所示的切伦科夫荧光图像中，四个圆形白色区域即是切伦科夫荧光信号，在图像中白色亮斑即是γ射线信号。

3.仿真结果

从图2b)中可见，本发明处理的图像中γ射线的光斑基本完全消除，切伦科夫荧光信号的边缘较清晰，从图2a)中可见，CCD相机自身配备的宇宙射线模块处理的图像中明显有大量γ射线光斑存在，切伦科夫荧光信号的边缘也较模糊。

从以上仿真效果可以看出，采用切伦科夫荧光成像中γ射线的去除方法，能够基本完全去除切伦科夫荧光图像中的γ射线，切伦科夫荧光信号的边缘较清晰，从而提高了切伦科夫荧光成像的图像质量。

Claims (4)

1.切伦科夫荧光成像中γ射线的去除方法，包括如下步骤：

(1)获取图像数据

1a)在避光条件下用CCD相机捕获切伦科夫荧光图像；

1b)在避光条件下用CCD相机捕获背景噪声图像；

(2)图像预处理

将采集的荧光图像减去背景噪声图像，得到被干扰图像；

(3)根据如下公式计算图像掩膜的宽度：

w＝2t+1

其中，w表示图像掩膜的宽度大小，t表示迭代的次数；

(4)获取局部方差图像

4a)使用平滑线性滤波器，对被干扰图像进行平滑线性滤波，获得图像中每个像素的局部均值；

4b)利用图像局部方差公式，计算被干扰图像中每个像素的局部方差，获取局部方差图像；

(5)确定被干扰像素

5a)采用最大类间方差法处理局部方差图像，获取判定阈值；

5b)将局部方差图像中方差大于判定阈值的像素标记为被干扰像素；

(6)构造信号判定矩阵

利用下面公式构造信号判定矩阵：

$S(i,j)=\left\{\begin{array}{cc}1& D(i,j)\∉C\\ 0& D(i,j)\∈C\end{array}\right.$

其中，S(i，j)表示构造的信号判定矩阵，(i，j)表示信号判定矩阵第i行第j列对应的坐标值，D(i，j)表示被干扰图像在第i行第j列对应坐标的像素值，C表示被干扰像素的集合；

(7)修复被干扰像素，获取初步修复图像

利用如下公式修复被干扰图像中被γ射线干扰的像素，获取初步修复图像：

$p(i,j)=\frac{\underset{(m,n)\∈A}{\mathrm{\Σ}}S(m,n)D(m,n)}{\underset{(m,n)\∈A}{\mathrm{\Σ}}S(m,n)}$

其中，p(i，j)表示修复后的被干扰像素，(i，j)表示修复后的被干扰像素对应的坐标值，(m，n)表示以坐标(i，j)为中心的图像掩膜范围内的图像像素对应的坐标值，A表示以坐标(i，j)为中心的图像掩膜范围内的图像像素的集合，∑·表示求和运算符，S(m，n)表示信号判定矩阵第m行第n列对应的值，D(m，n)表示被干扰图像第m行第n列对应坐标的像素值；

(8)判断初步修复图像是否满足处理要求

人工观察判断初步修复图像中是否存在被干扰的像素，如果是，将初步修复图像作为被干扰图像，返回步骤(3)；否则，执行下一步骤；

(9)对初步修复图像抑噪

采用非截断小波抑噪公式，对初步修复图像进行抑噪，得到去除干扰图像，去除γ射线的过程结束。

2.根据权利要求1中所述的切伦科夫荧光成像中γ射线的去除方法，其特征在于，步骤4b)所述的图像的局部方差根据如下公式计算得来：

$V(i,j)=\underset{(k,l)\∈A}{\mathrm{\Σ}}{[D(k,l)-M(k,l)]}^2$

其中，V(i，j)表示以坐标(i，j)为中心图像掩膜范围内图像的局部方差，(k，l)表示以坐标(i，j)为中心的图像掩膜范围内的图像像素的坐标值，A表示以坐标(i，j)为中心的图像掩膜范围内的图像像素的集合，∑·表示求和运算符，D(k，l)表示被干扰图像第k行第l列对应坐标的像素值，M(k，l)表示以像素坐标(k，l)为中心的图像掩膜范围内图像的局部均值。

3.根据权利要求1中所述的切伦科夫荧光成像中γ射线的去除方法，其特征在于，步骤5a)所述的最大类间方差法，按如下步骤进行：

第一步，任意选定局部方差图像中一个像素点，将该像素点的灰度值作为初始阈值；

第二步，利用初始阈值将局部方差图像分为背景和前景两部分；

第三步，按照下式计算局部方差图像的类间方差：

g＝w₀w₁(u₀-u₁)

其中，g表示局部方差图像的类间方差，w₀表示局部方差图像前景部分像素点总数占整幅图像的比值，w₁表示局部方差图像背景部分像素点总数占整幅图像的比值，u₀表示局部方差图像前景部分所有像素的灰度平均值，u₁表示局部方差图像背景部分所有像素的灰度平均值；

第四步，遍历局部方差图像中所有的灰度值，得到使局部方差图像的类间方差最大的阈值。

4.根据权利要求1中所述的切伦科夫荧光成像中γ射线的去除方法，其特征在于，步骤(9)所述的非截断小波抑噪公式如下：

$F=\mathrm{R\θ}(\mathrm{Dp},\stackrel{\‾}{D}p)$

其中，F表示抑噪后的图像，R表示非截断小波变换重建矩阵，θ表示逐点门限函数组，〔·〕表示矩阵合并运算符，D表示非截断小波变换分解矩阵，p表示初步修复图像，

表示非截断小波变换分解矩阵的共轭矩阵。

Images