CN104899827A - 基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法 - Google Patents

Abstract

基于固定分辨率条件下的离散Radon和Mojette投影转换方法，属于图像处理领域的计算机成像技术领域。其特征是将实际成像条件下的稀疏角度离散Radon投影转化为在离散域精确重建的Mojette投影。以计算机层析成像算法为基础，通过构造合理的成像系统来逼近Radon和Mojette投影场景，并在详细分析了Mojette投影和Radon投影的关系之后，给出了在固定分辨率前提下各投影角度下的Radon投影转化为对应的投影矢量下的Mojette投影的具体算法。本发明的效果和益处是在不同的稀疏投影角度下，不需要改变Radon投影图像的分辨率大小，克服了Mojette投影域分辨率随着投影矢量变化和Radon投影分辨率不变之间的转换障碍，搭建起这两种投影之间的转换桥梁。

Description

基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法

技术领域

本发明属于图像处理领域里的计算机层析成像(Computed Tomography，CT)技术领域，涉及在图像投影变换域的稀疏采样及高效复原问题。

背景技术

CT技术是一种基于计算机的三维重构技术，能够在不损坏物体表层和内部构造的前提下，重现出物体不可见的内部切片结构。在CT成像算法上，按照成像原理的不同，分为解析类和迭代类重建算法；按照重建精确性的不同，又分为近似重建和精确重建算法。在解析算法中，奥地利数学学家Radon提出的变换模型利用低维的压缩投影数据来重构高维度的物体结构，这一数学模型的提出对CT的图像重建技术起到了关键性的启蒙和指导作用。在这一理论基础的指引下，滤波反投影等解析算法相继被提出，该算法在完备的投影角度下，能精确的重建出图像断层，且能抑制在重建过程中产生的各种伪影，在CT成像领域具有里程碑式的意义。

但由于成像系统是数字离散系统，因而连续的投影模型无法适用于离散的图像，故而离散Radon变换随机应运而生。离散Radon变换及其逆变换解决了从模拟域向数字域转换的问题，但基于Radon逆变换的重建算法需要大量的采样样本和投影数量，才能重建出较好的重建断层，这种重建需求必然会导致过高的辐射剂量和过长的重建时间，出于对医学成像中病人身体健康和病灶动态变化会导致伪影的考虑，在短时间、低剂量和稀疏角度下进行高质量成像就显得尤为重要，同时又充满挑战。

在Radon变换的基础上，Guedon J P等提出了Mojette变换的概念，该变换 是离散Radon变换的一种特殊形式，Mojette变换通过改变不同投影矢量下的采样率，能够在最大程度上避免像素点重复和冗余采样，以避免出现一部分像素点过采样，而另一部分像素点欠采样而带来的重建效果上的不精确。并通过充分利用已重建出来的像素点的信息，从而大大减小了重建断层所需的投影角度和投影射线条数。

但由于Mojette投影系统要求探测器分辨率随着投影角度的变化而变化，且在投影角度之间的步进角也不是一个固定的常值，这与现行的采样模式大相径庭，实际的投影系统符合Radon变换的采样模式，探测器分辨率不随投影角度而变化，因此通过Radon变换得到的投影值不可能直接参与到Mojette逆变换当中来，而是先要将Radon投影转换到相关的Mojette投影上去，再利用Mojette逆变换来进行反重建。

完成从Radon投影到Mojette投影的转换，其中需要解决几个实际的问题：1.投影矢量和投影角度之间的匹配转换，由于Mojette投影系统中采用一对互质的整数来表达投影方向，称之为投影矢量，而Radon投影系统中采用的是极坐标中的旋转角度，故而需要通过简单的换算将Mojette投影矢量表达的方向转换成Radon投影中的旋转角；2.Mojette变换中可变分辨率与Radon变换中固定分辨率之间的矛盾，由于Mojette变换中的探测器分辨率或者说是投影射线条数为随着投影矢量而变化的函数，而Radon投影的探测器分辨率是不随着投影角度而变化的，这就使得采样系统中得到的Radon投影需要通过本发明中提出的算法进行转换，才能得到可以高效重建图像的Mojette投影。

发明内容

本发明的目的是提供了一种基于固定分辨率采集的Radon投影转化为 Mojette投影的算法，该方法在详细分析了Mojette投影和Radon投影的关系之后，解决了固定分辨率条件下的两投影之间的转化难题，给出了Radon投影转化为对应投影矢量下的Mojette投影的充要条件，并基于此给出了具体的转化算法，使得实际成像条件下采集的Radon投影能够转化为Mojette域的投影。

本发明的技术方案的原理及步骤如下：

一、技术方案的原理是：

1.建立投影空间坐标系：设投影空间坐标系用x-o-y表示，如图1所示，坐标原点o为待重建物体的几何中心，x轴正向与探测器像元索引号增加方向一致，y轴正向与射线传播方向一致，x-o-y平面内的一点的空间坐标为(x,y)。光源与探测器在x-o-y平面内逆时针旋转，旋转角度为θ，设旋转后的空间坐标系用x_r-o-y_r表示，像素点的旋转坐标为(x_r,y_r)。在连续域的二维图像空间内，将物体所在区域离散化为M×N的离散小块，设ObjSize为待重建矩形区域的行方向上的边长，则每一个重建像素的物理尺寸大小为ObjPixel＝ObjSize/M，如图2所示。其中，(i,j)为标记离散小块行列位置的索引坐标，即若利用二维矩阵存储该离散图像，则i代表在矩阵中的行数，j代表在矩阵中的列数，通常将图像左上角第一个小块的位置记为索引坐标起始位置(1,1)。

明确重建断层像素在x-o-y平面内的空间坐标(x,y)和在旋转坐标系x_r-o-y_r平面内旋转坐标(x_r,y_r)之间的关系，如式(1)所示。

$\left[\begin{array}{c}{\text{x}}_r\\ y_r\end{array}\right]=\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]---\left(1\right)$

明确重建断层像素在x-o-y平面内的空间坐标(x,y)和离散域索引坐标(i,j)之间的关系，通常，以待重建断层的左上角像素为起点，记该像素的索引坐标为(1,1)，索引坐标为(1,1)的像素中心点的空间坐标(x,y)为：

((1-(M+1)/2)·ObjPixel,(1-(N+1)/2)·ObjPixel)

索引坐标为(i,j)的像素中心点的空间坐标(x,y)为：

((i-(M+1)/2)·ObjPixel,(j-(N+1)/2)·ObjPixel)

基于以上两种坐标的对应，则得到索引坐标为(i,j)的像素中心点在旋转坐标系中的旋转坐标(x_r,y_r)为：

x_r＝[(j-(N+1)/2)·ObjPixel·sin(θ)+(i-(M+1)/2)·ObjPixel·cos(θ)]

y_r＝[-(i-(M+1)/2)·ObjPixel·sin(θ)+(j-(N+1)/2)·ObjPixel·cos(θ)]

2.离散Radon变换：对分辨率大小为M×N的离散图像f进行离散Radon变换，其过程用式(2)来表述：

${\mathrm{RFT}}_{\mathrm{\θ}}\left(x_r\right)=\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}f(i,j){\mathrm{\α}}_{i,j,\mathrm{\θ},x_r}---\left(2\right)$

其中，RFT_θ表示在投影角度θ下的Radon投影，RFT_θ(x_r)表示在投影角度θ下打在探测器x_r位置处的Radon投影值。f(i,j)代表待重建的图像切片上索引坐标为(i,j)的一点的灰度值。若以(θ,x_r)标记投影角度θ下打在探测器上的物理偏移量为x_r的一条投影射线，表示在投影射线(θ,x_r)和象素点(i,j)之间的权重核函数，通常为δ函数或0次样条函数(即将投影贡献权值等比例的视为穿过该象素点的线段长值)。

设在所有投影角度θ下穿过像素中心点的射线的线段权值为1，而对穿过该像素其它位置的的射线，则需要根据这些射线与中心射线之间的垂直距离|Δx_r|，计算出其线段权重值其具体的计算公式如式(3)所示。

$\left\{\begin{array}{c}\mathrm{if}\left|\mathrm{\&Delta;}{x}_r\right|\&le;\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|-|\sin \mathrm{\&theta;}\left|\right)}{2},\\ {\mathrm{\&alpha;}}_{i,j,\mathrm{\&theta;},x_r}=1\\ \mathrm{if}\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|-|\sin \mathrm{\&theta;}\left|\right)}{2}<\left|\mathrm{\&Delta;}{x}_r\right|\&le;\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)}{2}\\ {\mathrm{\&alpha;}}_{i,j,\mathrm{\&theta;},x_r}=\frac{-2\&CenterDot;\left|\mathrm{\&Delta;}{x}_r\right|+\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)\&CenterDot;\mathrm{ObjPixel}}{2\&CenterDot;\mathrm{ObjPixel}\&CenterDot;\min \left(\right|\cos \mathrm{\&theta;},\sin \mathrm{\&theta;}\left|\right)}\\ \mathrm{if}\left|\mathrm{\&Delta;}{x}_r\right|>\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)}{2}\\ {\mathrm{\&alpha;}}_{i,j,\mathrm{\&theta;},x_r}=0\end{array}\right.---\left(3\right)$

其中，ObjPixel代表每一个重建像素的物理尺寸。

3.Mojette变换：Mojette变换中用一对互质的整数(p,q)来表达投影方向，一般有p∈Z,q∈Z⁺，p代表了图像列方向上的整数位移，q代表了图像行方向上的整数位移，投影矢量(p,q)表达的投影角度θ＝tan^-1(q/p)，对分辨率大小为M×N的离散图像f进行Mojette变换，其过程用式(4)来表述：

${\mathrm{Moj}}_{p,q}\left(\mathrm{bin}\right)=\underset{i=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}f(i,j)\mathrm{\&delta;}(\mathrm{bin}-p(i-1)-q(j-1)+\mathrm{po})---\left(4\right)$

其中，Moj_p,q表示在投影矢量(p,q)下的Mojette投影，Moj_p,q(bin)表示在投影矢量(p,q)下打在探测器bin上Mojette投影值。f(i,j)代表待重建的图像切片上索引坐标为(i,j)的一点的灰度值，po为一个探测器像元位置矫正量，当投影矢量p>0时，po＝0；当投影矢量p<0时，po＝(N-1)·p。

Mojette变换与Radon变换相比最大的不同在于，Radon变换中探测器分辨率为固定的常值，穿过物体切片的Radon投影射线之间距离是固定的，因而不能保证每条Radon投影射线都能穿过像素的中心点，穿过像素中心点的线段权值记为1，而不穿过像素中心点的线段权值要根据具体线段长度而定，计算方法如式(3)所示；而Mojette变换中每条投影射线都只穿过像素的中心点，因此不同投影矢量(p,q)下探测器分辨率B(M,N,p,q)不同，探测器分辨率的值是由 图像重建分辨率M×N和投影矢量(p,q)共同决定的，即投影矢量(p,q)下的Mojette探测器分辨率为B(P,Q,p,q)＝(Q-1)|p|+(P-1)|q|+1。这意味着在覆盖相同直径范围内的切片时，不同投影角度下的Mojette投影射线之间的间距h_i不同：

$h_i=\frac{1}{\sqrt{p^2+q^2}}---\left(5\right)$

此外，Radon变换在模拟域重建是精确的，在离散域重建是近似的。而Mojette变换在离散域的重建是精确的。

4.Mojette投影矩阵的构造方法：

基于Mojette正变换原理，下面阐述一下Mojette投影矩阵的构造方法。

1)初始化参数,设重建断层分辨率为M×N，投影矢量为(p,q)，投影矩阵为 其中，B＝(N-1)|p_i|+(M-1)|q_i|+1，V＝M·N。待重建图像列矢量化后的一维向量为x，投影矢量(p,q)下的Mojette投影为Moj_p,q，则有System_Moj_p,q·x＝Moj_p,q；

2)在每个投影矢量(p,q)下，遍历每一个像素点(i,j)，计算其打在探测器上的位置bin：

$\mathrm{po}=\left\{\begin{array}{ccc}0& \mathrm{if}& p>0\\ (N-1)\&CenterDot;p& \mathrm{if}& p<0\end{array}\right.---\left(6\right)$

bin＝p·(i-1)+q·(j-1)-po+1

3)将通过投影矢量(p_i,q_i)建立起来的断层上的一点(i,j)与探测器像元bin之间的投影映射关系，存储到线性投影矩阵System_Moj_p,q中，即：

System_Moj_p,q(bin,(i-1)·N+j)＝1；

4)当遍历完所有像素点后，投影矩阵System_Moj_p,q构造完成。

投影矩阵System_Moj_p,q中每一行代表着一条Mojette投影射线，每一列代表着穿过一个像素点的所有投影射线的线段权值。投影方程中的一点System_Moj_p,q(bin,(i-1)·N+j)意味着列矢量化后的第(i-1)·N+j个像素在第bin条投影射线中的贡献程度，在Mojette投影模型中，若第bin条投影射线 穿过第(i-1)·N+j个像素的中心点，则System_Moj_p,q(bin,(i-1)·N+j)＝1，否则为0。

5.Radon投影矩阵的构造方法：

在Radon投影当中，线性投影矩阵的构造方法是由探测器和射线源与待重建物体之间的几何关系决定的。

1)设重建断层分辨率为M×N，包含待重建断层的最小正方形区域的边长为ObjSize，每一个重建像素的实际物理尺寸ObjPixel＝ObjSize/M。设探测器分辨率为DetRowNum，其中每一个探测器像元的物理尺寸为DetCCDSize，投影矩阵为System_rad_θ∈R^J×V，J＝DetRowNum，V＝M·N。待重建图像列矢量化后的一维向量为x，投影角度θ下的Radon投影为RFT_θ，则有System_rad_θ·x＝RFT_θ；

2)在投影角度θ下，遍历每一个像素点(i,j)，计算其打在探测器上的绝对偏移量x_r和相对探测器像元位置radbin：

其中，为待重建图像的对角线像素点个数，radbin表示索引坐标为(i,j)的像素点在投影角度θ_i下打在第几个像元上。

3)将断层上的一点(i,j)与探测器像元radbin之间的投影映射关系存储到投影矩阵System_Rad_θ中，即有：

$\mathrm{System}\_{\mathrm{rad}}_{\mathrm{\&theta;}}(\mathrm{radbin},(i-1)\&CenterDot;N+j)={\mathrm{\&alpha;}}_{i,j,\mathrm{\&theta;},x_r}---\left(8\right)$

其中，为小于等于1的线段权值系数，其计算方法如式(3)所示。

4)当遍历完所有像素点后，投影矩阵System_Rad_θ构造完成。

投影矩阵System_Rad_θ中每一行代表着一条Radon投影射线，投影方程中的一点System_rad_θ(radbin,(i-1)·N+j)意味着列矢量化后的第(i-1)·N+j个像 素在第radbin条投影射线中的贡献程度，在Radon投影模型中，一般将投影贡献权值等比例的视为穿过该象素点的线段长值。

6.Radon系统和Mojette系统相互转化的可行性。

通过对比System_Rad_θ和System_Moj_p,q的构造方式得到，当Radon投影系统中的投影角度为θ＝tan^-1(q/p)，探测器像元物理尺寸为时，Radon投影射线之间的距离和Mojette投影射线之间的距离完全一致，Radon投影系统与Mojette投影系统的对比如图3所示，这时Radon投影射线的投影路径与Mojette投影射线路径一致，但与Mojette投影不同之处在于，Radon投影中还包括了穿过像素区域但没有穿过中心点的像素值的贡献，因此将Radon投影看作是由一组权重系数不同的Mojette投影线性累加构成，图3中的投影关系用代数形式来表达，有：

$\left\{\begin{array}{c}{\mathrm{PFT}}_{\mathrm{\&theta;}}\left(1\right)={\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(1,1)\\ {\mathrm{PFT}}_{\mathrm{\&theta;}}\left(2\right)={\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2,2)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,1}\right)\\ {\mathrm{PFT}}_{\mathrm{\&theta;}}\left(3\right)={\mathrm{Moj}}_{p,q}\left(3\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(3,3)+{\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(3,2)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{3,1}\right)\\ {\mathrm{PFT}}_{\mathrm{\&theta;}}\left(4\right)={\mathrm{Moj}}_{p,q}\left(4\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(4,4)+{\mathrm{Moj}}_{p,q}\left(3\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(4,3)+{\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{4,2}\right)\\ {\mathrm{PFT}}_{\mathrm{\&theta;}}\left(5\right)={\mathrm{Moj}}_{p,q}\left(5\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(5,5)+{\mathrm{Moj}}_{p,q}\left(4\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(5,4)+{\mathrm{Moj}}_{p,q}\left(3\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{5,3}\right)\\ {\mathrm{PFT}}_{\mathrm{\&theta;}}\left(6\right)={\mathrm{Moj}}_{p,q}\left(6\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(6,6)+{\mathrm{Moj}}_{p,q}\left(5\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(6,5)+{\mathrm{Moj}}_{p,q}\left(4\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(6,4)\\ {\mathrm{PFT}}_{\mathrm{\&theta;}}\left(7\right)={\mathrm{Moj}}_{p,q}\left(7\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(7,7)+{\mathrm{Moj}}_{p,q}\left(6\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(7,6)+{\mathrm{Moj}}_{p,q}\left(5\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{7,5}\right)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(k\right)={\mathrm{Moj}}_{p,q}\left(k\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k,k)+{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(k-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k-j,k-j)\end{array}\right.---\left(9\right)$

其中，α_p,q(i,j)表示Radon投影中的第i个分量和Mojette投影中第j个分量之间的依存关系。

但在实际的Radon投影中，探测器像元尺寸的值往往是固定的，并不能在每个投影角度下都满足的条件，若即Radon投影射线之间的距离大于Mojette投影射线之间的距离，此时的Radon 投影系统与Mojette投影系统的对比如图4所示，这时Radon投影射线的投影路径与Mojette投影射线路径不一致，Radon投影更为稀疏，图4中的投影关系用代数形式来表达，有：

$\left\{\begin{array}{c}{\mathrm{RFT}}_{\mathrm{\&theta;}}\left(1\right)={\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(2\right)={\mathrm{Moj}}_{p,q}\left(3\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,3}\right){\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(3\right)={\mathrm{Moj}}_{p,q}\left(5\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{3,5}\right)+{\mathrm{Moj}}_{p,q}\left(4\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{3,4}\right)+{\mathrm{Moj}}_{p,q}\left(3\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{3,3}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(4\right)={\mathrm{Moj}}_{p,q}\left(7\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{4,7}\right)+{\mathrm{Moj}}_{p,q}\left(6\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{4,6}\right)+{\mathrm{Moj}}_{p,q}\left(5\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{4,5}\right)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(k\right)={\mathrm{Moj}}_{p,q}(2k-1)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k,k)+{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(2k-1-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2k-1-j,2k-1-j)\end{array}\right.---\left(10\right)$

若Radon投影系统与Mojette投影系统的对比如图5所示。这时Radon投影射线的投影路径与Mojette投影射线路径不一致，Radon投影更为紧凑，虽然此时的Radon投影射线并不一定穿过每一组构成Mojette投影中像素点集的中心，但在这种更紧致条件下的Radon投影能穿过构成全部Mojette投影中的像素点集。图5中的投影关系用代数形式来表达，有：

$\left\{\begin{array}{c}{\mathrm{RFT}}_{\mathrm{\&theta;}}\left(1\right)={\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(2\right)={\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{{\mathrm{\&alpha;}}_{p,q}^{\&prime;}\left(\mathrm{1,1}\right)}_{}\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(3\right)={\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{3,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{3,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(4\right)={\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{4,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{4,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(5\right)={\mathrm{Moj}}_{p,q}\left(3\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{5,3}\right)+{\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{5,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{5,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(6\right)={\mathrm{Moj}}_{p,q}\left(3\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{6,3}\right)+{\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{6,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{6,1}\right)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}(2k-1)={\mathrm{Moj}}_{p,q}\left(k\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2k-1,k)+{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(k-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2k-1-j,k-j)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(2k\right)={\mathrm{Moj}}_{p,q}\left(k\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2k,k)+{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(k-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2k-j,k-j)\\ .\\ .\\ .\end{array}\right.---\left(11\right)$

设两投影之间转换矩阵为A_p,q∈R^J×B,使得RFT_θ＝A_p,q·Moj_p,q，其中，J＝DetRowNum代表Radon投影射线条数，也代表了探测器像元数量，B＝(N-1)|p|+(M-1)|q|+1代表了需要转换的Mojette投影总个数。该矩阵中的一个 分量α_p,q(i,j)表示Radon投影中的第i个分量和Mojette投影中第j个分量之间的关系，即若α_p,q(i,j)≠0，则说明RFT_θ(i)投影值中有Moj_p,q(j)中像素和的贡献；该矩阵中的一行表示该投影角度下Radon投影的第i个分量与所有Mojette投影之间的关系，即若α_p,q(i,j')≠0，则说明RFT_θ(i)投影值中有Moj_p,q(j')中像素和的贡献。

由RFT_θ＝A_p,q·Moj_p,q可知，当转换矩阵A_p,q可逆时，显然Radon投影是能够直接转换为Mojette投影的，但这并不意味着只有当转换矩阵A_p,q可逆时，Radon投影才能转化为Mojette投影。即|A_p,q|≠0是Radon投影转化为Mojette投影的充分条件，却不是转换的必要条件。

例如，在探测器分辨率较小的情况下，即时，构成同一个Mojette投影的像素点集内穿过了多条Radon射线，例如RFT_θ(1)＝Moj_p,q(1)·α_p,q(1,1)，RFT_θ(2)＝Moj_p,q(1)·α'_p,q(1,1)。这使得转换矩阵A_p,q中有多行线性相关，但由于存在用来转化为全部Mojette投影的完备Radon投影集，即存在满足下列条件的Radon投影：需要穿过构成Mojette投影的像素点集，且至少需要穿过两个相邻Mojette投影路径上的像素点的并集，第一条投影射线除外，如，RFT_θ(2)＝Moj_p,q(2)·α_p,q(2,2)+Moj_p,q(1)·α_p,q(2,1)，则RFT_θ(2)穿过了构成Mojette投影Moj_p,q(2)，Moj_p,q(1)中的像素点集，且穿过了2个彼此相邻的Mojette投影路径上的像素点的并集，当Moj_p,q(1)由RFT_θ(1)求出后，根据RFT_θ(2)的值求出Moj_p,q(2)。故而此处依然能够通过对线性方程的化简，将Radon投影转化为Mojette投影。即，当转换矩阵A_p,q经过线性变换化简至最简形式时，去掉其中的非零行，若其剩余代数余子式的行列式不为零，则该投影模型下的Radon投影能够转换为Mojette投影。

用数学条件来更准确地描述的话，即转换矩阵A_p,q中线性不相关的行矢量构 成的满秩矩阵的主对角线元素非零，或者说，当转换矩阵需要满足rank(A)≥(M-1)|p_i|+(N-1)|q_i|+1时，Radon投影满足转换为Mojette投影的完备性要求。

若探测器分辨率较大，即时，则采集的Radon投影较稀疏，根本不能满足上述的转换条件。故而，在投影矢量(p,q)的投影方向下，当时，对应投影角度θ＝tan^-1q/p,下的Radon投影都满足转换完备性的要求。

二、技术方案的步骤是：

S1设置投影参数，即探测器分辨率和重建图像分辨率；设探测器分辨率为B×B，重建图像分辨率为M×N，其中，重建图像分辨率和投影矢量(p,q)的选择需要满足max{(N-1)|p_i|+(M-1)|q_i|+1}<B。

S2选择投影矢量(p,q)，计算出相应的投影角度θ＝tan^-1q/p，并在选定的投影角度下采集Radon投影。

S3在投影矢量(p,q)下建立相应的0阶B样条插值投影系统System_rad_θ，需要明确两点，其一为穿过每个像素点的投影射线的投影位置，其二为穿过每个像素点的投影射线的线段权值。

S3.1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影射线打在探测器上的有效范围。

设重建像素的实际物理尺寸为ObjPixel，探测器像元的实际物理尺寸为DetCCDSize，当DetCCDSize<ObjPixel时，说明在一个矩形区域形状的离散像素点范围内，有不止一条射线穿过该像素点。遍历每一个像素(i,j)，其中有1<i<M,1<j<N，并计算出穿过该像素点的射线打在探测器上的绝对物理偏移值，此处用Xr_up和Xr_down来标记这个有效范围的最小值和最大值，并计算 出与Xr_up相对应的探测器像元标号最小值bin_up和与Xr_down相对应的探测器像元标号最大值bin_down，投影射线打在探测器上的有效范围的计算公式如式(12)所示：

$\left\{\begin{array}{c}\mathrm{Xr}\_\mathrm{up}=[((i-1)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+((j-1)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\\ \mathrm{Xr}\_\mathrm{down}[(\left(i\right)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+(\left(j\right)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\end{array}\right.,\mathrm{\&theta;}\&le;90$

$\left\{\begin{array}{c}\mathrm{Xr}\_\mathrm{up}=[(\left(i\right)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+((j-1)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\\ \mathrm{Xr}\_\mathrm{down}[((i-1)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+(\left(j\right)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\end{array}\right.,\mathrm{\&theta;}\&le;90---\left(12\right)$

其中，((i-1-M/2)·ObjPixel,(j-1-N/2)·ObjPixel)为离散化后的最小单位块的左上角那一临界点的空间坐标，((i-M/2)·ObjPixel,(j-N/2)·ObjPixel)为单位块的右下角那一临界点的空间坐标，((i-M/2)·ObjPixel,(j-1-N/2)·ObjPixel)为单位块的左下角那一临界点的空间坐标，((i-1-M/2)·ObjPixel,(j-N/2)·ObjPixel)为单位块的右上角那一临界点的空间坐标，为待重建图像的对角线像素点个数，如图6所示。

穿过像素中心点的射线打在探测器上的具体位置的计算公式，如式(13)所示。

x_r＝[((i-1/2)-M/2)·cos(θ)+((j-1/2)-N/2)·sin(θ)]·ObjPixel

                                                      (13) 

bin＝(x_r+D/2·ObjPixel)/DetCCDSize+1

其中，为待重建图像的对角线像素点个数。

S3.2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长度。设在所有投影角度下穿过像素中心点的射线的线段权值为1，而对穿过该像素其它位置的的射线，则需要根据这些射线与中心射线之间的垂直距离|ΔXr|，计算出其线段权重值α_p,q，其具体的计算公式如式(14)所示。

$\left\{\begin{array}{c}\mathrm{if}\left|\mathrm{\&Delta;Xr}\right|\&le;\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|-|\sin \mathrm{\&theta;}\left|\right)}{2},\\ {\mathrm{\&alpha;}}_{p,q}=1\\ \mathrm{if}\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|-|\sin \mathrm{\&theta;}\left|\right)}{2}<\left|\mathrm{\&Delta;}{x}_r\right|\&le;\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)}{2}\\ {\mathrm{\&alpha;}}_{p,q}=\frac{-2\&CenterDot;\left|\mathrm{\&Delta;Xr}\right|+\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)\&CenterDot;\mathrm{ObjPixel}}{2\&CenterDot;\mathrm{ObjPixel}\&CenterDot;\min \left(\right|\cos \mathrm{\&theta;},\sin \mathrm{\&theta;}\left|\right)}\\ \mathrm{if}\left|\mathrm{\&Delta;Xr}\right|>\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)}{2}\\ {\mathrm{\&alpha;}}_{p,q}=0\end{array}\right.---\left(14\right)$

其中，ObjPixel·(|cosθ|-|sinθ|)/2为一个临界距离值，当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离|ΔXr|小于等于这一个临界值时，穿过像素点的线段长度是最长的，将其记为1；与此类似，ObjPixel·(|cosθ|+|sinθ|)/2也是一个临界距离值，它标记着当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离|ΔXr|大于等于这一个临界值时，穿过像素点的线段长度是该投影角度下最短，已不在该像素点的范围之内，将其记为0。

并将计算出的权重值存入到系统矩阵System_rad_θ中相应的位置bin处，即  $\mathrm{System}\_{\mathrm{rad}}_{\mathrm{\&theta;}}(\mathrm{radbin},(i-1)\&CenterDot;N+j)={\mathrm{\&alpha;}}_{i,j,\mathrm{\&theta;},x_r},$ 其中N代表图像的列数。System_rad_θ中每一行代表一条具体的投影射线穿过图像中所有点的线段权值矢量，每一列代表所有投影射线穿过图像中一点的线段权值矢量。

S4根据已求得的系统矩阵System_rad_θ，来进行投影转换。

此处需要注意，由于一个像素里通常都有超过1条射线穿过，具体的射线条数取决于ρ＝ObjPixel/DetCCDSize，即像素物理大小和探测器实际大小的比值，故而对于Mojette采样来说，Radon采样就存在大量冗余，例如，仅穿过一个像素的Radon采样射线就不止一条，而此处只需取其中一条进行计算，即能够求出对应的Mojette投影。

因为系统方程System_rad_θ的一行代表该投影角度下的一条投影射线，选择射线其实就是在选择系统方程中的行矢量，选择这样的射线需要满足以下两个条件：1.该行矢量非零；2.该行矢量中存在未求出的像素值，即该行矢量包含已求出的像素之外的新像素，从而避免无用的重复计算。

S4.1导入由S3中生成的系统矩阵System_rad_θ，找出系统方程中第一个非零行，且该行矢量中只有一个非零分量，即意味着该射线只穿过了一个像素，投射在探测器上，将探测器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Moj_p,q(1)利用下式计算：

Moj_p,q(1)＝RFT_θ(1)/α_p,q(1,1)

其中，RFT_θ(1)表示Moj_p,q(1)的对应位置的Radon投影值，α_p,q(1,1)表示投影值为RFT_θ(1)的射线穿过的对应像素的线段权重值。

S4.2寻找用来求出下一个Mojette投影的Radon投影射线。在S4.1已求出的Mojette投影Moj_p,q(1)的基础上，若想求出下一个Mojette投影Moj_p,q(2)，需要找到形如RFT_θ(2)＝Moj_p,q(2)·α_p,q(2,2)+Moj_p,q(1)·α_p,q(2,1)的Radon投影，即，由且仅由已求得的Mojette投影Moj_p,q(1)和待求Mojette投影Moj_p,q(2)构成的Radon投影，若所有的Radon投影都不满足这一条件，则该待求的Mojette投影Moj_p,q(2)无法被求出。

S4.2.1首先，找到待求Mojette投影Moj_p,q(2)中穿过的一个像素点的索引坐标(i,j)，并计算出其在列矢量化后的向量x中的列标l＝(i-1)·N+j。

S4.2.2找到系统矩阵System_rad_θ的第l列，遍历System_rad_θ(:,l)，找出其中非零值对应的行数范围，这相当于遍历所有的投影射线，并找出穿过(i,j)这一点，权重系数值不为零的投影射线，设这些投影矢量在投影矩阵中的行标最小值为bin_up，最大值为bin_down。在这一范围内，遍历每一条投影射线在投影矩 阵System_rad_θ中对应的横行，即在bin_up≤bin≤bin_down范围内遍历System_rad_θ(bin,:)，若系统矩阵行System_rad_θ(bin,:)中的非零系数处对应的像素点集，由且仅由Moj_p,q(1)和Moj_p,q(2)中穿过的像素点构成，即RFT_θ(2)＝Moj_p,q(2)·α_p,q(2,2)+Moj_p,q(1)·α_p,q(2,1)，则该投影射线即为候选的更新射线。

S4.2.3在满足这一要求的射线当中，选择α_p,q(2,2)值最大的Radon投影射线作为下一个用来更新Mojette投影的射线。

S4.3利用求得的Mojette投影来更新Radon投影，即在该Radon投影中减去已求得的Mojette投影，再除以线段权值，就得到了新的Mojette投影：Moj_p,q(2)＝(RFT_θ(2)-Moj_p,q(1)·α_p,q(2,1))/α_p,q(2,2)。

重复这一步骤，直到将所有的Mojette投影求出。 

即以式(15)为标准，寻找用来更新Mojette投影的Radon投影：

$\left\{\begin{array}{c}{\mathrm{RFT}}_{\mathrm{\&theta;}}\left(1\right)={\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(2\right)={\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2,1)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(k\right)={\mathrm{Moj}}_{p,q}\left(k\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k,k)+{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(k-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k-j,k-j)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(\mathrm{DetRN}\right)={\mathrm{\&Sigma;}}_{k=1}^{k\&le;\mathrm{DetRN}}{\mathrm{Moj}}_{p,q}\left(k\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k,k)\end{array}\right.---\left(15\right)$

其中， ${\mathrm{Moj}}_{p,q}\left(k\right)=\underset{k=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{l=0}{\overset{Q-1}{\mathrm{\&Sigma;}}}\mathrm{pixel}(i,j)\mathrm{\&delta;}(k-p_{i}i+q_{i}j).$

基于已经求出的Mojette投影，从式(15)中推导出各Mojette投影的转换公式，即如式(16)所示：

$\left\{\begin{array}{c}{\mathrm{Moj}}_{p,q}\left(1\right)={\mathrm{RFT}}_{\mathrm{\&theta;}}\left(1\right)/{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{Moj}}_{p,q}\left(2\right)=({\mathrm{RFT}}_{\mathrm{\&theta;}}\left(2\right)-{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(2,1))/{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,2}\right)\\ .\\ .\\ .\\ {\mathrm{Moj}}_{p,q}\left(k\right)=({\mathrm{RFT}}_{\mathrm{\&theta;}}\left(k\right)-{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(k-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k-j,k-j))/{\mathrm{\&alpha;}}_{p,q}(k,k)\\ .\\ .\\ .\end{array}\right.---\left(16\right)$

至此，已将投影角度θ＝tan^-1q/p下的Radon投影转化为投影矢量(p,q)下的Mojette投影。

本发明的效果和益处是在不同的稀疏投影角度下，不需要改变Radon投影图像的分辨率大小，克服了Mojette投影域分辨率随着投影矢量变化和Radon投影分辨率不变之间的转换障碍，搭建起这两种投影之间的转换桥梁。

附图说明

图1是待重建物体所处的连续域空间坐标图。图中，1为固定坐标系，2为旋转坐标系。重建物体所在的二维固定坐标系用x-o-y表示，坐标原点o为物体的几何中心。模拟光源与探测器在x-o-y平面内逆时针旋转，其旋转坐标系用x_r-o-y_r表示。

图2是离散域索引坐标图，即将待重建物体所处的区域划分成等间距的正方行像素块，离散化后的重建区域中每一个像素点的行列位置通过索引坐标(i,j)来标记，i代表行数，j代表列数。

图3是当探测器分辨率时，Radon投影与Mojette投影的关系示意图，图中，(a)为Radon投影过程；(b)为Mojette投影过程。

图4是当探测器分辨率时，Radon投影与Mojette投影的关系示意图，图中，(a)为Radon投影过程；(b)为Mojette投影过程。

图5是当探测器分辨率时，Radon投影与Mojette投影的关系示意图，图中，(a)为Radon投影过程；(b)为Mojette投影过程。

图6是穿过一个像素点的全部Radon投影射线示意图，图中，(a)为投影角度小于90°时，穿过该像素点的射线打在探测器上的有效范围。图中，(b)为投影角度大于90°时，穿过该像素点的射线打在探测器上的有效范围。

具体实施方式

下面结合技术方案和附图详细说明本发明的具体实施方式：

S1设置投影参数,设在投影矢量(1,2)下，对一个3×3的小块进行投影，设探测器像元个数B＝10>(3-1)|1|+(3-1)|2|+1＝7；

S2投影矢量为(1,2)时，相应的投影角度θ＝tan^-11/2＝26.5651°；

S3在投影矢量(1,2)下建立相应的0阶B样条插值投影系统System_rad_26.6，此处设DetCCDSize＝0.4mm。

S3.1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影射线打在探测器上的有效范围。

设重建像素的实际物理尺寸为ObjPixel，ObjPixel＝1mm；探测器像元的实际物理尺寸为DetCCDSize，DetCCDSize＝0.4mm。

遍历每一个像素{Pixel(i,j)|i,j∈Z,1<i<3,1<j<3}，并计算出穿过每一个像素点的射线打在探测器上的有效范围，根据式(6)计算得出，穿过(1,1)该点的探测器像元有效范围为1～4，即在(1,2)的投影矢量下，穿过第一个像素小块的各射线分别打在第1～4个探测器像元上；穿过(1,2)该点的射线分别打在第4～6个探测器像元上；穿过(1,3)该点的射线分别打在第6～8个探测器像元上；穿过(2,1)该点的射线分别打在第2～5个探测器像元上，依次类推。

S3.2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长度。根据式(7)计算得出各权重的值，并存储至System_rad_26.6相应的位置处。

经过计算得到的系统矩阵如下所示：

$\left[\begin{array}{ccccccccc}0.2038& 0& 0& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0.098& 0.& 0& 0& 0& 0\\ 1& 0& 0& 0.993& 0& 0& 0& 0& 0\\ 0.113& 0.887& 0& 1& 0& 0& 0.887& 0& 0\\ 0& 1& 0& 0.218& 0.782& 0& 1& 0& 0\\ 0& 0.324& 0.676& 0& 1& 0& 0.324& 0.676& 0\\ 0& 0& 1& 0& 0.430& 0.570& 0& 1& 0\\ 0& 0& 0.535& 0& 0& 1& 0& 0.535& 0.465\\ 0& 0& 0& 0& 0& 0.641& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0.746\end{array}\right]\&CenterDot;x=\mathrm{RFT}$

S4根据已求得的系统矩阵System_rad_θ，来进行投影转换。

S4.1导入由S3中生成的系统矩阵System_rad_θ，找出系统方程中第一个非零行，且该行矢量中只有一个非零分量，意味着该射线只穿过了一个像素，投射在探测器上，将探测器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Moj_p,q(1)利用下式计算：

Moj_p,q(1)＝RFT_θ(1)/0.2038

S4.2寻找用来求出下一个Mojette投影的Radon投影射线。 

S4.2.1首先，找到构成Moj_p,q(2)的像素点的索引坐标位置，此处为(2,1)，并计算出其在列矢量化后的向量中的位置l＝(2-1)·3+1＝4。

S4.2.2找到系统矩阵System_rad_θ的第4列，遍历System_rad_θ(:,4)，找出其中非零值所在的行数范围，设最小值为bin_up，最大值为bin_down，此处，bin_up＝2，bin_down＝5，在2≤bin≤5范围内遍历System_rad_θ(bin,:)，观察其投影路径上穿过的像素点的位置。若该系统矩阵行中的非零系数对应位置处的像素点，由且仅由Mojette投影射线Moj_p,q(1)和Moj_p,q(2)中穿过的像素点构成，此处，点(1,1)(在矩阵中第1列)构成Moj_p,q(1)，点(2,1)(在矩阵中第4列)构成Moj_p,q(2)，若有RFT_θ(bin)＝Moj_1,2(1)·α_1,2(bin,1)+Moj_1,2(2)·α_1,2(bin,2)，则该投影射线 System_rad_θ(bin,:)即为候选的更新射线。此处投影射线System_rad_θ(2,:)和System_rad_θ(3,:)均满足要求，而System_rad_θ(4,:)和System_rad_θ(5,:)中还包含了Mojette投影射线Moj_1,2(3)中穿过的像素点，因此不能作为求出Mojette投影Moj_p,q(2)的候选Radon射线。

S4.2.3在满足这一要求的射线当中，选择α_1,2(bin,2)值最大的Radon投影射线作为下一个用来更新Mojette投影的射线，即选择System_rad_θ(3,:)这条射线得到的Radon投影来进行求解，此时的α_1,2(3,2)＝0.9927。

S4.3利用求得的Mojette投影来更新Radon投影，即在该Radon投影中减去已求得的Mojette投影，再除以线段权值，就得到了新的Mojette投影Moj_p,q(2)＝(RFT_θ(3)-Moj_p,q(1)·1)/0.9927。

重复这一步骤，直到将所有的Mojette投影求出。 

此处的转换矩阵A_1,2如下所示：

$\left[\begin{array}{cccccccc}0.2038& 0& 0& 0& 0& 0& 0& 0\\ 1& 0.0982& 0& 0& 0& 0& 0& 0\\ 1& 0.9927& 0& 0& 0& 0& 0& 0\\ 0.1129& 1& 0.8871& 0& 0& 0& 0& 0\\ 0& 0& 0.3241& 1& 0.6759& 0& 0& 0\\ 0& 0& 0& 0.4296& 1& 0.5704& 0& 0\\ 0& 0& 0& 0& 0.5352& 1& 0.4648& 0\\ 0& 0& 0& 0& 0& 0& 0.6408& 1\end{array}\right]\&CenterDot;\left[\begin{array}{c}\mathrm{Moj}\left(1\right)\\ \mathrm{Moj}\left(1\right)\\ \mathrm{Moj}\left(2\right)\\ \mathrm{Moj}\left(3\right)\\ \mathrm{Moj}\left(4\right)\\ \mathrm{Moj}\left(5\right)\\ \mathrm{Moj}\left(6\right)\\ \mathrm{Moj}\left(7\right)\end{array}\right]=\mathrm{RFT}$

各Mojette投影的求解过程如下：

$\left\{\begin{array}{c}{\mathrm{Moj}}_{p,q}\left(1\right)={\mathrm{RFT}}_{\mathrm{\&theta;}}\left(1\right)/0.2038\\ {\mathrm{Moj}}_{p,q}\left(2\right)=({\mathrm{RFT}}_{\mathrm{\&theta;}}\left(3\right)-{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;1)/0.9927\\ {\mathrm{Moj}}_{p,q}\left(3\right)=({\mathrm{RFT}}_{\mathrm{\&theta;}}\left(4\right)-{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;0.1129-{\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;1)/0.8871\\ .\\ .\\ .\\ {\mathrm{Moj}}_{p,q}\left(7\right)=({\mathrm{RFT}}_{\mathrm{\&theta;}}\left(7\right)-{\mathrm{Moj}}_{p,q}\left(6\right)\&CenterDot;0.6408)/1\\ .\\ .\\ .\end{array}\right.$

至此，已将投影角度θ＝tan^-1q/p下的Radon投影转化为投影矢量(p,q)下的Mojette投影。

Claims (1)

1.基于固定分辨率条件下的离散Radon和Mojette投影转换方法，以计算机层析成像算法为基础，通过构造合理的成像系统来逼近Radon和Mojette投影场景，并在详细分析了Mojette投影和Radon投影的关系之后，给出了在固定分辨率前提下各投影角度下的Radon投影转化为对应的投影矢量下的Mojette投影的充要条件和具体的转化算法，其特征是，步骤如下：

S1设置投影参数，即探测器分辨率和重建图像分辨率；设探测器分辨率为B×B，重建图像分辨率为M×N，其中，重建图像分辨率和投影矢量(p,q)的选择需要满足max{(N-1)|p_i|+(M-1)|q_i|+1}<B；

S2选择投影矢量(p,q)，计算出相应的投影角度θ＝tan^-1q/p，并在选定的投影角度下采集Radon投影；

S3在投影矢量(p,q)下建立相应的0阶B样条插值投影系统System_rad_θ，需要明确两点，其一为穿过每个像素点的投影射线的投影位置，其二为穿过每个像素点的投影射线的线段权值；

S3.1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影射线打在探测器上的有效范围；

设重建像素的实际物理尺寸为ObjPixel，探测器像元的实际物理尺寸为DetCCDSize，当DetCCDSize<ObjPixel时，说明在一个矩形区域形状的离散像素点范围内，有不止一条射线穿过该像素点；遍历每一个像素(i,j)，其中有1<i<M,1<j<N，并计算出穿过该像素点的射线打在探测器上的绝对物理偏移值，此处用Xr_up和Xr_down来标记这个有效范围的最小值和最大值，并计算出与Xr_up相对应的探测器像元标号最小值bin_up和与Xr_down相对应的探测器像元标号最大值bin_down，投影射线打在探测器上的有效范围的计算公式如式(12)所示：

$\left\{\begin{array}{c}\mathrm{Xr}\_\mathrm{up}=[((i-1)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+((j-1)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\\ \mathrm{Xr}\_\mathrm{down}=[(\left(i\right)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+(\left(j\right)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\end{array}\right.,\mathrm{\&theta;}\&le;90$

$\left\{\begin{array}{c}\mathrm{Xr}\_\mathrm{up}=[(\left(i\right)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+((j-1)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\\ \mathrm{Xr}\_\mathrm{down}=[((i-1)-M/2)\&CenterDot;\cos \left(\mathrm{\&theta;}\right)+(\left(j\right)-N/2)\&CenterDot;\sin \left(\mathrm{\&theta;}\right)]\&CenterDot;\mathrm{ObjPixel}\end{array}\right.\mathrm{\&theta;}>90---\left(12\right)$

其中，((i-1-M/2)·ObjPixel,(j-1-N/2)·ObjPixel)为离散化后的最小单位块的左上角那一临界点的空间坐标，((i-M/2)·ObjPixel,(j-N/2)·ObjPixel)为单位块的右下角那一临界点的空间坐标，((i-M/2)·ObjPixel,(j-1-N/2)·ObjPixel)为单位块的左下角那一临界点的空间坐标，((i-1-M/2)·ObjPixel,(j-N/2)·ObjPixel)为单位块的右上角那一临界点的空间坐标，为待重建图像的对角线像素点个数，如图6所示；

穿过像素中心点的射线打在探测器上的具体位置的计算公式，如式(13)所示：

x_r＝[((i-1/2)-M/2)·cos(θ)+((j-1/2)-N/2)·sin(θ)]·ObjPixel   (13)

bin＝(x_r+D/2·ObjPixel)/DetCCDSize+1

其中，为待重建图像的对角线像素点个数；

S3.2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长度；设在所有投影角度下穿过像素中心点的射线的线段权值为1，而对穿过该像素其它位置的的射线，则需要根据这些射线与中心射线之间的垂直距离|ΔXr|，计算出其线段权重值α_p,q，其具体的计算公式如式(14)所示：

$\left\{\begin{array}{c}\mathrm{if}\left|\mathrm{\&Delta;Xr}\right|\&le;\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|-|\sin \mathrm{\&theta;}\left|\right)}{2},\\ {\mathrm{\&alpha;}}_{p,q}=1\\ \mathrm{if}\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|-|\sin \mathrm{\&theta;}\left|\right)}{2}<\left|\mathrm{\&Delta;Xr}\right|\&le;\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)}{2}\\ {\mathrm{\&alpha;}}_{p,q}=\frac{-2\&CenterDot;\left|\mathrm{\&Delta;Xr}\right|+\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)\&CenterDot;\mathrm{ObjPixel}}{2\&CenterDot;\mathrm{ObjPixel}\&CenterDot;\min \left(\right|\cos \mathrm{\&theta;},\sin \mathrm{\&theta;}\left|\right)}\\ \mathrm{if}\left|\mathrm{\&Delta;Xr}\right|>\frac{\mathrm{ObjPixel}\&CenterDot;\left(\right|\cos \mathrm{\&theta;}|+|\sin \mathrm{\&theta;}\left|\right)}{2}\\ {\mathrm{\&alpha;}}_{p,q}=0\end{array}\right.---\left(14\right)$

其中，ObjPixel·(|cosθ|-|sinθ|)/2为一个临界距离值，当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离|ΔXr|小于等于这一个临界值时，穿过像素点的线段长度是最长的，将其记为1；与此类似，ObjPixel·(|cosθ|+|sinθ|)/2也是一个临界距离值，它标记着当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离|ΔXr|大于等于这一个临界值时，穿过像素点的线段长度是该投影角度下最短，已不在该像素点的范围之内，将其记为0；

并将计算出的权重值存入到系统矩阵System_rad_θ中相应的位置bin处，即 ${\mathrm{System}\_\mathrm{rad}}_{\mathrm{\&theta;}}(\mathrm{radbin},(i-1)\&CenterDot;N+j)={\mathrm{\&alpha;}}_{i,j,\mathrm{\&theta;},x_r},$ 其中N代表图像的列数；System_rad_θ中每一行代表一条具体的投影射线穿过图像中所有点的线段权值矢量，每一列代表所有投影射线穿过图像中一点的线段权值矢量；

S4根据已求得的系统矩阵System_rad_θ，来进行投影转换；

此处需要注意，由于一个像素里通常都有超过1条射线穿过，具体的射线条数取决于ρ＝ObjPixel/DetCCDSize，即像素物理大小和探测器实际大小的比值，故而对于Mojette采样来说，Radon采样就存在大量冗余，例如，仅穿过一个像素的Radon采样射线就不止一条，而此处只需取其中一条进行计算，即求出了对应的Mojette投影；

因为系统方程System_rad_θ的一行代表该投影角度下的一条投影射线，选择射线其实就是在选择系统方程中的行矢量，选择这样的射线需要满足以下两个条件：1.该行矢量非零；2.该行矢量中存在未求出的像素值，即该行矢量包含已求出的像素之外的新像素，从而避免无用的重复计算；

S4.1导入由S3中生成的系统矩阵System_rad_θ，找出系统方程中第一个非零行，且该行矢量中只有一个非零分量，即意味着该射线只穿过了一个像素，投射在探测器上，将探测器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Moj_p,q(1)利用下式计算：

Moj_p,q(1)＝RFT_θ(1)/α_p,q(1,1)

其中，RFT_θ(1)表示Moj_p,q(1)的对应位置的Radon投影值，α_p,q(1,1)表示投影值为RFT_θ(1)的射线穿过的对应像素的线段权重值；

S4.2寻找用来求出下一个Mojette投影的Radon投影射线；在S4.1已求出的Mojette投影Moj_p,q(1)的基础上，若想求出下一个Mojette投影Moj_p,q(2)，需要找到形如RFT_θ(2)＝Moj_p,q(2)·α_p,q(2,2)+Moj_p,q(1)·α_p,q(2,1)的Radon投影，即，由且仅由已求得的Mojette投影Moj_p,q(1)和待求Mojette投影Moj_p,q(2)构成的Radon投影，若所有的Radon投影都不满足这一条件，则该待求的Mojette投影Moj_p,q(2)无法被求出；

S4.2.1首先，找到待求Mojette投影Moj_p,q(2)中穿过的一个像素点的索引坐标(i,j)，并计算出其在列矢量化后的向量x中的列标l＝(i-1)·N+j；

S4.2.2找到系统矩阵System_rad_θ的第l列，遍历System_rad_θ(:,l)，找出其中非零值对应的行数范围，这相当于遍历所有的投影射线，并找出穿过(i,j)这一点，权重系数值不为零的投影射线，设这些投影矢量在投影矩阵中的行标最小值为bin_up，最大值为bin_down；在这一范围内，遍历每一条投影射线在投影矩阵System_rad_θ中对应的横行，即在bin_up≤bin≤bin_down范围内遍历System_rad_θ(bin,:)，若系统矩阵行System_rad_θ(bin,:)中的非零系数处对应的像素点集，由且仅由Moj_p,q(1)和Moj_p,q(2)中穿过的像素点构成，即RFT_θ(2)＝Moj_p,q(2)·α_p,q(2,2)+Moj_p,q(1)·α_p,q(2,1)，则该投影射线即为候选的更新射线；

S4.2.3在满足这一要求的射线当中，选择α_p,q(2,2)值最大的Radon投影射线作为下一个用来更新Mojette投影的射线；

S4.3利用求得的Mojette投影来更新Radon投影，即在该Radon投影中减去已求得的Mojette投影，再除以线段权值，就得到了新的Mojette投影：

Moj_p,q(2)＝(RFT_θ(2)-Moj_p,q(1)·α_p,q(2,1))/α_p,q(2,2)；

重复这一步骤，直到将所有的Mojette投影求出；

即以式(15)为标准，寻找用来更新Mojette投影的Radon投影：

$\left\{\begin{array}{c}{\mathrm{RFT}}_{\mathrm{\&theta;}}\left(1\right)={\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(2\right)={\mathrm{Moj}}_{p,q}\left(2\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,1}\right)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(k\right)={\mathrm{Moj}}_{p,q}\left(k\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k,k)+{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(k-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k-j,k-j)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\mathrm{\&theta;}}\left(\mathrm{DetRN}\right)={\mathrm{\&Sigma;}}_{k=1}^{k\&le;\mathrm{DetRN}}{\mathrm{Moj}}_{p,q}\left(k\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k,k)\end{array}\right.---\left(15\right)$

其中， ${\mathrm{Moj}}_{p,q}\left(k\right)=\underset{k=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{l=0}{\overset{Q-1}{\mathrm{\&Sigma;}}}\mathrm{pixel}(i,j)\mathrm{\&delta;}(k-p_{i}i+q_{i}j);$

基于已经求出的Mojette投影，从式(15)中推导出各Mojette投影的转换公式，即如式(16)所示：

$\left\{\begin{array}{c}{\mathrm{Moj}}_{p,q}\left(1\right)={\mathrm{RFT}}_{\mathrm{\&theta;}}\left(1\right)/{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{Moj}}_{p,q}\left(2\right)=({\mathrm{RFT}}_{\mathrm{\&theta;}}\left(2\right)-{\mathrm{Moj}}_{p,q}\left(1\right)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,1}\right))/{\mathrm{\&alpha;}}_{p,q}\left(\mathrm{2,2}\right)\\ .\\ .\\ .\\ {\mathrm{Moj}}_{p,q}\left(k\right)=({\mathrm{RFT}}_{\mathrm{\&theta;}}\left(k\right)-{\mathrm{\&Sigma;}}_{j=1}^{j\&le;k-1}{\mathrm{Moj}}_{p,q}(k-j)\&CenterDot;{\mathrm{\&alpha;}}_{p,q}(k-j,k-j))/{\mathrm{\&alpha;}}_{p,q}(k,k)\\ .\\ .\\ .\end{array}\right.---\left(16\right)$

至此，已将投影角度θ＝tan^-1q/p下的Radon投影转化为投影矢量(p,q)下的Mojette投影。