CN102779350B - 一种锥束ct迭代重建算法投影矩阵构建方法 - Google Patents

Abstract

本发明涉及一种锥束CT迭代重建算法，特别是涉及一种锥束CT迭代重建算法投影矩阵构建方法。本发明针对锥束CT迭代重建算法投影矩阵的高精度刻画问题，提出了基于有限元模型和Radon算子的投影矩阵刻画方法。结合射线覆盖模型和基函数模型各自的特点，从一幅连续三维自然图像的数学刻画出发，按照射线投影规律，充分考虑对投影各物理过程的数学刻画，提出一种新的投影矩阵刻画方法，对投影过程进行了更为充分的刻画。实验结果表明，本发明有效提高了模型的刻画精度和重建的质量。

Description

一种锥束CT迭代重建算法投影矩阵构建方法

技术领域

本发明涉及一种锥束CT迭代重建算法，特别是涉及一种锥束CT迭代重建算法投影矩阵构建方法。

背景技术

随着计算机断层成像技术（computerized tomography，CT）在医疗、工业和安防等领域的广泛应用，CT重建算法的研究成为热点问题之一。CT重建算法在数学上可以分为解析重建（Analytic Reconstruction,AR）和迭代重建（Iterative Reconstruction,IR）两类。迭代重建算法相对解析重建算法的优势是能在数据缺失或低信噪比条件下获得更好的成像质量，这对解决医疗、工业和安防中低剂量成像、不完全数据成像等问题具有重要的实际意义。而迭代重建算法的投影矩阵直接联系着已知的投影图像和待求的重建图像，其刻画精度对图像的重建质量起着至关重要的作用。

影响投影矩阵的因素包括焦斑形状、系统几何、探测器响应以及CT系统其他物理参数。目前的研究中，投影矩阵的刻画模型主要分为两类。

一类是对射线束与重建物体体素作用强度的刻画，称为射线覆盖模型，常见的包括点模型、线模型和面积模型，Mueller等指出点模型可能会引起锯齿状伪影，他们的研究中采用线模型进行改进取得了更好的效果。此后，Ziegler等的研究表明面积模型比线模型能更好的抑制伪影。张顺利对点、线、面模型作了实验对比分析。2004年，GE公司CT系统与应用实验室的De Man和Basu提出了一种距离驱动的方法，可归结为线模型的一种解决方案；

另一类是从插值的思想出发形成的基函数模型，常见的插值核或基函数包括高斯函数、双三次样条函数和Kaiser-Bessel函数等。Joseph使用高斯基函数，Lewitt等使用Kaiser-Bessel基函数,Kohler等使用三线性插值基函数，莫会云等使用双三次样条基函数分别进行了投影矩阵基函数模型的讨论。Danielsso等采用射线基函数、体素基函数等的融合进行了研究。2010年飞利浦公司的Ziegler等在专利（US 7672424）“采用与体素相关插值的图像重建方法”中，将放大比、标准化等因素融合进了投影矩阵基函数模型。

射线覆盖模型和基函数模型有着各自的特点。

射线束与重建物体体素立方体作用过程中，由于入射角度、入射点的不同，体素立方体内物质对射线束的衰减也不同，从而引起测量投影值的差异性，射线覆盖模型即是单独刻画这种射线束与重建物体体素作用强度的投影矩阵模型。根据对射线的数学抽象方式不同，射线覆盖模型常见的有点模型、线模型和面积模型等。

投影矩阵的元素w_ij表示第j个重建体素对第i条射线的贡献。射线覆盖模型中w_ij反映的几何关系如图1（为表示方便，以二维为例）所示，射线的宽度为τ，体素各向同性，宽度为Δ。

点模型中，认为射线宽度τ=Δ，体素值集中在体素中心，射线穿过重建物体体素对投影测量值的影响为1，不穿过为零，即，

点模型对投影矩阵的刻画是一种最简的形式，对射线与体素作用过程的刻画也较为粗略，在实际应用中可能会引起锯齿状伪影。但由于其容易计算机实现，是实际中用得较多的一种模型。较点模型刻画得精确一点的模型为线模型，该模型认为射线宽度τ＝0，体素值均匀分布在体素立方体内，投影测量值与射线穿过体素立方体的长度相关，记，

w_ij＝i号射线交j号体素的长度（2）

线模型较点模型更符合成像过程的实际情况，能在一定程度上改善锯齿状伪影。与此类似的还有面模型，该模型认为τ＞0（常取τ=Δ），体素值均匀分布在体素立方体内，探测器检测结果与射线穿过体素立方体的面积相关，记，

文献《ART算法几种重建模型的研究和比较》（航空计算技术，2005(2):39-41，张顺利）对这三个模型进行了研究和比较，并通过仿真实验表明，从重建质量上来看面模型略优于线模型，线模型优于点模型，但是从重建时间上来看，面模型耗时大于线模型，线模型耗时大于点模型。

从射线束与体素立方体的作用强度出发对投影矩阵进行刻画，形成了独立的射线覆盖模型。也有学者根据插值、有限元的思想从体素立方体内物质分布不均匀性的描述出发对投影矩阵进行刻画，形成独立的基函数模型。

基函数模型，根据物理意义对基函数的解释与选择。

迭代重建算法的CT重建模型将重建物体离散化的过程中，将重建物体划分为有限个小的立方体。实际上，每个小的立方体内物质的分布可能是不均匀的。并且，小立方体区域投射到探测器上，可能影响多个投影矩阵的像素。基函数模型即是单独刻画这种体素立方体内物质不均匀性的投影矩阵模型，基函数投影示意图如图4所示。根据基函数构造方法的不同，常见的基函数投影矩阵模型包括Joseph方法、Siddon方法、Koehler方法、高斯基函数模型、双三次样条基函数模型和Kaiser-Bessel函数模型等。Joseph方法，Siddon方法和Koehler方法认为体素立方体是可以彼此相交的，二维示意图如图3a、图3b、图3c所示，基函数区域大小依次为Δ³、和(2Δ)³，

在该区域内当前体素的体素值是均匀分布的，重建物体为各基函数的线性组合。

Joseph方法，Siddon方法和Koehler方法从基函数的区域大小上突破了体素立方体对区域的限制，有利于刻画单个体素对多个投影像素的影响。高斯基函数模型、双三次样条基函数模型和Kaiser-Bessel函数模型等认为在体素立方体内物质的密度分布服从基函数所描述的分布，对定义区域内体素值的不均匀性进行了刻画。

射线覆盖模型和基函数模型分别从不同角度出发单独对投影矩阵进行了刻画，反映了成像过程中不同的几何与物理关系。是否存在一个更为精确的投影矩阵刻画方法，能综合反映成像过程中各几何与物理关系，使投影矩阵数学模型尽可能逼近真实系统，满足实际成像中对投影矩阵的高精度要求，重建出更高质量的重建图像？

尽管CT迭代重建投影矩阵刻画模型的研究已取得以上成果，但随着迭代重建在低剂量，低信噪比以及不完全数据成像等应用领域的优势得到进一步发挥，对其重建质量和模型刻画精度也提出越来越高的要求，目前的模型均为根据几何关系构造的模型或插值模型，其对于影响投影矩阵各因素的刻画较为单一，易引起伪影。

发明内容

为了进一步提高锥束CT迭代重建算法投影矩阵的刻画精度，重建出更高质量的图像，本发明在充分分析常见投影矩阵模型的基础上，对投影矩阵的数学抽象方法做了进一步的探索和研究，提出了一种基于有限元与Radon算子的投影矩阵刻画方法。

本发明所采用的技术方案：

一种锥束CT迭代重建算法投影矩阵构建方法，结合探元灵敏度模型和基函数模型各自的特点，从一幅连续三维自然图像的数学刻画出发，按照射线投影规律，充分考虑对投影各过程的数学刻画，提出一种新的投影矩阵模型，该投影矩阵构建方法包括下述步骤：

（a）三维有限元模型对重建图像的刻画

投影矩阵是反映重建图像和投影图像关系的量，为了得到更为精确的投影矩阵模型，首先研究连续的重建物体的数学刻画。实际上，CT成像描述的是重建物体的密度分布值，而重建物体的密度分布是三维连续的，为了便于重建密度分布值的计算机表达与运算，需要将三维连续的密度分布值转换为离散的数字图像，三维有限元模型提供了一种可能的途径。

有限元方法是对真实系统理想化的数学抽象，属于一种解决实际数学物理问题的数值分析方法，能将一个连续的无限自由度问题转化为离散的有限自由度问题，其基本思想是：将求解区域离散化为一组有限个单元的单元组合体，各单元内可以定义各自的近似函数，通过各单元函数的线性组合逼近求解域上待求的场函数。单元数目越多、近似函数的精度越高，解的近似程度越高。

在重建图像的刻画中引入三维有限元模型，设连续的重建图像（密度分布值）为

其中

为空间点的坐标(x,y,z)，将重建图像区域划分为J＝N×N×N个小立方体，在第j个小立方体内定义局部基函数

实际上，

刻画的是第j个小立方体内物体密度值分布的不均匀性，则重建函数

可表示为这一系列基函数的线性组合，

$f\left(\stackrel{\→}{r}\right)\≈\tilde{f}\left(\stackrel{\→}{r}\right)=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{x}_j{b}_j\left(\stackrel{\→}{r}\right)---\left(1\right)$

其中，是有限元模型对

的近似解，X＝{x_j|j＝1,2,3,…,J}即为重建物体的数字图像，用来近似刻画体素小立方体内物体密度分布的不均匀性。

（b）Radon算子对投影过程的刻画

在CT系统中，X射线从光源S发出，穿过物体O投射到探测器D上，第i条射线穿过第j个体素，投影为p_ij，如图2所示。

由Radon定理可知投影图像为重建图像沿射线的线积分，投影过程可以由Radon算子R_i作用于

进行刻画，

$R_{i}f\left(\stackrel{\→}{r}\right)={\∫s}_i\left(\stackrel{\→}{r}\right)f\left(\stackrel{\→}{r}\right)d\stackrel{\→}{r}---\left(2\right)$

射线与体素立方体作用的过程中，不同入射角度、不同入射点引起射线衰减的程度是不一样的，由其引起的投影值也是不一样的，式中

即是刻画这一现象的参数，称为射线覆盖参数。

（c）投影图像的刻画

Radon算子R_i作用于为第i条射线穿过物体得到的投影值p_i，

$p_i=R_{i}f\left(\stackrel{\→}{r}\right)$

可得，

$p_i=\∫s_i\left(\stackrel{\→}{r}\right)f\left(\stackrel{\→}{r}\right)d\stackrel{\→}{r}\≈{\∫s}_i\left(\stackrel{\→}{r}\right)\left[\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{x}_j{b}_j\left(\stackrel{\→}{r}\right)\right]d\stackrel{\→}{r}---\left(3\right)$

将求和放到积分外，

$p_i\≈\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\left[{\∫s}_i\left(\stackrel{\→}{r}\right)b_j\left(\stackrel{\→}{r}\right)d\stackrel{\→}{r}\right]x_j$

上式即为CT重建离散模型P＝WX的第i个方程，与迭代重建中建立的CT重建问题模型是一致的。

（d）投影矩阵的刻画

考察式上式与CT重建离散模型P＝WX表达的区别，上式中的

即为问题中的投影系数w_ij，定义投影系数，

$w_{\mathrm{ij}}={\∫s}_i\left(\stackrel{\→}{r}\right)b_j\left(\stackrel{\→}{r}\right)d\stackrel{\→}{r}---\left(4\right)$

则投影矩阵，即为本发明的投影矩阵模型。新的投影矩阵模型引入了三维有限元模型和Radon算子，对投影系数进行了更为精细的表达。

对于一个重建规模、投影规模、系统尺寸、投影角度等系统参数确定的系统，投影矩阵的值也是确定的，理论上可以事先计算出来。但由于矩阵维度（(M×M×Θ)×(N×N×N)）巨大，如，M＝256，Θ＝360，N＝256，采用32bit浮点数表示时，其存储空间达到1440TB，带来存储和检索上的困难，因此一般情况下投影矩阵不事先计算和存储，而是在重建中使用时再由式（4）计算投影系数。

新的投影矩阵模型中，

是刻画射线与体素立方体作用程度的量，

是刻画体素立方体内部不均匀性的量。考察以往的模型，射线覆盖模型仅刻画了射线与体素立方体作用强度，可以由本发明模型中

得到；基函数模型仅刻画了体素立方体内部的不均匀性，可以由模型中

得到。可见，本发明投影矩阵模型涵盖了目前常见的模型。

射线覆盖参数

反映射线覆盖体素立方体的程度，本专利以射线覆盖体素立方体的长度即长度模型为例说明，射线覆盖系数的线积分

为第i条射线穿过第j个体素立方体的长度，用表示第i条射线的斜截式方程，则参数

为，

$s_i\left(\stackrel{\→}{r}\right)=\mathrm{\δ}({\stackrel{\→}{k}}_i\·\stackrel{\→}{r}-{\mathrm{\τ}}_i)---\left(5\right)$

其几何关系如图1（以二维为例）所示，射线的宽度为0，体素小方格宽度为Δ。

基函数反映体素立方体内部的不均匀性，本发明以Kaiser-Bessel函数作为基函数

为例进行说明，其在三维空间内是一种球状函数，定义为：

$b_{m,a,\mathrm{\α}}\left(r_{b,j}\right)=\left\{\begin{array}{cc}\frac{I_m\left[\mathrm{\α}\sqrt{1-{(r_{b,j}/a)}^2}\right]}{I_m\left(\mathrm{\α}\right)}{\sqrt{1-{(r_{b,j}/a)}^2}}^m& 0\≤r_{b,j}\≤a\\ 0& \mathrm{else},\end{array}\right.---\left(6\right)$

其中，r_b,j指某点到第j个基函数球心的距离，I_m为m阶的修正Bessel函数，a为球状基函数邻域半径，α为控制函数形状的参数。一组标准参数为：m＝2，a/Δ＝2.00（即基函数的定义域包括两个体素长度），α＝10.4。在时域内一维示意图如图5所示。

确定了模型的一组参数，得到本专利新的投影矩阵模型

的一种实现为，

$w_{\mathrm{ij}}=\∫\mathrm{\δ}({\stackrel{\→}{k}}_i\·\stackrel{\→}{r}-{\mathrm{\τ}}_i)b_{\mathrm{2,2}\mathrm{\Δ},10.4}\left(r_{b,j}\right)d\stackrel{\→}{r}---\left(7\right)$

每次计算投影系数w_ij时，均需计算射线穿过半径为2Δ的球状基函数b_{2，2Δ，10.4}(r_b,j)的连续积分值。设第i条射线与第j个体素中心（即b_{2，2Δ，10.4}(r_b,j)的球心）距离为r_ij，当r_ij值恒定时，射线穿过球体的弦长是恒定的，射线穿过球体的基函数值分布是恒定的，因此，射线穿过球体的连续积分值由r_ij唯一确定，是r_ij的一维函数w_ij＝w(r_ij)。为了便于算法的计算机实现，提高运算效率。本专利对球状的Kaiser-Bessel基函数b_{2，2Δ，10.4}(r_b,j)进行离散化，采用离散线性差值的方法求取投影系数值。

对球状基函数b_{2，2Δ，10.4}(r_b,j)沿半径方向L等分，每一分点l(l＝0,1,2,...,L)距离球心的距离为r_b,j(l)（r_b,j(l)∈[0,2Δ]）。以分点l为中点，过分点l生成球体的任一根弦（得到弦所在直线方程为

），作为当前射线与球体的交线。由（7）式计算每一分点l处的投影系数值，存储为一维的查找表T(l)。重建时仅需计算第j个体素距第i条射线的距离r_ij，调用查找表中的值进行线性插值求取投影系数，令

（

为向下取整），则投影系数值为：

$w_{\mathrm{ij}}=w\left(r_{\mathrm{ij}}\right)=\left\{\begin{array}{cc}T\left(l_{\mathrm{ij}}\right)+\mathrm{\&delta;}*[T(l_{\mathrm{ij}}+1)-T\left(l_{\mathrm{ij}}\right)]& 0\&le;r_{\mathrm{ij}}<2\mathrm{\&Delta;}\\ 0& r_{\mathrm{ij}}\&GreaterEqual;2\mathrm{\&Delta;}\end{array}\right.---\left(8\right)$

本发明的有益效果

本发明针对锥束CT迭代重建算法投影矩阵的高精度刻画问题，提出了基于有限元模型和Radon算子的投影矩阵刻画方法。结合射线覆盖模型和基函数模型各自的特点，从一幅连续三维自然图像的数学刻画出发，按照射线投影规律，充分考虑对投影各物理过程的数学刻画，提出一种新的投影矩阵刻画方法，对投影过程进行了更为充分的刻画。实验结果表明，本发明有效提高了模型的刻画精度和重建的质量。

本发明锥束CT迭代重建算法投影矩阵构建方法，通过有限元模型对体素立方体内物质的不均匀性进行了估计，通过Radon算子对射线与体素立方体作用强度进行了估计，证明了目前常见的模型是本发明投影矩阵模型在特定参数下的特例，并采用离散线性差值的方法对新模型的实现进行了优化。

附图说明

图1：射线覆盖模型中射线与体素关系；

图2：锥束CT投影示意图；

图3a：插值基函数示意图之Siddon方法；

图3b：插值基函数示意图之Joseph方法；

图3c：插值基函数示意图之Koehler方法；

图4：基函数的投影；

图5：Kaiser-Bessel基函数一维示意图；

图6a：De Man方法60°时投影图像精度比较；

图6b：Ziegler方法60°时投影图像精度比较；

图6c：本专利方法60°时投影图像精度比较；

图7a～图7d：不同投影矩阵下重建图像切片图比较，其中图7a为Man方法，图7b为Ziegler方法，图7c为本专利方法，图7d为Shepp-Logan头模；

图8：不同投影矩阵下重建图像剖线（Z=64，Y=64）灰度图。

具体实施方式

实施例一：

本发明“锥束CT迭代重建算法投影矩阵构建方法”，结合射线覆盖模型和基函数模型各自的特点，从一幅连续三维自然图像的数学刻画出发，按照射线投影规律，充分考虑对投影过程的数学刻画，提出一种新的投影矩阵模型，该投影矩阵构建方法包括下述步骤：

1）三维有限元模型对重建图像的刻画

在重建图像的刻画中引入三维有限元模型，设连续的重建图像（密度分布值）为其中

为空间点的坐标(x,y,z)，将重建图像区域划分为J＝N×N×N个小立方体，在第j个小立方体内定义局部基函数

实际上，

刻画的是第j个小立方体内物体密度值分布的不均匀性，则重建函数

可表示为这一系列基函数的线性组合，

$f\left(\stackrel{\&RightArrow;}{r}\right)\&ap;\tilde{f}\left(\stackrel{\&RightArrow;}{r}\right)=\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{x}_j{b}_j\left(\stackrel{\&RightArrow;}{r}\right)---\left(1\right)$

其中，

是有限元模型对

的近似解，X＝{x_j|j＝1,2,3,…,J}即为重建物体的数字图像，

用来近似刻画体素小立方体内物体密度分布的不均匀性。

2）Radon算子对投影过程的刻画；

由Radon定理可知投影图像为重建图像沿射线的线积分，投影过程可以由Radon算子R_i作用于

进行刻画，

$R_{i}f\left(\stackrel{\&RightArrow;}{r}\right)={\&Integral;s}_i\left(\stackrel{\&RightArrow;}{r}\right)f\left(\stackrel{\&RightArrow;}{r}\right)d\stackrel{\&RightArrow;}{r}---\left(2\right)$

射线与体素立方体作用的过程中，不同入射角度、不同入射点引起射线衰减的程度是不一样的，由其引起的投影值也是不一样的，式中

即是刻画这一现象的参数，称为射线覆盖参数。

3）投影图像的刻画：

Radon算子R_i作用于

为第i条射线穿过物体得到的投影值p_i，即有，

$p_i=R_{i}f\left(\stackrel{\&RightArrow;}{r}\right)=\&Integral;s_i\left(\stackrel{\&RightArrow;}{r}\right)f\left(\stackrel{\&RightArrow;}{r}\right)d\stackrel{\&RightArrow;}{r}\&ap;{\&Integral;s}_i\left(\stackrel{\&RightArrow;}{r}\right)\left[\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{x}_j{b}_j\left(\stackrel{\&RightArrow;}{r}\right)\right]d\stackrel{\&RightArrow;}{r}---\left(3\right)$

4）投影矩阵的刻画：

由上式，定义投影系数，

$w_{\mathrm{ij}}={\&Integral;s}_i\left(\stackrel{\&RightArrow;}{r}\right)b_j\left(\stackrel{\&RightArrow;}{r}\right)d\stackrel{\&RightArrow;}{r}---\left(4\right)$

则投影矩阵W＝{w_ij}，即为本专利新的投影矩阵模型，其中，

是刻画射线与体素立方体作用程度的量，

是刻画体素立方体内部不均匀性的量。

实施例二：

本实施例的锥束CT迭代重建算法投影矩阵构建方法，与实施例一不同的是：射线覆盖参数

反映射线覆盖体素立方体的程度，可抽象成点、线或面。如采用射线覆盖体素立方体的长度即长度模型，射线覆盖系数的线积分

为第i条射线穿过第j个体素立方体的长度，用

表示第i条射线的斜截式方程，则参数

为，

$s_i\left(\stackrel{\&RightArrow;}{r}\right)=\mathrm{\&delta;}({\stackrel{\&RightArrow;}{k}}_i\&CenterDot;\stackrel{\&RightArrow;}{r}-{\mathrm{\&tau;}}_i)---\left(5\right)$

实施例三：

本实施例的锥束CT迭代重建算法投影矩阵构建方法，与实施例一或实施例二不同的是：基函数反映体素立方体内部的不均匀性，可采用的函数包括Joseph方法、Siddon方法、Koehler方法、高斯基函数、双三次样条基函数和Kaiser-Bessel函数等。如Kaiser-Bessel函数作为基函数时，其在三维空间内是一种球状函数，定义为：

$b_{m,a,\mathrm{\&alpha;}}\left(r_{b,j}\right)=\left\{\begin{array}{cc}\frac{I_m\left[\mathrm{\&alpha;}\sqrt{1-{(r_{b,j}/a)}^2}\right]}{I_m\left(\mathrm{\&alpha;}\right)}{\sqrt{1-{(r_{b,j}/a)}^2}}^m& 0\&le;r_{b,j}\&le;a\\ 0& \mathrm{else},\end{array}\right.---\left(6\right)$

其中，r_b,j指某点到第j个基函数球心的距离，I_m为m阶的修正Bessel函数，a为球状基函数邻域半径，α为控制函数形状的参数。一组标准参数为：m＝2，a/Δ＝2.00（即基函数的定义域包括两个体素长度），α＝10.4。

实施例四：

本实施例锥束CT迭代重建算法投影矩阵构建方法，包括下述步骤：

（a）三维有限元模型对重建图像的刻画

投影矩阵是反映重建图像和投影图像关系的量，为了得到更为精确的投影矩阵模型，首先研究连续的重建物体的数学刻画。实际上，CT成像描述的是重建物体的密度分布值，而重建物体的密度分布是三维连续的，为了便于重建密度分布值的计算机表达与运算，需要将三维连续的密度分布值转换为离散的数字图像，三维有限元模型提供了一种可能的途径。

有限元方法是对真实系统理想化的数学抽象，属于一种解决实际数学物理问题的数值分析方法，能将一个连续的无限自由度问题转化为离散的有限自由度问题，其基本思想是：将求解区域离散化为一组有限个单元的单元组合体，各单元内可以定义各自的近似函数，通过各单元函数的线性组合逼近求解域上待求的场函数。单元数目越多、近似函数的精度越高，解的近似程度越高。

在重建图像的刻画中引入三维有限元模型，设连续的重建图像（密度分布值）为

其中

为空间点的坐标(x,y,z)，将重建图像区域划分为J＝N×N×N个小立方体，在第j个小立方体内定义局部基函数实际上，

刻画的是第j个小立方体内物体密度值分布的不均匀性，则重建函数可表示为这一系列基函数的线性组合，

$f\left(\stackrel{\&RightArrow;}{r}\right)\&ap;\tilde{f}\left(\stackrel{\&RightArrow;}{r}\right)=\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{x}_j{b}_j\left(\stackrel{\&RightArrow;}{r}\right)---\left(7\right)$

其中，

是有限元模型对

的近似解，X＝{x_j|j＝1,2,3,…,J}即为重建物体的数字图像，

用来近似刻画体素小立方体内物体密度分布的不均匀性。

（b）Radon算子对投影过程的刻画

在CT系统中，X射线从光源S发出，穿过物体O投射到探测器D上，第i条射线穿过第j个体素，投影为p_ij，如图2所示。

由Radon定理可知投影图像为重建图像沿射线的线积分，投影过程可以由Radon算子R_i作用于

进行刻画，

$R_{i}f\left(\stackrel{\&RightArrow;}{r}\right)={\&Integral;s}_i\left(\stackrel{\&RightArrow;}{r}\right)f\left(\stackrel{\&RightArrow;}{r}\right)d\stackrel{\&RightArrow;}{r}---\left(8\right)$

射线与体素立方体作用的过程中，不同入射角度、不同入射点引起射线衰减的程度是不一样的，由其引起的投影值也是不一样的，式中

即是刻画这一现象的参数，称为射线覆盖参数。

（c）投影图像的刻画

Radon算子R_i作用于

为第i条射线穿过物体得到的投影值p_i，

$p_i=R_{i}f\left(\stackrel{\&RightArrow;}{r}\right)---\left(9\right)$

可得，

$p_i=\&Integral;s_i\left(\stackrel{\&RightArrow;}{r}\right)f\left(\stackrel{\&RightArrow;}{r}\right)d\stackrel{\&RightArrow;}{r}\&ap;{\&Integral;s}_i\left(\stackrel{\&RightArrow;}{r}\right)\left[\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{x}_j{b}_j\left(\stackrel{\&RightArrow;}{r}\right)\right]d\stackrel{\&RightArrow;}{r}---\left(10\right)$

将求和放到积分外，

$p_i\&ap;\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}\left[{\&Integral;s}_i\left(\stackrel{\&RightArrow;}{r}\right)b_j\left(\stackrel{\&RightArrow;}{r}\right)d\stackrel{\&RightArrow;}{r}\right]x_j---\left(11\right)$

上式即为CT重建离散模型P＝WX的第i个方程，与迭代重建中建立的CT重建问题模型是一致的。

（d）投影矩阵的刻画

考察式上式与CT重建离散模型P＝WX表达的区别，上式中的

即为问题中的投影系数w_ij，定义投影系数，

$w_{\mathrm{ij}}={\&Integral;s}_i\left(\stackrel{\&RightArrow;}{r}\right)b_j\left(\stackrel{\&RightArrow;}{r}\right)d\stackrel{\&RightArrow;}{r}---\left(12\right)$

则投影矩阵W＝{w_ij}，即为本文新的投影矩阵模型。新的投影矩阵模型引入了三维有限元模型和Radon算子，对投影系数进行了更为精细的表达。

对于一个重建规模、投影规模、系统尺寸、投影角度等系统参数确定的系统，投影矩阵的值也是确定的，理论上可以事先计算出来。但由于矩阵维度（(M×M×Θ)×(N×N×N)）巨大，如，M＝256，Θ＝360，N＝256，采用32bit浮点数表示时，其存储空间达到1440TB，带来存储和检索上的困难，因此一般情况下投影矩阵不事先计算和存储，而是在重建中使用时再由式（12）计算投影系数。

新的投影矩阵模型中，

是刻画射线与体素立方体作用程度的量，

是刻画体素立方体内部不均匀性的量。考察以往的模型，射线覆盖模型仅刻画了射线与体素立方体作用强度，可以由本发明矩阵模型中

得到；基函数模型仅刻画了体素立方体内部的不均匀性，可以由本发明模型中

得到。可见，本发明矩阵模型是投影矩阵的更为一般的模型，涵盖了目前常见的模型。

射线覆盖参数

反映射线覆盖体素立方体的程度，本专利以射线覆盖体素立方体的长度即长度模型为例说明，射线覆盖系数的线积分为第i条射线穿过第j个体素立方体的长度，用

表示第i条射线的斜截式方程，则参数

为，

$s_i\left(\stackrel{\&RightArrow;}{r}\right)=\mathrm{\&delta;}({\stackrel{\&RightArrow;}{k}}_i\&CenterDot;\stackrel{\&RightArrow;}{r}-{\mathrm{\&tau;}}_i)---\left(13\right)$

其几何关系如图1（以二维为例）所示，射线的宽度为0，体素小方格宽度为Δ。

基函数反映体素立方体内部的不均匀性，本专利以Kaiser-Bessel函数作为基函数

为例进行说明，其在三维空间内是一种球状函数，定义为：

$b_{m,a,\mathrm{\&alpha;}}\left(r_{b,j}\right)=\left\{\begin{array}{cc}\frac{I_m\left[\mathrm{\&alpha;}\sqrt{1-{(r_{b,j}/a)}^2}\right]}{I_m\left(\mathrm{\&alpha;}\right)}{\sqrt{1-{(r_{b,j}/a)}^2}}^m& 0\&le;r_{b,j}\&le;a\\ 0& \mathrm{else},\end{array}\right.---\left(14\right)$

其中，r_b,j指某点到第j个基函数球心的距离，I_m为m阶的修正Bessel函数，a为球状基函数邻域半径，α为控制函数形状的参数。一组标准参数为：m＝2，a/Δ＝2.00（即基函数的定义域包括两个体素长度），α＝10.4。在时域内一维示意图如图5所示。

确定了模型的一组参数，得到本专利新的投影矩阵模型

的一种实现为，

$w_{\mathrm{ij}}=\&Integral;\mathrm{\&delta;}({\stackrel{\&RightArrow;}{k}}_i\&CenterDot;\stackrel{\&RightArrow;}{r}-{\mathrm{\&tau;}}_i)b_{\mathrm{2,2}\mathrm{\&Delta;},10.4}\left(r_{b,j}\right)d\stackrel{\&RightArrow;}{r}---\left(15\right)$

每次计算投影系数w_ij时，均需计算积分式（15，即计算射线

穿过半径为2Δ的球状基函数b_{2，2Δ，10.4}(r_b,j)的连续积分值。设第i条射线与第j个体素中心（即b_2，2Δ,10.4(r_b,j)的球心）距离为r_ij，当r_ij值恒定时，射线穿过球体的弦长是恒定的，射线穿过球体的基函数值分布是恒定的，因此，射线穿过球体的连续积分值由r_ij唯一确定，是r_ij的一维函数w_ij＝w(r_ij)。为了便于算法的计算机实现，提高运算效率。

本实施例通过对球状的Kaiser-Bessel基函数b_{2，2Δ，10.4}(r_b,j)进行离散化，采用离散线性差值的方法求取投影系数值。

对球状基函数b_{2，2Δ，10.4}(r_b,j)沿半径方向L等分，每一分点l(l＝0,1,2,...,L)距离球心的距离为r_b,j(l)（r_b,j(l)∈[0,2Δ]）。以分点l为中点，过分点l生成球体的任一根弦（得到弦所在直线方程为），作为当前射线与球体的交线。由（15）式计算每一分点l处的投影系数值，存储为一维的查找表T(l)。重建时仅需计算第j个体素距射第i条射线的距离r_ij，调用查找表中的值进行线性插值求取投影系数，令

（

为向下取整），则投影系数值为：

$w_{\mathrm{ij}}=w\left(r_{\mathrm{ij}}\right)=\left\{\begin{array}{cc}T\left(l_{\mathrm{ij}}\right)+\mathrm{\&delta;}*[T(l_{\mathrm{ij}}+1)-T\left(l_{\mathrm{ij}}\right)]& 0\&le;r_{\mathrm{ij}}<2\mathrm{\&Delta;}\\ 0& r_{\mathrm{ij}}\&GreaterEqual;2\mathrm{\&Delta;}\end{array}\right.---\left(16\right)$

实验与结果

采用Shepp-Logan头模通过不同投影矩阵模型下的投影实验和迭代重建实验讨论和验证模型刻画精度和重建质量。实验平台CPU为Pentium(R)4，主频3.00GHz，内存3GB，GPU为GeForce G100，开发平台为自研的MatrixCloud（sourceforge.net/projects/matrixcloud/）。实验参数重建规模为128×128×128，投影规模为128×128，包含圆轨迹均匀分布的36个投影角度。

实验设置两组对比组，第一组为2004年GE公司Man等的方法，该方法属于线模型的一种解决方案，可由本发明模型

实现。第二组为2010年飞利浦公司Ziegler等在专利中提到的方法，该方法属于基函数模型的一种解决方案，可由本发明模型中

实现。本发明模型参数采用

$b_j\left(\stackrel{\&RightArrow;}{r}\right)=b_{m,a,\mathrm{\&alpha;}}\left(r_{b,j}\right).$

在投影实验中，对标准Shepp-Logan头模在不同的投影矩阵模型下进行投影，取60°时投影图像精度进行比较，实验结果如图6所示。

取60°时投影图像九宫格的中间部分为感兴趣区，计算感兴趣区信噪比（均值比方差）如表1所示。信噪比越大，说明图像质量越好，从而投影矩阵模型对实际过程的刻画越精确。

表1投影实验信噪比比较

从三种方法投影图像可以看出，Man方法存在一定线状伪影。这种伪影在Ziegler方法中得到一定程度改善，但依然较为明显。本专利的模型从很大程度上抑制了这种伪影。通过感兴趣区信噪比的结果，也可看出本专利模型在改善投影精度方面的优势。

在迭代重建实验中，重建算法采用SART算法。如图7a～图7d所示，分别得到不同投影矩阵下重建图像Z=64的切片图。为便于比较，取四幅重建图像Z=64，Y=64的剖线图如图8所示。

为了客观的比较不同投影矩阵的重建图像效果，分别用归一化均方根距离d、归一化平均绝对距离r和最坏情况距离e度量图像质量。这三种距离越小，说明重建图像和原始头像越接近，在重建算法一致的条件下，投影矩阵模型刻画越精确。

表2重建质量比较

通过重建质量的对比分析，特别是从局部放大的剖线灰度图可见，本专利方法最接近于理想值，重建结果优于线模型和基函数模型。本专利方法由于综合考虑了射线与体素立方体作用强度和体素立方体的不均匀性，对于投影矩阵的刻画更符合实际作用过程，得到质量更优的投影和重建图像。

Claims (3)

1.一种锥束CT迭代重建算法投影矩阵构建方法，结合射线覆盖模型和基函数模型各自的特点，从一幅连续三维自然图像的数学刻画出发，按照射线投影规律，充分考虑对投影过程的数学刻画，提出一种新的投影矩阵模型，其特征是：该投影矩阵构建方法包括下述步骤： 

1）三维有限元模型对重建图像的刻画： 

在重建图像的刻画中引入三维有限元模型，设连续的重建图像为

其中

为空间点的坐标(x,y,z)，将重建图像区域划分为J=N′N′N个小立方体，在第j个小立方体内定义局部基函数实际上，

刻画的是第j个小立方体内物体密度值分布的不均匀性，则重建函数可表示为这一系列基函数的线性组合， 

其中，

是有限元模型对

的近似解，X={x_j|j=1,2,3,...,J}即为重建物体的数字图像，用来近似刻画体素小立方体内物体密度分布的不均匀性； 

基函数反应体素立方体内部的不均匀性，可采用的函数包括Joseph方法、Siddon方法、Koehler方法、高斯基函数、双三次样条基函数和Kaiser-Bessel函数；采用Kaiser-Bessel函数作为基函数

时，其在三维空间内是一种球状函数，定义为： 

其中，r_b,j指某点到第j个基函数球心的距离，I_m为m阶的修正Bessel函数，a为球状基函数邻域半径，α为控制函数形状的参数；一组标准参数为：m=2，a/Δ=2.00，即基函数的定义域包括两个体素长度，α=10.4； 

2）Radon算子对投影过程的刻画： 

由Radon定理可知投影图像为重建图像沿射线的线积分，投影过程可以由Radon算子R_i作用于

进行刻画， 

射线与体素立方体作用的过程中，不同入射角度、不同入射点引起射线衰减的程度是不一样的，由其引起的投影值也是不一样的，式中

即是刻画这一现象的参数，称为射线覆盖参数； 

3）投影图像的刻画： 

Radon算子R_i作用于

为第i条射线穿过物体得到的投影值p_i，即有 

4）投影矩阵的刻画： 

由上式，定义投影系数， 

则投影矩阵W={w_ij}，即为投影矩阵模型，其中，

是刻画射线与体素立方体作用程度的量，

是刻画体素立方体内部不均匀性的量。 

2.根据权利要求1所述的锥束CT迭代重建算法投影矩阵构建方法，其特征是：射线覆盖参数

反应射线覆盖体素立方体的程度，可抽象成点、线或面；对于采用射线覆盖体素立方体的长度即长度模型，射线覆盖系数的线积分

为第i条射线穿过第j个体素立方体的长度，用表示第i条射线的斜截式方程，则参数为， 

3.根据权利要求1所述的锥束CT迭代重建算法投影矩阵构建方法，其特征是：对球状的Kaiser-Bessel基函数b_2,2Δ,10.4(r_b,j)进行离散化，采用离散线性差值的方法求取投影系数值： 

对球状基函数b_2,2Δ,10.4(r_b,j)沿半径方向L等分，每一分点l(l=0,1,2,...,L)距离球心的距离为r_b,j(l)（r_b,j(l)∈[0,2Δ]），以分点l为中点，过分点l生成球体的任一根弦，得到弦所在直线方程为

作为当前射线与球体的交线；计算每一分点l处的投影系数值，存储为一维的查找表T(l)；重建时计算第j个体素距射第i条射线的距离r_ij，调用查找表中的值进行线性插值求取投影系数，令

（

为向下取整），则投影系数值为： 

Images