CN103679706A - 一种基于图像各向异性边缘检测的ct稀疏角度重建方法 - Google Patents

Abstract

本发明提供一种基于图像各向异性边缘检测的CT稀疏角度重建方法，该方法基于图像各向异性的边缘检测信息，在图像重建过程中更新TV各加权分量的加权权重，从而达到增强图像边缘，提高重建图像的分辨率，改善重建图像质量的目的。本发明可以在更少的投影幅数下重建图像，更加有效的降低CT成像的辐射剂量。

Description

一种基于图像各向异性边缘检测的CT稀疏角度重建方法

技术领域

本发明针对计算机断层扫描仪的低剂量扫描和重建，主要涉及数字图像处理和CT图像重建方法，具体涉及一种基于图像各向异性边缘检测的CT稀疏角度重建方法。

背景技术

1967年Hounsfield发明了世界上第一台计算机断层扫描(ComputedTomography,CT)设备。CT技术通过采集X射线穿过物体后多个角度的投影数据，利用计算机重建图像，进而获得物体内部结构的断层图像。由于CT设备利用X射线的穿透性，在不破坏物体的情况下可检测物体的内部结构，至今已广泛应用于医学诊断、工业检测、考古研究、珠宝鉴定等多个领域。

然而，医学研究表明，X射线对人体有损伤，电离辐射可从分子水平和细胞水平对生物体造成损伤，照射X光可诱发癌症、白血病或其他遗传性疾病。而医学CT是利用X射线进行的一种高辐射剂量的检查，其辐射危害对于胎儿和儿童尤为严重。据新英格兰医学杂志2009年针对95万人的统计显示，在所有能够产生辐射的医学成像设备中，检查量仅占总数21%的CT及核医学成像，产生的辐射量却占总量的75%以上，特别是CT，产生的辐射量接近总量的50%。因此如何在保持图像质量的前提下有效的降低辐射剂量引起越来越多的关注，国外厂商西门子、GE等都在积极的研究可能的方案。

提高探测器检测效率或使用准直器等硬件方案以减少辐射剂量非常有效，但是会增加额外的成本，且某些方案技术难度较大。在不改变硬件配置的前提下，降低辐射剂量通常采用的方式是减少投影角度或降低电流。前者面临着CT投影数据的不完全和稀疏角度投影的重建困难，后者会带来图像的高噪声。为改善图像质量，研究者针对两种情况分别提出了很多改进的算法。其中利用压缩感知(Compressed Sensing,CS)理论，结合人体组织结构分片连续的特点，可较好的实现稀疏角度的重建。压缩感知理论指出，只要信号是可压缩的或在某个变换域是稀疏的，就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，通过求解一个优化问题从这个低维信号的少量投影中以高概率重构出原信号，该投影包含了重构信号的足够信息，并在理论上证明了该技术的错误率是无穷小，从而奠定了CS技术的理论基础，引发了采样技术的革命性变化，成为近年来最热门的研究技术之一。CS的思想是从尽量少的数据中提取尽量多的信息，从而实现信号在较低采样率下的采集和复原，其前提条件是图像本身是稀疏的或者图像是可压缩的。

CT图像经过梯度变换在不同程度上满足稀疏化条件。根据CS理论，重构稀疏信号采用求解其l₀范数最小化的方法，而实际中往往转换为求其l₁范数最小化，研究者提出采用CT图像梯度向量的模的l₁范数（即CT图像的总变差Total Variation，TV）最小化方法重建CT图像。传统的TV计算采用图像像素梯度向量的模相加得到，而梯度向量的模为图像像素沿x和y两个方向的一阶导数(即梯度的分量)的平方和相加后的平方根（式(1)），其在x和y方向的一阶导数的权重因子是相同的（即各向同性），均为1，未考虑图像梯度的方向性，因此在图像求解过程中不可避免的带来图像边缘的模糊和分辨率的下降。

${\stackrel{\→}{u}}_{\mathrm{TV}}=\underset{i,j}{\mathrm{\Σ}}\sqrt{{(u_{i,j}-u_{i-1,j})}^2+{(u_{i,j}-u_{i,j-1})}^2}---\left(1\right)$

中国专利CN 103136773A提出的“一种稀疏角度X射线CT成像方法”中对首次扫描的全角度数据重建图像与当前稀疏角度重建的图像采用加权平均滤波处理获得先验图像，由此先验图像构建稀疏角度CT图像重建模型，再对稀疏角度重建的图像优化求解获得最终的重建图像。此专利虽然可以消除由于先前扫描的CT图像与当前重建的CT图像之间因成像位置不匹配引起的运动伪影，但是需要有全角度数据重建的先验图像用以构建稀疏角度重建模型，且其TV计算仍然采用传统的梯度的相加的方法。

如上指出的，传统的TV最小化方法用于CT稀疏角度重建时，采用梯度分量等权重加和的方法计算各像素梯度的模，是一种各向同性的梯度模稀疏化方法，并且重建过程中梯度分量的加权权重保持不变，会导致重建图像的边缘模糊，影响图像中较小组织结构的分辨。

发明内容

本发明的目的在于提供一种基于图像各向异性边缘检测的CT稀疏角度重建方法。

为达到上述目的，本发明采用了以下技术方案。

该CT稀疏角度重建方法将CT图像特征与CT图像重建相结合，利用图像各向异性边缘检测的边缘信息改变图像总变差计算中各分量的加权权重，并采用图像总变差最小化的方法重建图像。

所述CT稀疏角度重建方法具体包括以下步骤：

1）设置初始迭代次数k=0，初始化重建图像u(0)，图像大小为L×L；同时，初始化图像总变差计算中用于各向异性加权的权重矩阵g；

2）增加迭代次数k=k+1，利用代数迭代法重建得到更新的重建图像；

3）根据步骤2）中得到的更新的重建图像，求解满足图像总变差最小化的图像u(k)；

4）对图像u(k)进行各向异性的边缘检测，得到x和y两个方向的边缘图像EdgeX以及EdgeY；

5）根据边缘图像EdgeX与EdgeY更新权重矩阵为g_new；

6）使用更新后的权重矩阵g_new计算图像u(k)的图像总变差；

重复步骤2）至步骤6）直至满足迭代截止条件，并输出图像u(k)。

所述权重矩阵g对于四邻域时为L×L*4的三维矩阵，每个像素点在图像总变差计算中对应有四个权重，权重初始值均设为π/4。

所述图像总变差的计算方法为：

${\stackrel{\→}{u}}_{\mathrm{TV}}^{\mathrm{aniso}}=\underset{\mathrm{\α}}{\mathrm{\Σ}}{g}_{\mathrm{\α}}\left(D_{\mathrm{\α}}\stackrel{\→}{u}\right)=\mathrm{\Σ}{g}_{(i,j)\_(s,t)}|u_{i,j}-u_{s,t}|$

其中，

表示各向异性加权的图像总变差，g_{(i,j)_(s,t)}表示

中各个分量的权重，(s,t)为像素点(i,j)的邻域像素，

计算中分量的个数由像素点的邻域α决定，u_i,j表示图像u中像素点(i,j)处的像素值，u_s,t表示图像u中像素点(s,t)处的像素值，

表示图像u沿像素的邻域方向的导数的绝对值。

所述边缘检测的具体步骤为：

a）对图像进行高斯滤波；

b）对滤波后的图像分别沿x、y方向求导，得到对应的梯度分量图像Dx与Dy；

c）根据梯度分量图像Dx与Dy将图像像素点(i,j)按照其梯度矢量（Dx_ij,Dy_ij）所属区域分为x向边缘点以及y向边缘点；x向边缘指边缘切线方向偏向x轴，y向边缘指边缘切线方向偏向y轴；

d）对步骤c）分类后的边缘点沿梯度方向进行梯度模的非极大值抑制，选择梯度模的局部最大值点位置；

e）对步骤d）找到的梯度模的局部最大值点的梯度模进行高、低双阈值约束，从而找到大于等于高阈值的点域Ω和位于高、低阈值之间的点域Ψ；

f）对Ω中的点以Ψ中的点作为联通邻域补充进行形态学连接和边缘细化。

权重矩阵g_new根据边缘图像EdgeX与EdgeY得到，当EdgeX中像素点(i,j)的值为逻辑真时，更新权重g_{(i,j)_(i+1,j)}以及g_{(i+1,j)_(i,j)}为零，否则保持原值不变；同样，当EdgeY中像素点(i,j)的值为逻辑真时，更新权重g_{(i,j)_(i,j+1)}以及g_{(i,j+1)_(i,j})为零，否则保持原值不变。

本发明的有益效果体现在：

本发明将图像本身的边缘特性与CT图像的重建相结合，基于图像各向异性的边缘检测信息更新TV计算中各分量的权重，用于后续的最优化计算过程，从而达到增强图像边缘，提高重建图像的分辨率，改善重建图像质量的目的，本发明可以在更少的投影幅数下重建图像，可以更加有效的降低CT成像的辐射剂量。

附图说明

图1为各向异性边缘检测导引的TV最小化重建流程图；

图2为四邻域和八邻域示意图，其中，(a)为像素点(i,j)的四邻域表示，(i-1,j),(i+1,j),(i,j-1),(i,j+1)为像素点(i,j)的四邻域；(b)为像素点(i,j)的八邻域表示，其中(i-1,j),(i+1,j)(i,j-1)(i,j+1)为(i,j)最近邻的八邻域，(i-1,j-1),(i-1,j+1)(i+1,j-1)(i+1,j+1)为(i,j)较远的八邻域；

图3为边缘方向的划分方法，像素点梯度矢量位于1、1’、4、4’区域时其边缘属于x向；像素点梯度矢量位于2、2’、3、3’区域时其边缘属于y向边缘；

图4为各向异性边缘检测结果示意图，其中，(a)为原始图像，(b)为各向异性边缘检测后边缘属于x向的边缘，(c)为各向异性边缘检测后边缘属于y向的边缘；

图5为修正的shepp-logan的15副投影重建图，其中，(a)为被投影原始图像，(b)为传统TV最小化方法重建图，(c)为本发明中方法重建图，其中黑色框中区域放大图置于各图的右下角；

图6为修正的shepp-logan的20副投影加噪重建图，其中，(a)为传统TV最小化方法重建图，(b)为本发明中方法重建图，其中黑色框中区域放大图置于各图的右下角。

具体实施方式

下面结合附图对本发明作详细说明。

针对由于图像总变差采用梯度模加和计算带来的边缘模糊，本发明提供一种利用图像各向异性边缘检测的边缘信息改变TV各加权分量的加权权重的CT稀疏角度重建方法，简称为各向异性边缘检测导引的TV最小化重建，本发明可以增强图像边缘，提高重建图像的分辨率，同时可以在更少的角度下重建图像。

图像边缘是图像稀疏化的一种直观表示，也反映了图像的大部分信息。而图像边缘的切线方向平行于图像对x和y的一阶导数的比值，反映了图像梯度向量的方向信息。本发明将图像本身的边缘特性与CT图像重建相结合，基于图像各向异性的边缘检测信息，改变TV计算中各分量的权重，用于后续的最优化计算过程。本发明为CT重建方法和稀疏角度重建提供了一种新的思路，其基本步骤如下：

（1）获取CT稀疏角度投影数据p；

（2）对投影数据进行稀疏角度重建，具体步骤为：

a）设置初始迭代次数k=0，初始化重建图像u(0)，图像大小为L×L；

b）初始化图像总变差计算中用于各向异性加权的权重矩阵g；

c）增加迭代次数k=k+1；利用代数迭代法重建得到更新的重建图像；

d）利用步骤c）的结果求解满足各向异性加权的图像TV最小化的图像u(k)；

e）对u(k)进行各向异性的边缘检测，得到x和y两个方向的边缘图像EdgeX以及EdgeY；

f）根据边缘图像EdgeX与EdgeY更新权重矩阵为g_new；

g）使用权重矩阵g_new计算图像u(k)的各向异性加权的图像总变差；

重复步骤c)-g)直至满足迭代截止条件，并输出图像u(k)。

本发明的具体实施步骤如下：

首先获取CT稀疏角度的投影数据p。具体可以通过计算机仿真或者CT设备采集得到。可以直接采集360°范围内均匀分布的稀疏角度投影，也可以对正常扫描的投影均匀采样得到。投影幅数根据噪声的实际情况设定，实验扫描数据往往减少为正常投影数的1/4或者更少，计算机仿真在不含噪声的情况下可以减至20副投影以下。另外重建所需投影幅数也和被扫描对象的图像灰度分布有关系。p为M×N的矩阵，其中，M为稀疏角度投影幅数，N为探测器单元数。

其次利用循环迭代的思想对投影数据进行稀疏角度重建。其重建流程图如图1所示。首先设定初始的二维图像u(0)，在没有额外先验图像的条件下，可以将u(0)设为零矩阵，矩阵大小为L×L。传统的图像总变差计算公式如式(1)，本发明按照式(2)计算各向异性加权的图像总变差

式(2)中，g_{(i,j)_(s,t)}代表

中各个分量的权重因子，(s,t)为像素点(i,j)的邻域像素。

计算中分量的个数由像素点的邻域α决定，u_i,j表示图像u中像素点(i,j)处的像素值，u_s,t表示图像u中像素点(s,t)处的像素值，

表示图像u沿像素的邻域方向的导数的绝对值。本发明采用四邻域计算，同时给出图像像素点的八邻域表示(图2)。在四邻域中，权重初始值设为π/4，八邻域中，最近邻的四邻域权重初始值设为π/8，较远的四邻域权重初始值设为

本发明中权重矩阵g构造为L×L*4的三维矩阵，计算

时每个像素点有四个权重因子和其对应。

${\stackrel{\→}{u}}_{\mathrm{TV}}^{\mathrm{aniso}}=\underset{\mathrm{\α}}{\mathrm{\Σ}}{g}_{\mathrm{\α}}\left(D_{\mathrm{\α}}\stackrel{\→}{u}\right)=\mathrm{\Σ}{g}_{(i,j)\_(s,t)}|u_{i,j}-u_{s,t}|---\left(2\right)$

求解式(3)所表达的最优化问题。式(3)条件约束部分表达了重建图像需满足CT投影方程和图像像素值的非负性。式中H为CT投影的系统矩阵，描述了每个图像像素点对各个探测器单元的贡献权重。式(3)具体求解步骤为：采用代数迭代方法得到本步骤中的初始图像，然后计算其各向异性加权的图像总变差，通过求解其最小值重建新的图像，并使其满足非负约束，此时重建图像记做u(k)。本发明中采用最速下降法求解TV最小化问题，也可以采用其它的优化方法。

$\begin{array}{cc}\underset{\stackrel{\→}{u}}{\min }\left({\stackrel{\→}{u}}_{\mathrm{TV}}^{\mathrm{aniso}}\right)& s.t.H\stackrel{\→}{u}=p,u_{i,j}\≥0---\left(5\right)\end{array}$

对u(k)进行各向异性的边缘检测。边缘检测采用改进的亚像素级canny检测方法。具体如下：

a）对图像进行高斯滤波，高斯核中的标准偏差随噪声大小决定，当噪声较大时，增大标准偏差；否则减少标准偏差。默认值取1.0。

b）对滤波后的图像分别沿x、y方向求导，得到梯度分量图像Dx与Dy。

$\mathrm{Dx}=\∂\stackrel{\→}{u}/\∂x,\mathrm{Dy}=\∂\stackrel{\→}{u}/\∂y---\left(4\right)$

c）根据梯度分量图像Dx与Dy将图像像素点(i,j)按照其梯度矢量（Dx_ij,Dy_ij）所属区域分为x向边缘（即边缘切线方向偏向x轴）点、y向边缘（边缘切线方向偏向y轴）点。如图3所示，当像素点梯度方向位于1、1’、4、4’区域时其边缘属于x向；而2、2’、3、3’属于y向边缘。

d）对以上分类后的边缘点沿梯度方向进行梯度模的非极大值抑制，选择梯度模的局部最大值点位置。

对步骤d）中找到的点的梯度模进行高低双阈值约束，高阈值默认为非边缘点占总像素点70%时对应的梯度模值。具体实现过程为根据图像梯度模值对应的直方图累加，累加数目达到总像素数目的70%时，对应的图像梯度模值设置为高阈值，70%为默认值，可以根据实际情况改变。低阈值默认为高阈值的40%，也可以根据实际情况改变。从中找到大于等于高阈值的点域Ω和位于高低阈值之间的点域Ψ（所述之间表示小于高阈值，且大于等于低阈值）。

对Ω中的点以Ψ中的点作为联通邻域补充进行形态学连接和边缘细化。最终得到沿x向和y向的两幅边缘图像EdgeX，EdgeY。对一个圆和矩形各向异性边缘检测后如图4所示。

根据图像边缘信息更新权重矩阵g。当EdgeX中像素点(i,j)的值为逻辑真时，更新权重g_{(i,j)_(i+1,j})以及g_{(i+1,j)_(i,j)}为零，否则保持原值不变；同样，当EdgeY中像素点(i,j)的值为逻辑真时，更新权重g_{(i,j)_(i,j+1)}以及g_{(i,j+1)_(i,j)}为零，否则保持原值不变。

重复以上步骤，求解式(3)，直至满足迭代终止条件。本发明中采用迭代次数终止迭代过程。一般迭代两次图像就可以有明显的边缘增强效果，对于噪声较强的图像可以增加迭代次数，但是一般不超过5次。

下面通过计算机仿真实验说明本发明方法的有效性。对修正的256×256大小的shepp-logan模型(图5(a))进行计算机仿真投影，其中源到探测器的距离为800mm，源到旋转中心的距离为400mm。探测器像素大小为0.806641mm，512个探测单元。此时视野半径为200mm。共产生360副投影数据，从中均匀选取15副投影作为重建所需投影数据p，其矩阵大小为20×512。传统TV重建如图5(b)，按照本发明方法重建图像如图5(c)所示。

对投影均匀抽取20个投影角度，并加上0.1%的高斯噪声，再做稀疏角度重建。重建结果如图6所示，其中图6(a)为传统方法重建，图6(b)为本发明方法重建效果。

由以上实例重建效果可以看出，本发明所提出的方法可以更好的恢复图像的边缘，提高图像分辨率，且在少角度下可以得到更好的图像质量。

Claims (6)

1.一种基于图像各向异性边缘检测的CT稀疏角度重建方法，其特征在于：该CT稀疏角度重建方法包括以下步骤：将CT图像特征与CT图像重建相结合，利用图像各向异性边缘检测的边缘信息改变图像总变差计算中各分量的加权权重，并采用图像总变差最小化的方法重建图像。

2.根据权利要求1所述一种基于图像各向异性边缘检测的CT稀疏角度重建方法，其特征在于：所述CT稀疏角度重建方法具体包括以下步骤：

1）设置初始迭代次数k=0，初始化重建图像u(0)，图像大小为L×L；同时，初始化图像总变差计算中用于各向异性加权的权重矩阵g；

2）增加迭代次数k=k+1，利用代数迭代法重建得到更新的重建图像；

3）根据步骤2）中得到的更新的重建图像，求解满足图像总变差最小化的图像u(k)；

4）对图像u(k)进行各向异性的边缘检测，得到x和y两个方向的边缘图像EdgeX以及EdgeY；

5）根据边缘图像EdgeX与EdgeY更新权重矩阵为g_new；

6）使用更新后的权重矩阵g_new计算图像u(k)的图像总变差；

重复步骤2）至步骤6）直至满足迭代截止条件，并输出图像u(k)。

3.根据权利要求2所述一种基于图像各向异性边缘检测的CT稀疏角度重建方法，其特征在于：所述权重矩阵g对于四邻域时为L×L*4的三维矩阵，每个像素点在图像总变差计算中对应有四个权重，权重初始值均设为π/4。

4.根据权利要求1或2所述一种基于图像各向异性边缘检测的CT稀疏角度重建方法，其特征在于：所述图像总变差的计算方法为：

${\stackrel{\→}{u}}_{\mathrm{TV}}^{\mathrm{aniso}}=\underset{\mathrm{\α}}{\mathrm{\Σ}}{g}_{\mathrm{\α}}\left(D_{\mathrm{\α}}\stackrel{\→}{u}\right)=\mathrm{\Σ}{g}_{(i,j)\_(s,t)}|u_{i,j}-u_{s,t}|$

其中，

表示各向异性加权的图像总变差，g_{(i,j)_(s,t)}表示中各个分量的权重，(s,t)为像素点(i,j)的邻域像素，

计算中分量的个数由像素点的邻域α决定，u_i，j表示图像u中像素点(i,j)处的像素值，u_s，t表示图像u中像素点(s,t)处的像素值，表示图像u沿像素的邻域方向的导数的绝对值。

5.根据权利要求1或2所述一种基于图像各向异性边缘检测的CT稀疏角度重建方法，其特征在于：所述边缘检测的具体步骤为：

a）对图像进行高斯滤波；

b）对滤波后的图像分别沿x、y方向求导，得到对应的梯度分量图像Dx与Dy；

c）根据梯度分量图像Dx与Dy将图像像素点(i,j)按照其梯度矢量（Dx_ij,Dy_ij）所属区域分为x向边缘点以及y向边缘点；x向边缘指边缘切线方向偏向x轴，y向边缘指边缘切线方向偏向y轴；

d）对步骤c）分类后的边缘点沿梯度方向进行梯度模的非极大值抑制，选择梯度模的局部最大值点位置；

e）对步骤d）找到的梯度模的局部最大值点的梯度模进行高、低双阈值约束，从而找到大于等于高阈值的点域Ω和位于高、低阈值之间的点域Ψ；

f）对Ω中的点以Ψ中的点作为联通邻域补充进行形态学连接和边缘细化。

6.根据权利要求2所述一种基于图像各向异性边缘检测的CT稀疏角度重建方法，其特征在于：权重矩阵g_new根据边缘图像EdgeX与EdgeY得到，当EdgeX中像素点(i,j)的值为逻辑真时，更新权重g_{(i,j)_(i+1,j)}以及g_{(i+1,j)_(i,j})为零，否则保持原值不变；同样，当EdgeY中像素点(i,j)的值为逻辑真时，更新权重g_{(i,j)_(i,j+1)}以及g_{(i,j+1)_(i,j})为零，否则保持原值不变。

Images