CN104361618A - 基于分形和压缩感知的ct重建方法 - Google Patents

Abstract

一种基于分形和压缩感知的CT重建方法，包括步骤：(1)已知CT投影数据b和相应的投影角度θ，初始化参数β、λ、λ_n，CT图像u＝0；(2)根据投影角度θ，计算出相应的投影矩阵A；(3)使用式(5)计算CT图像u_ART；(4)对u_ART进行分形编码处理得到Φ(u)，使用式(12)得到α；(5)对Φ(u)做解码处理得到Φ^T(u)；6)分别求解Φ^T(u)Φ(u)、Φ^T(u)α、A^TA、A^Tb；(7)使用共轭梯度法求解式(15)，得到CT图像；8)检查是否满足迭代结束的条件，是转至步骤(9)，否则步骤(3)；(9)结束，输出CT图像。方法使用分形作为稀疏变换以获得更稀疏的系数表示，使用迭代连续迭代算法进行求解，实现在少投影的情况下重建出高质量的CT图像。

Description

基于分形和压缩感知的CT重建方法

技术领域

本发明涉及断层影像重建技术领域，特别是用于少剂量CT重建算法。

背景技术

自1973年问世以来，CT(computed tomography)广泛应用于医疗诊断、工业无损检测等领域，尤其是在医学诊断中发挥着举足轻重的作用。然而CT所使用的光源是对人体有害的X光，有可能引发基因突变、癌症等疾病，所以尽可能减少X光辐射剂量势在必行。减少X光辐射剂量有两个基本思路，第一，减少X光的强度。第二，较少X光投影的个数即减少投影个数。采用第一种方法会导致CT图像的信噪比降低，第二种方法重建出的图像将出现较多的伪影，但随着压缩感知理论广泛应用于CT图像重建中，人们往往采用第二种方法减少X光辐射剂量。即在少投影的条件下重建出质量高的CT图像，这也是CT领域中的研究热点。

2006年，Candes、Tao和Donoho等人提出的压缩感知理论的主要思想是大部分信号在适当的正交变换域中，如梯度变换、小波变换下是“稀疏”的，即信号的大部分系数接近于0或者等于0。对可压缩的(稀疏)信号在远低于奈奎斯特采样频率进行数据采样后，仍能够精确恢复原始信号。而CT的采样信号则是典型的“稀疏”信号，因为在检测物体(如人体)的相邻区域内，X射线的衰减系数相近或者相等。所以将压缩感知技术应用于CT领域，使得在“稀疏”投影的条件下也能重建出质量较高的图像。

在压缩感知技术中很重要的一环便是稀疏变换，传统的稀疏变换如梯度变换、小波变换都不能最稀疏的、很好的表示图像，而分形则可以更稀疏、更好的表示图像。

分形的思想基础是基于大自然界的物体图像的局部和全局或者局部和局部具有极强的自相似性，如树木、河流、湖泊及海岸线等都具有极强的自相似性。所以可以利用这种自相似性使用较少的数据量描述图像，再通过迭代函数系统重建出图像。使用分形的图像压缩技术可以获得比传统压缩方法更高的压缩比，如使用小波的图像压缩方法可以到达1：20的压缩比，而采用分形的图像压缩方法可以获得1：100甚至更高的压缩比，个别图像甚至可以达到1：10000的压缩比。所以采用分形作为压缩感知的稀疏变换则可以获得比传统稀疏变换更稀疏的稀疏变换。

发明内容

本发明为了解决减少CT扫描的辐射剂量即为基于少投影的CT重建，提出基于分形和压缩感知的CT重建算法，使用分形作为稀疏变换以获得更稀疏的系数表示，使用迭代连续迭代算法(Alternating and Continuation Algorithm，ACA)进行求解，以实现在少投影的情况下重建出高质量的CT图像。

本发明的技术方案如下：

一种基于分形和压缩感知的CT重建方法，所述方法使用分形作为稀疏变换以获得更稀疏的系数表示，使用迭代连续迭代算法进行求解，重建出高质量的CT图像，步骤如下：

(1)已知CT投影数据b和相应的投影角度θ，初始化公式(12)中的阈值参数β，压缩感知中的正则化参数λ、ART算法的迭代参数λ_n，初始化CT图像为0即u＝0；

(2)根据投影角度θ，计算出相应的投影矩阵A；

(3)使用式公式(5)计算CT图像u_ART

${u_j}^{(n+1)}={u_j}^{\left(n\right)}+{\mathrm{\λ}}_n\frac{a_{\mathrm{ij}}}{{\left|\right|A^i\left|\right|}^2}(b^i-A^i{u_j}^{\left(n\right)}),n=\mathrm{0,1},....---\left(5\right)$

(4)对u_ART进行分形编码处理得到Φ(u)，使用式(12)得到α

$\mathrm{\α}=\max \left(\right|\mathrm{\Φ}\left(u\right)|-\frac{1}{\mathrm{\β}},0)\frac{\mathrm{\Φ}\left(u\right)}{\left|\mathrm{\Φ}\left(u\right)\right|}---\left(12\right)$

其中Φ(u)为分形稀疏变换，α为软阈值方法的求解值；

(5)对Φ(u)做解码处理得到Φ^T(u)

(6)分别求解Φ^T(u)Φ(u)、Φ^T(u)α、A^TA、A^Tb；

(7)使用共轭梯度法求解式(15)，得到CT图像u

(βΦ^T(u)Φ(u)+λA^TA)u＝βΦ^T(u)α+λA^Tb    (15)

(8)检查是否满足迭代结束的条件，是转至步骤(9)，否则转至步骤(3)；

(9)结束，输出CT图像

本发明由于采用了上述的方法步骤，同现有技术相比，具有如下的有益效果：

1、本方法使用分形作为稀疏表示，使得图像的稀疏表示更加稀疏，进而在求解压缩感知最优化的问题上可以获得精度更高的CT图像。

2、本方法使用ACA算法避免了求解压缩感知最优化过程中数据重排的问题，从而提高的CT图像的成像质量。

附图说明

图1是CT的原理图。

图2是投影矩阵的求解示意图。

具体实施方式

以下进一步说明本发明原理是实现方法：

一、CT原理与压缩感知

CT的原理如图1所示，X光源出射的X光强度为I₀，因为X光的能量很高可以穿透被扫描的物体，则出射的X光强度为I₁，根据lambert-beer定律我们有

I₁＝I₀e^-∫μds    (1)

其中μ为X光路上被扫描物体对经过的X光的衰减所对应的衰减系数，不同的物质具有不同的衰减系数，而CT正是将物体的衰减系数映射到灰度值上便得到了CT图像u。

对(1)式处理可得

$\ln \frac{I_0}{I_1}=\∫\mathrm{\μds}---\left(2\right)$

其中为投影数据，即投影数据为X光光路上的X光衰减系数的线积分，将其离散化我们可得

$\ln \frac{I_0}{I_1}=\mathrm{\Σ}\mathrm{\μ}---\left(3\right)$

将物体的μ映射到灰度值上便得到了CT图像u，由此我们可以建立CT的数学模型，即为线性方程组:

Au＝b(4)

其中b＝(b¹,b²,...b^M)∈R^M为已知的投影数据，将所有的未知像素排列成一列向量即u＝(u₁,u₂,...u_M×M)∈R^N为未知的重建图像，则A＝(a_ij)为已知的投影矩阵即为每个像素对所在的X光路上对X光衰减所占的比重，其求解方法如图2所示，每个像素对应着不同的衰减系数和对该X光束衰减的贡献即投影系数。将像素的权重集中在像素的中心，若X光束经过该像素，则其权重为1即投影矩阵对应的值为1，反之则为0。

CT重建算法分为解析算法和迭代算法，其中解析算法的代表为滤波反投影(Filter Back-projection,FBP)。迭代算法有代数重建技术(Algebraic ReconstructionTechnique,ART)，同时代数重建技术(Simultaneous Algebraic ReconstructionTechnique,SART)，最大似然重建算法(Expectation-maximization,EM)等。其中ART算法的迭代公式如下

${u_j}^{(n+1)}={u_j}^{\left(n\right)}+{\mathrm{\λ}}_n\frac{a_{\mathrm{ij}}}{{\left|\right|A^i\left|\right|}^2}(b^i-A^i{u_j}^{\left(n\right)}),n=\mathrm{0,1},....---\left(5\right)$

在“稀疏”投影的情况下，利用传统的CT重建算法将很难重建出较高质量的图像，尤其是对数据完备性要求很高的解析算法。而对于传统的迭代重建算法，其重建出的图像也存在较多的伪影。而基于压缩感知的CT重建算法可以在少投影的情况下重建出质量高的CT图像。于是，CT重建问题则转化为了有约束的最优化问题

$\underset{u}{\min }{\left|\right|\mathrm{\Φ}\left(u\right)\left|\right|}_1^{1}s.t.\mathrm{Au}=b---\left(6\right)$

其中正则化函数，通常用梯度变换，小波变换等的L1范数表示。为了简化式(2)，利用最优化方法中的惩罚函数法，将其转化为非约束的最优化问题：

$u=\underset{u}{\arg \min }{\left|\right|\mathrm{\Φ}\left(u\right)\left|\right|}_1^1+\mathrm{\λ}{\left|\right|b-\mathrm{Au}\left|\right|}_2^2---\left(7\right)$

其中λ为权重系数。由于传统的梯度变换、小波变换不能很稀疏的表示图像。本文采用分形作为稀疏变换以获得更稀疏的稀疏表示。

二、分型压缩

分形稀疏变换利用图像中不同区域间存在跨尺度的自相似性，把图像建模成分形体来实现图像的稀疏表示的。目前，分形编码以其新颖的思想、高压缩比、分辨率无关性等优点受到学术界广泛的关注。

分形编码的思想为自然图像中包括许多复杂类似于分形的物体，把这些图像设法描述为分形，而不是“光滑”的近似表示它们。使用一种图像近似不变的压缩仿射变换，即迭代函数系统的量化参数来表示图像，这样就大大的减少了数据量，即可以很稀疏的表示图像。

目前常用的分形编码算法为基于分块的迭代函数系统(Iteration FunctionSystem，IFS)的分形块编码算法，即首先将图像分为定义域块和值域块，其中定义域块大小一般为值域块的两倍，再把值域块通过对定义域的匹配获取IFS码，最后存储使用到的定义域块和值域块相应的IFS码就完成了分形编码。

三、本发明的基于分形和压缩感知的CT重建方法

本方法是使用分形作为压缩感知的稀疏变换，即

$u=\underset{u}{\arg \min }{\left|\right|\mathrm{\Φ}\left(u\right)\left|\right|}_1^1+\frac{\mathrm{\λ}}{2}{\left|\right|b-\mathrm{Au}\left|\right|}_2^2---\left(8\right)$

其中Φ(u)为分形稀疏变换。为了解决式(8)，我们采用迭代连续迭代算法进行求解。该方法在求解最小化图像的L1范数拥有独特的优势。为了利用ACA算法，式(8)等价于

$u=\underset{u}{\arg \min }{\left|\right|\mathrm{\α}\left|\right|}_1^1+\frac{\mathrm{\λ}}{2}{\left|\right|b-\mathrm{Au}\left|\right|}_2^{2}s.t.-\mathrm{\α}=\mathrm{\Φ}\left(u\right)---\left(9\right)$

写成无约束形式：

$u=\underset{u}{\arg \min }{\left|\right|\mathrm{\α}\left|\right|}_1^1+\frac{\mathrm{\β}}{2}{\left|\right|\mathrm{\α}-\mathrm{\Φ}\left(u\right)\left|\right|}_2^2+\frac{\mathrm{\λ}}{2}{\left|\right|b-\mathrm{Au}\left|\right|}_2^2---\left(10\right)$

其中β为权重参数。当β固定时，可以通过交替计算得到：

固定u，求解

$\mathrm{\α}=\underset{u}{\arg \min }{\left|\right|\mathrm{\α}\left|\right|}_1^1+\frac{\mathrm{\β}}{2}{\left|\right|\mathrm{\α}-\mathrm{\Φ}\left(u\right)\left|\right|}_2^2---\left(11\right)$

由软阈值方法得到：

$\mathrm{\α}=\max \left(\right|\mathrm{\Φ}\left(u\right)|-\frac{1}{\mathrm{\β}},0)\frac{\mathrm{\Φ}\left(u\right)}{\left|\mathrm{\Φ}\left(u\right)\right|}---\left(12\right)$

固定α_j，求解

$u=\underset{u}{\arg \min }\frac{\mathrm{\β}}{2}{\left|\right|\mathrm{\α}-\mathrm{\Φ}\left(u\right)\left|\right|}_2^2+\frac{\mathrm{\λ}}{2}{\left|\right|b-\mathrm{Au}\left|\right|}_2^2---\left(13\right)$

对u求导得：0＝βΦ^T(u)(α-Φ(u))+λA^T(b-Au)    (14)

整理得：(βΦ^T(u)Φ(u)+λA^TA)u＝βΦ^T(u)α+λA^Tb    (15)

其中式(15)使用共轭梯度法求解便可得到CT图像u。

本方法的算法流程包括如下步骤：

1)已知CT投影数据b和相应的投影角度θ，初始化参数β、λ、λ_n，CT图像u＝0；

2)根据投影角度θ，计算出相应的投影矩阵A；

3)使用式(5)计算CT图像u_ART；

4)对u_ART进行分形编码处理得到Φ(u)，使用式(12)得到α；

5)对Φ(u)做解码处理得到Φ^T(u)

6)分别求解Φ^T(u)Φ(u)、Φ^T(u)α、A^TA、A^Tb；

7)使用共轭梯度法求解式(15)，得到CT图像

8)检查是否满足迭代结束的条件，是转至9)，否则转至3)；

9)结束，输出CT图像

Claims (1)

1.一种基于分形和压缩感知的CT重建方法，所述方法使用分形作为稀疏变换以获得更稀疏的系数表示，使用迭代连续迭代算法进行求解，重建出高质量的CT图像，步骤如下：

(1)已知CT投影数据b和相应的投影角度θ，初始化阈值参数β、压缩感知中的正则化参数λ、ART算法的迭代参数λ_n，初始化CT图像为0即u＝0；

(2)根据投影角度θ，计算出相应的投影矩阵A；

(3)使用式公式(5)计算CT图像u_ART

${u_j}^{(n+1)}={u_j}^{\left(n\right)}+{\mathrm{\λ}}_n\frac{a_{\mathrm{ij}}}{{\left|\right|A^i\left|\right|}^2}(b^i-A^i{u_j}^{\left(n\right)}),n=\mathrm{0,1},....---\left(5\right)$

(4)对u_ART进行分形编码处理得到Φ(u)，使用式(12)得到α

$\mathrm{\α}=\max \left(\right|\mathrm{\Φ}\left(u\right)|-\frac{1}{\mathrm{\β}},0)\frac{\mathrm{\Φ}\left(u\right)}{\left|\mathrm{\Φ}\left(u\right)\right|}---\left(12\right)$

其中Φ(u)为分形稀疏变换，α为ACA算法求解过程中的中间变量，软阈值方法的求解值；

(5)对Φ(u)做解码处理得到Φ^T(u)

(6)分别求解Φ^T(u)Φ(u)、Φ^T(u)α、A^TA、A^Tb；

(7)使用共轭梯度法求解式(15)，得到CT图像u

(βΦ^T(u)Φ(u)+λA^TA)u＝βΦ^T(u)α+λA^Tb   (15)

(8)检查是否满足迭代结束的条件，是转至步骤(9)，否则转至步骤(3)；

(9)结束，输出CT图像