基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法
基于固定分辨率条件下的离散Radon投影和Mojette投影转换方法
【技术领域】
[0001] 本发明属于图像处理领域里的计算机层析成像（ComputedTomography，CT)技术 领域，涉及在图像投影变换域的稀疏采样及高效复原问题。
【背景技术】
[0002] CT技术是一种基于计算机的三维重构技术，能够在不损坏物体表层和内部构造的 前提下，重现出物体不可见的内部切片结构。在CT成像算法上，按照成像原理的不同，分为 解析类和迭代类重建算法；按照重建精确性的不同，又分为近似重建和精确重建算法。在解 析算法中，奥地利数学学家Radon提出的变换模型利用低维的压缩投影数据来重构高维度 的物体结构，这一数学模型的提出对CT的图像重建技术起到了关键性的启蒙和指导作用。 在这一理论基础的指引下，滤波反投影等解析算法相继被提出，该算法在完备的投影角度 下，能精确的重建出图像断层，且能抑制在重建过程中产生的各种伪影，在CT成像领域具 有里程碑式的意义。
[0003] 但由于成像系统是数字离散系统，因而连续的投影模型无法适用于离散的图像， 故而离散Radon变换随机应运而生。离散Radon变换及其逆变换解决了从模拟域向数字域 转换的问题，但基于Radon逆变换的重建算法需要大量的采样样本和投影数量，才能重建 出较好的重建断层，这种重建需求必然会导致过高的辐射剂量和过长的重建时间，出于对 医学成像中病人身体健康和病灶动态变化会导致伪影的考虑，在短时间、低剂量和稀疏角 度下进行高质量成像就显得尤为重要，同时又充满挑战。
[0004] 在Radon变换的基础上，GuedonJP等提出了Mojette变换的概念，该变换是离 散Radon变换的一种特殊形式，Mojette变换通过改变不同投影矢量下的采样率，能够在最 大程度上避免像素点重复和冗余采样，以避免出现一部分像素点过采样，而另一部分像素 点欠采样而带来的重建效果上的不精确。并通过充分利用已重建出来的像素点的信息，从 而大大减小了重建断层所需的投影角度和投影射线条数。
[0005] 但由于Mojette投影系统要求探测器分辨率随着投影角度的变化而变化，且在投 影角度之间的步进角也不是一个固定的常值，这与现行的采样模式大相径庭，实际的投影 系统符合Radon变换的采样模式，探测器分辨率不随投影角度而变化，因此通过Radon变换 得到的投影值不可能直接参与到Mojette逆变换当中来，而是先要将Radon投影转换到相 关的Mojette投影上去，再利用Mojette逆变换来进行反重建。
[0006] 完成从Radon投影到Mojette投影的转换，其中需要解决几个实际的问题：1.投 影矢量和投影角度之间的匹配转换，由于Mojette投影系统中采用一对互质的整数来表达 投影方向，称之为投影矢量，而Radon投影系统中采用的是极坐标中的旋转角度，故而需要 通过简单的换算将Mojette投影矢量表达的方向转换成Radon投影中的旋转角；2.Mojette 变换中可变分辨率与Radon变换中固定分辨率之间的矛盾，由于Mojette变换中的探测器 分辨率或者说是投影射线条数为随着投影矢量而变化的函数，而Radon投影的探测器分辨 率是不随着投影角度而变化的，这就使得采样系统中得到的Radon投影需要通过本发明中 提出的算法进行转换，才能得到可以高效重建图像的Mojette投影。
【发明内容】
[0007] 本发明的目的是提供了一种基于固定分辨率采集的Radon投影转化为Mojette 投影的算法，该方法在详细分析了Mojette投影和Radon投影的关系之后，解决了固定分辨 率条件下的两投影之间的转化难题，给出了Radon投影转化为对应投影矢量下的Mojette 投影的充要条件，并基于此给出了具体的转化算法，使得实际成像条件下采集的Radon投 影能够转化为Mojette域的投影。
[0008] 本发明的技术方案的原理及步骤如下：
[0009] -、技术方案的原理是：
[0010] 1.建立投影空间坐标系：设投影空间坐标系用x-o-y表示，如图1所示，坐标原 点〇为待重建物体的几何中心，X轴正向与探测器像元索引号增加方向一致，y轴正向与 射线传播方向一致，x-o-y平面内的一点的空间坐标为（x，y)。光源与探测器在x-o-y平 面内逆时针旋转，旋转角度为0，设旋转后的空间坐标系用\-o_h表示，像素点的旋转坐 标为在连续域的二维图像空间内，将物体所在区域离散化为MXN的离散小块， 设ObjSize为待重建矩形区域的行方向上的边长，则每一个重建像素的物理尺寸大小为 ObjPixel=ObjSiZe/M，如图2所示。其中，（i，j)为标记离散小块行列位置的索引坐标， 即若利用二维矩阵存储该离散图像，则i代表在矩阵中的行数，j代表在矩阵中的列数，通 常将图像左上角第一个小块的位置记为索引坐标起始位置（1，1)。
[0011] 明确重建断层像素在x-〇-y平面内的空间坐标（X，y)和在旋转坐标系Xj-o-yj?平 面内旋转坐标0^，^)之间的关系，如式（1)所示。
[0013] 明确重建断层像素在x-0-y平面内的空间坐标（x，y)和离散域索引坐标（i，j)之 间的关系，通常，以待重建断层的左上角像素为起点，记该像素的索引坐标为（1，1)，索引坐 标为（1，1)的像素中心点的空间坐标（x，y)为：
[0014] ((l-(M+l)/2) ?ObjPixel，（l-(N+l)/2) ?ObjPixel)
[0015] 索引坐标为（i，j)的像素中心点的空间坐标（x，y)为：
[0016] ((i-(M+l)/2) ?ObjPixel, (j-(N+l)/2) ?ObjPixel)
[0017] 基于以上两种坐标的对应，则得到索引坐标为（i，j)的像素中心点在旋转坐标系 中的旋转坐标为：
[0018] xr= [(j-(N+l)/2) ?ObjPixel?sin( 0 ) + (i-(M+l)/2) ?ObjPixel?cos( 0 )]
[0019] yr= [-(i-(M+l)/2) ?ObjPixel?sin( 0 ) + (j-(N+l)/2) ?ObjPixel?cos( 0 )]
[0020] 2.离散Radon变换：对分辨率大小为MXN的离散图像f进行离散Radon变换，其 过程用式（2)来表述：
[0022] 其中，1^1\表示在投影角度0下的Radon投影，RFT0 〇〇表示在投影角度0下 打在探测器x,位置处的Radon投影值。f(i，j)代表待重建的图像切片上索引坐标为（i，j) 的一点的灰度值。若以（0，\)标记投影角度0下打在探测器上的物理偏移量为&的一 条投影射线，％，从^表示在投影射线（0，\)和象素点（i，j)之间的权重核函数，通常为 S函数或0次样条函数（即将投影贡献权值等比例的视为穿过该象素点的线段长值）。
[0023] 设在所有投影角度0下穿过像素中心点的射线的线段权值为1，而对穿过该像素 其它位置的的射线，则需要根据这些射线与中心射线之间的垂直距离IAX」，计算出其线 段权重值％，从，&，其具体的计算公式如式（3)所示。
[0025] 其中，ObjPixel代表每一个重建像素的物理尺寸。
[0026] 3.Mojette变换：Mojette变换中用一对互质的整数（p,q)来表达投影方向，一般 有pGZ，qeZ+，p代表了图像列方向上的整数位移，q代表了图像行方向上的整数位移， 投影矢量（P，q)表达的投影角度0 =tarTHq/p)，对分辨率大小为MXN的离散图像f进 行Mojette变换，其过程用式（4)来表述：
[0028] 其中，1\1〇几(1表示在投影矢量（p,q)下的Mojette投影，Moj^(bin)表示在投影矢 量（P，q)下打在探测器bin上Mojette投影值。f(i，j)代表待重建的图像切片上索引坐标 为（i，j)的一点的灰度值，P〇为一个探测器像元位置矫正量，当投影矢量P>〇时，P〇=〇; 当投影矢量P〈〇时，P〇= (N-1) ?p。
[0029] Mojette变换与Radon变换相比最大的不同在于，Radon变换中探测器分辨率为 固定的常值，穿过物体切片的Radon投影射线之间距离是固定的，因而不能保证每条Radon 投影射线都能穿过像素的中心点，穿过像素中心点的线段权值记为1，而不穿过像素中心点 的线段权值要根据具体线段长度而定，计算方法如式（3)所示；而Mojette变换中每条投影 射线都只穿过像素的中心点，因此不同投影矢量（P，q)下探测器分辨率B(M，N，p，q)不同， 探测器分辨率的值是由图像重建分辨率MXN和投影矢量（p，q)共同决定的，即投影矢量 (p,q)下的Mojette探测器分辨率为B(P,Q,p,q) = (Q-1) |p| + (P-1) |q|+1。这意味着在覆 盖相同直径范围内的切片时，不同投影角度下的Mojette投影射线之间的间距hi不同：
[0031] 此外，Radon变换在模拟域重建是精确的，在离散域重建是近似的。而Mojette变 换在离散域的重建是精确的。
[0032] 4.Mojette投影矩阵的构造方法：
[0033] 基于Mojette正变换原理，下面阐述一下Mojette投影矩阵的构造方法。
[0034] 1)初始化参数，设重建断层分辨率为MXN，投影矢量为（p，q)，投影矩阵为
'，其中，B= (N-l)|Pi| + (M-l)|qi|+l，V=M.N。待重建图像列矢 量化后的一维向量为X，投影矢量（p,q)下的Mojette投影为Mojj^，则有Systen^Mojp^*x=M〇jp；q；
[0035] 2)在每个投影矢量（p，q)下，遍历每一个像素点（i，j)，计算其打在探测器上的位 置bin:
[0037] bin=p? (i_l)+q? (j_l)-po+1
[0038] 3)将通过投影矢量（Pi，qi)建立起来的断层上的一点（i，j)与探测器像元bin之 间的投影映射关系，存储到线性投影矩阵SyStem_M〇jp，q中，即：
[0039] System_Mojp；q(bin, (i~l) ?N+j) = 1 ；
[0040] 4)当遍历完所有像素点后，投影矩阵System_Mojp;(1构造完成。
[0041] 投影矩阵Systen^Moj^if每一行代表着一条Mojette投影射线，每一列 代表着穿过一个像素点的所有投影射线的线段权值。投影方程中的一 点SyStem_ Mojp，q(bin，（i-1) *N+j)意味着列矢量化后的第（i-1) *N+j个像素在第bin条投影射线中 的贡献程度，在Mojette投影模型中，若第bin条投影射线穿过第（i-1) *N+j个像素的中 心点，则System-MoUbin，（i-1) .N+j) = 1，否则为 0。
[0042] 5.Radon投影矩阵的构造方法：
[0043] 在Radon投影当中，线性投影矩阵的构造方法是由探测器和射线源与待重建物体 之间的几何关系决定的。
[0044] 1)设重建断层分辨率为MXN，包含待重建断层的最小正方形区域的边长为 ObjSize，每一个重建像素的实际物理尺寸ObjPixel=ObjSize/M。设探测器分辨率 为DetRowNum，其中每一个探测器像元的物理尺寸为DetCXDSize，投影矩阵为System_ rad0GRJ=DetRowNum，V=M*N。待重建图像列矢量化后的一维向量为x，投影角度 Q下的Radon投影为RFT0，贝丨」有System_rad0 ?x=RFT0;
[0045] 2)在投影角度0下，遍历每一个像素点（i，j)，计算其打在探测器上的绝对偏移 量Xj?和相对探测器像元位置radbin:
[0046]
[0047] 其中，
^待重建图像的对角线像素点个数，radbin表示索引坐标 为（i，j)的像素点在投影角度9 打在第几个像元上。
[0048] 3)将断层上的一点（i，j)与探测器像元radbin之间的投影映射关系存储到投影 矩阵System_Rad0中，S卩有：
[0050] 其中，气/H.为小于等于1的线段权值系数，其计算方法如式⑶所示。
[0051] 4)当遍历完所有像素点后，投影矩阵System_Rad0构造完成。
[0052] 投影矩阵SyStem_Rad0中每一行代表着一条Radon投影射线，投影方程中的一点 System_rad0 (radbin, (i-1)，N+j)意味着列矢量化后的第（i-1)，N+j个像素在第radbin 条投影射线中的贡献程度，在Radon投影模型中，一般将投影贡献权值等比例的视为穿过 该象素点的线段长值。
[0053] 6.Radon系统和Mojette系统相互转化的可行性。
[0054] 通过对比System_Rad0和System_Moj^的构造方式得到，当Radon投影系统中的 投影角度为0 =tarT^q/p)，探测器像元物理尺寸为
时，Radon投 影射线之间的距离和Mojette投影射线之间的距离完全一致，Radon投影系统与Mojette投影系统的对比如图3所示，这时Radon投影射线的投影路径与Mojette投影射线路径一 致，但与Mojette投影不同之处在于，Radon投影中还包括了穿过像素区域但没有穿过中心 点的像素值的贡献，因此将Radon投影看作是由一组权重系数不同的Mojette投影线性累 加构成，图3中的投影关系用代数形式来表达，有：
[0056] 其中，apq(i，j)表示Radon投影中的第i个分量和Mojette投影中第j个分量之 间的依存关系。
[0057] 但在实际的Radon投影中，探测器像元尺寸的值往往是固定的，并不能在每个投 影角度下都满足
，即Radon投 影射线之间的距离大于Mojette投影射线之间的距离，此时的Radon投影系统与Mojette投影系统的对比如图4所示，这时Radon投影射线的投影路径与Mojette投影射线路径不 一致，Radon投影更为稀疏，图4中的投影关系用代数形式来表达，有：
[0059]
，Radon投影系统与Mojette投影系统的对比如图5所 示。这时Radon投影射线的投影路径与Mojette投影射线路径不一致，Radon投影更为紧 凑，虽然此时的Radon投影射线并不一定穿过每一组构成Mojette投影中像素点集的中心， 但在这种更紧致条件下的Radon投影能穿过构成全部Mojette投影中的像素点集。图5中 的投影关系用代数形式来表达，有：
[0061] 设两投影之间转换矩阵为Ap，qGRTXB，使得RFT0 =Ap，q ?Mojp，q，其 中，J=DetRowNum代表Radon投影射线条数，也代表了探测器像元数量，B= (N-l)|p| + (M_l)|q|+l代表了需要转换的Mojette投影总个数。该矩阵中的一个分量 aRq(i，j)表示Radon投影中的第i个分量和Mojette投影中第j个分量之间的关系， 即若ap，q(i，j)辛0,则说明RFT0⑴投影值中有Mojp，q(j)中像素和的贡献；该矩阵中 的一行表示该投影角度下Radon投影的第i个分量与所有Mojette投影之间的关系，即
，若aP，q(i，j'）# 〇,则说明RFT0⑴投影值中有 MoUj'）中像素和的贡献。
[0062] 由RFT0 =A?Mojp^可知，当转换矩阵A可逆时，显然Radon投影是能够直接 转换为Mojette投影的，但这并不意味着只有当转换矩阵八^可逆时，Radon投影才能转化 为Mojette投影。即lAp，」乒0是Radon投影转化为Mojette投影的充分条件，却不是转 换的必要条件。
[0063] 例如，在探测器分辨率较小的情况下，S卩
时，构成同一个 Mojette投影的像素点集内穿过了多条Radon射线，例如RFT0 (1)zMoUl) ?ap，q(l, 1)，RFTJ2) =Mojp，q(l) ?a'pq(l，l)。这使得转换矩阵Ap，q中有多行线性相关，但由于存在 用来转化为全部Mojette投影的完备Radon投影集，即存在满足下列条件的Radon投影： 需要穿过构成Mojette投影的像素点集，且至少需要穿过两个相邻Mojette投影路径上的 像素点的并集，第一条投影射线除外，如，RFT0(2) =Mojp，q(2) ?apq(2,2)+M〇jp，q(l) ?ap ，(1(2，1)，贝1」1^1'0(2)穿过了构成舭知忧6投影舭叉，(1(2)肩〇叉，(1(1)中的像素点集，且穿过了 2个彼此相邻的Mojette投影路径上的像素点的并集，当Mojp，q(l)由RFT0 (1)求出后，根据 RFT0 (2)的值求出Mojp，q (2)。故而此处依然能够通过对线性方程的化简，将Radon投影转化 为Mojette投影。即，当转换矩阵Ap，q经过线性变换化简至最简形式时，去掉其中的非零行， 若其剩余代数余子式的行列式不为零，则该投影模型下的Radon投影能够转换为Mojette 投影。
[0064] 用数学条件来更准确地描述的话，即转换矩阵Ap，q中线性不相关的行矢量构成的 满秩矩阵的主对角线元素非零，或者说，当转换矩阵需要满足rank(A)彡（M-l) |Pi| + (N-l) qi| +1时，Radon投影满足转换为Mojette投影的完备性要求。
[0065] 若探测器分辨率较大，即
时，则采集的Radon投影 较稀疏，根本不能满足上述的转换条件。故而，在投影矢量（P，q)的投影方向下，当
、时，对应投影角度0itar^q/p，下的Radon投影都满足转换完 备性的要求。
[0066] 二、技术方案的步骤是：
[0067] S1设置投影参数，即探测器分辨率和重建图像分辨率；设探测器分辨率为 BXB，重建图像分辨率为MXN，其中，重建图像分辨率和投影矢量（p，q)的选择需要满足 max{(N_l) |pi| + (M-l) |qi|+l}〈B〇
[0068] S2选择投影矢量（p,q)，计算出相应的投影角度0 =tar^q/p，并在选定的投影 角度下采集Radon投影。
[0069] S3在投影矢量（p,q)下建立相应的0阶B样条插值投影系统System_rad0，需要 明确两点，其一为穿过每个像素点的投影射线的投影位置，其二为穿过每个像素点的投影 射线的线段权值。
[0070] S3. 1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影 射线打在探测器上的有效范围。
[0071] 设重建像素的实际物理尺寸为ObjPixel，探测器像元的实际物理尺寸为 DetCCDSize，当DetCCDSiZe〈ObjPiXel时，说明在一个矩形区域形状的离散像素点范围内， 有不止一条射线穿过该像素点。遍历每一个像素（i，j)，其中有l〈i〈M，l〈j〈N，并计算出穿 过该像素点的射线打在探测器上的绝对物理偏移值，此处用Xr_up和Xr_d〇wn来标记这个 有效范围的最小值和最大值，并计算出与Xr_up相对应的探测器像元标号最小值bin_up 和与Xr_down相对应的探测器像元标号最大值bin_down，投影射线打在探测器上的有效范 围的计算公式如式（12)所示：
[0076] 其中，（(i-l-M/2) ?ObjPixel，（j-l-N/2) ?ObjPixel)为离散化后的最小单位块 的左上角那一临界点的空间坐标，（(i-M/2) .ObjPixel，（j-N/2) .ObjPixel)为单位块的 右下角那一临界点的空间坐标，（(i-M/2) ?ObjPixel，（j-l-N/2) ?ObjPixel)为单位块的 左下角那一临界点的空间坐标，（(i-l-M/2) ?ObjPixel，（j-N/2) ?ObjPixel)为单位块的 右上角那一临界点的空间坐标，
为待重建图像的对角线像素点个数，如图6 所示。
[0077] 穿过像素中心点的射线打在探测器上的具体位置的计算公式，如式（13)所示。
[0078] xr= [((i-l/2)-M/2) ?cos( 0 ) + ((j-l/2)-N/2) *sin(0)] ?ObjPixel
[0079] (13)
[0080] bin= (xr+D/2 ?ObjPixel)/DetCCDSize+1
[0081] 其中，
为待重建图像的对角线像素点个数。
[0082] S3. 2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长 度。设在所有投影角度下穿过像素中心点的射线的线段权值为1，而对穿过该像素其它位置 的的射线，则需要根据这些射线与中心射线之间的垂直距离IAXr|，计算出其线段权重值 apq，其具体的计算公式如式（14)所示。
[0084] 其中，ObjPixel?( |cos9 |-|sin9 | )/2为一个临界距离值，当没有穿过中心点的 射线与穿过中心点的射线之间的垂直距离 IAXr|小于等于这一个临界值时，穿过像素点 的线段长度是最长的，将其记为1;与此类似，〇1^1^狀1*(|(3 〇80| + |81110|)/2也是一个临 界距离值，它标记着当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离IAXr 大于等于这一个临界值时，穿过像素点的线段长度是该投影角度下最短，已不在该像素点 的范围之内，将其记为0。
[0085] 并将计算出的权重值存入到系统矩阵System_rad0中相应的位置bin处，即
，其中N代表图像的列数。System_rad0 中每一行代表一条具体的投影射线穿过图像中所有点的线段权值矢量，每一列代表所有投 影射线穿过图像中一点的线段权值矢量。
[0086] S4根据已求得的系统矩阵System_rad0，来进行投影转换。
[0087] 此处需要注意，由于一个像素里通常都有超过1条射线穿过，具体的射线条数取 决于P=ObjPixel/DetCXDSize，即像素物理大小和探测器实际大小的比值，故而对于 Mojette采样来说，Radon采样就存在大量冗余，例如，仅穿过一个像素的Radon采样射线就 不止一条，而此处只需取其中一条进行计算，即能够求出对应的Mojette投影。
[0088] 因为系统方程SyStem_rad0的一行代表该投影角度下的一条投影射线，选择射 线其实就是在选择系统方程中的行矢量，选择这样的射线需要满足以下两个条件：1.该行 矢量非零；2.该行矢量中存在未求出的像素值，即该行矢量包含已求出的像素之外的新像 素，从而避免无用的重复计算。
[0089] S4. 1导入由S3中生成的系统矩阵System_rad0，找出系统方程中第一个非零行， 且该行矢量中只有一个非零分量，即意味着该射线只穿过了一个像素，投射在探测器上，将 探测器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Mojp，q(l)利用 下式计算：
[0090] Mojp；q(l) =RFTe(l)/ap；Q(l, 1)
[0091] 其中，RFT0 (1)表示Mojp，q(l)的对应位置的Radon投影值，apq(l，1)表示投影值 为RFT0 (1)的射线穿过的对应像素的线段权重值。
[0092] S4. 2寻找用来求出下一个Mojette投影的Radon投影射线。在S4. 1已求出的 Mojette投影MoUl)的基础上，若想求出下一个Mojette投影Mojp^O，需要找到形如RFTJ2) =Mojp，q(2) ?ap，q(2,2)+M〇jp，q(l) ?ap，q(2，l)的Radon投影，即，由且仅由已求 得的Mojette投影MoUl)和待求Mojette投影MoU〗）构成的Radon投影，若所有的 Radon投影都不满足这一条件，则该待求的Mojette投影Mo(2)无法被求出。
[0093] S4. 2. 1首先，找到待求Mojette投影Mojp，q(2)中穿过的一个像素点的索引坐标 (i，j)，并计算出其在列矢量化后的向量x中的列标1 = (i-1) ?N+j。
[0094] S4. 2. 2找到系统矩阵System_rad0的第1列，遍历System_rad0 (:，1)，找出其中 非零值对应的行数范围，这相当于遍历所有的投影射线，并找出穿过（i，j)这一点，权重系 数值不为零的投影射线，设这些投影矢量在投影矩阵中的行标最小值为bin_up，最大值为 bin_down。在这一范围内，遍历每一条投影射线在投影矩阵System_rad0*对应的横行， 即在bin_up<bin<bin_down范围内遍历System_rad0 (bin,:)，若系统矩阵行System_ rad0 (bin,:)中的非零系数处对应的像素点集，由且仅由Mojp，q(l)和Mojp，q(2)中穿过的像 素点构成，即RFT0(2) =Mojp，q(2) ?ap，q(2,2)+M〇jp，q(l) ? % (1(2，1)，则该投影射线即为 候选的更新射线。
[0095] S4. 2. 3在满足这一要求的射线当中，选择apq(2, 2)值最大的Radon投影射线作 为下一个用来更新Mojette投影的射线。
[0096] S4. 3利用求得的Mojette投影来更新Radon投影，S卩在该Radon投影中减去已求 得的Mojette投影，再除以线段权值，就得到了新的Mojette投影以〇1，(1(2) = (1^1'0(2)-M〇jP，q(l) ?ap，q(2，l))/ap，q(2,2)〇
[0097] 重复这一步骤，直到将所有的Mojette投影求出。
[0098] 即以式（15)为标准，寻找用来更新Mojette投影的Radon投影：
[0101] 基于已经求出的Mojette投影，从式（15)中推导出各Mojette投影的转换公式， 即如式（16)所示：
[0103] 至此，已将投影角度0 =tarriq/p下的Radon投影转化为投影矢量（p，q)下的 Mojette投影。
[0104] 本发明的效果和益处是在不同的稀疏投影角度下，不需要改变Radon投影图像的 分辨率大小，克服了Mojette投影域分辨率随着投影矢量变化和Radon投影分辨率不变之 间的转换障碍，搭建起这两种投影之间的转换桥梁。
【附图说明】
[0105] 图1是待重建物体所处的连续域空间坐标图。图中，1为固定坐标系，2为旋转坐 标系。重建物体所在的二维固定坐标系用x-o-y表示，坐标原点〇为物体的几何中心。模 拟光源与探测器在x-o-y平面内逆时针旋转，其旋转坐标系用\-o_h表示。
[0106] 图2是离散域索引坐标图，即将待重建物体所处的区域划分成等间距的正方行像 素块，离散化后的重建区域中每一个像素点的行列位置通过索引坐标（i，j)来标记，i代表 行数，j代表列数。
[0107] 图3是当探测器分辨率
1寸，Radon投影与Mojette投影的 关系示意图，图中，（a)为Radon投影过程；(b)为Mojette投影过程。
[0108] 图4是当探测器分辨率
时，Radon投影与Mojette投影的 关系示意图，图中，（a)为Radon投影过程；(b)为Mojette投影过程。
[0109] 图5是当探测器分辨率
_时，Radon投影与Mojette投影 的关系示意图，图中，（a)为Radon投影过程；(b)为Mojette投影过程。
[0110] 图6是穿过一个像素点的全部Radon投影射线示意图，图中，（a)为投影角度小于 90°时，穿过该像素点的射线打在探测器上的有效范围。图中，（b)为投影角度大于90° 时，穿过该像素点的射线打在探测器上的有效范围。
【具体实施方式】
[0111] 下面结合技术方案和附图详细说明本发明的【具体实施方式】：
[0112] S1设置投影参数，设在投影矢量（1，2)下，对一个3X3的小块进行投影，设探测 器像元个数B= 10> (3-1) 111 + (3-1) | 2 | +1 = 7 ;
[0113] S2投影矢量为（1，2)时，相应的投影角度0 =tan-il/2 = 26. 5651° ;
[0114] S3在投影矢量（1，2)下建立相应的0阶B样条插值投影系统SyStem_rad26.6，此处 设DetCCDSize= 0? 4mm。
[0115] S3. 1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影 射线打在探测器上的有效范围。
[0116] 设重建像素的实际物理尺寸为ObjPixel，ObjPixel= 1mm;探测器像元的实际物 理尺寸为DetCCDSize，DetCCDSize= 0? 4mm〇
[0117] 遍历每一个像素{Pixel(i,j) |i,jGZ,l〈i〈3,l〈j〈3}，并计算出穿过每一个像素 点的射线打在探测器上的有效范围，根据式（6)计算得出，穿过（1，1)该点的探测器像元有 效范围为1~4,即在（1，2)的投影矢量下，穿过第一个像素小块的各射线分别打在第1~ 4个探测器像元上；穿过（1，2)该点的射线分别打在第4~6个探测器像元上；穿过（1，3) 该点的射线分别打在第6~8个探测器像元上；穿过（2, 1)该点的射线分别打在第2~5 个探测器像元上，依次类推。
[0118] S3. 2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长 度。根据式（7)计算得出各权重的值，并存储至Systenuradi#应的位置处。
[0119] 经过计算得到的系统矩阵如下所示：
[0121] S4根据已求得的系统矩阵SyStem_rad0，来进行投影转换。
[0122] S4. 1导入由S3中生成的系统矩阵System_rad0，找出系统方程中第一个非零行， 且该行矢量中只有一个非零分量，意味着该射线只穿过了一个像素，投射在探测器上，将探 测器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Mojp，q(l)利用下 式计算：
[0123] Mojp；q(l) =RFT0 (l)/〇. 2038
[0124] S4. 2寻找用来求出下一个Mojette投影的Radon投影射线。
[0125] S4. 2. 1首先，找到构成Mojp，q(2)的像素点的索引坐标位置，此处为（2, 1)，并计算 出其在列矢量化后的向量中的位置1 = (2-1) ? 3+1 = 4。
[0126] S4. 2. 2找到系统矩阵System_rad0的第4列，遍历System_rad0 (:，4)，找出其中 非零值所在的行数范围，设最小值为bin_up，最大值为bin_down，此处，bin_up= 2,bin_ down= 5,在2 <bin< 5范围内遍历System_rad0 (bin,:)，观察其投影路径上穿过的像 素点的位置。若该系统矩阵行中的非零系数对应位置处的像素点，由且仅由Mojette投影 射线Mojp，q(l)和Mojp，q(2)中穿过的像素点构成，此处，点（1，1)(在矩阵中第1列）构成 Mojp，q(l)，点（2, 1)(在矩阵中第4列）构成Mojp，q(2)，若有RFT0 (bin)zMoUl) ?a12(b in,D+Mojy⑵?auODin, 2)，贝该投影射线System_rad0 (bin,:)即为候选的更新射线。 此处投影射线System_rad0 (2,：）和System_rad0 (3,：）均满足要求 ，而System_rad0 (4,：） 和System_rad0 (5,：）中还包含了Mojette投影射线MOU3)中穿过的像素点，因此不能 作为求出Mojette投影Mojp^(2)的候选Radon射线。
[0127] S4. 2. 3在满足这一要求的射线当中，选择a12(bin，2)值最大的Radon投影射 线作为下一个用来更新Mojette投影的射线，即选择System_rad0 (3,：）这条射线得到的 Radon投影来进行求解，此时的a12 (3, 2) =0? 9927。
[0128] S4. 3利用求得的Mojette投影来更新Radon投影，S卩在该Radon投影中减去已求 得的Mojette投影，再除以线段权值，就得到了新的Mojette投影Mojp，q(2) = (RFT0 (3)_M ojp；q(l) ? 1)/0.9927〇
[0129] 重复这一步骤，直到将所有的Mojette投影求出。
[0130] 此处的转换矩阵\2如下所示：
[0132] 各Mojette投影的求解过程如下：
[0134] 至此，已将投影角度0 =tar^q/p下的Radon投影转化为投影矢量（p,q)下的Mojette投影。
【主权项】
1.基于固定分辨率条件下的离散Radon和Mojette投影转换方法，以计算机层析成像 算法为基础，通过构造合理的成像系统来逼近Radon和Mo jette投影场景，并在详细分析了 Mojette投影和Radon投影的关系之后，给出了在固定分辨率前提下各投影角度下的Radon 投影转化为对应的投影矢量下的Mojette投影的充要条件和具体的转化算法，其特征是， 步骤如下： Sl设置投影参数，即探测器分辨率和重建图像分辨率；设探测器分辨率为BXB， 重建图像分辨率为MXN，其中，重建图像分辨率和投影矢量（p，q)的选择需要满足 max {(N-I) I Pi I + (M-I) I qi I +1}〈B ; S2选择投影矢量（p，q)，计算出相应的投影角度Θ itar^q/p，并在选定的投影角度 下采集Radon投影； S3在投影矢量（p, q)下建立相应的0阶B样条插值投影系统System_rad0，需要明确 两点，其一为穿过每个像素点的投影射线的投影位置，其二为穿过每个像素点的投影射线 的线段权值； S3. 1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影射线 打在探测器上的有效范围； 设重建像素的实际物理尺寸为ObjPixel，探测器像元的实际物理尺寸为DetCXDSize， 当DetCXDSizKObjPixel时，说明在一个矩形区域形状的离散像素点范围内，有不止一条 射线穿过该像素点；遍历每一个像素 （i, j)，其中有l〈i〈M, l〈j〈N，并计算出穿过该像素点 的射线打在探测器上的绝对物理偏移值，此处用Xr_up和Xr_down来标记这个有效范围的 最小值和最大值，并计算出与Xr_up相对应的探测器像元标号最小值bin_up和与Xr_down 相对应的探测器像元标号最大值bin_down，投影射线打在探测器上的有效范围的计算公式 如式（12)所示：其中，（(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)所示： xr= [ ((i-1/2)-M/2) · cos ( Θ ) + ((j-l/2)-N/2) ?sin(0)] · ObjPixel (13) bin = (xr+D/2 · ObjPixel)/DetCCDSize+1 其中，为待重建图像的对角线像素点个数； S3. 2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长度；设 在所有投影角度下穿过像素中心点的射线的线段权值为1，而对穿过该像素其它位置的的 射线，则需要根据这些射线与中心射线之间的垂直距离I ΔΧΓ|，计算出其线段权重值αρ(?， 其具体的计算公式如式（14)所示：其中，ObjPixel · (|cos Θ |-|sin Θ |)/2为一个临界距离值，当没有穿过中心点的射线 与穿过中心点的射线之间的垂直距离I ΔΧι·|小于等于这一个临界值时，穿过像素点的线 段长度是最长的，将其记为1;与此类似，ObjPixel · (Icos Θ | + | sin Θ |)/2也是一个临界距 离值，它标记着当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离I ΔΧΓ|大 于等于这一个临界值时，穿过像素点的线段长度是该投影角度下最短，已不在该像素点的 范围之内，将其记为〇; 并将计算出的权重值存入到系统矩阵System_rad 0中相应的位置bin处，即，其中N代表图像的列数；System_rad θ 中每一行代表一条具体的投影射线穿过图像中所有点的线段权值矢量，每一列代表所有投 影射线穿过图像中一点的线段权值矢量； S4根据已求得的系统矩阵Systenurad0，来进行投影转换； 此处需要注意，由于一个像素里通常都有超过1条射线穿过，具体的射线条数取决于 P = ObjPixel/DetCCDSize，即像素物理大小和探测器实际大小的比值，故而对于Mojette 采样来说，Radon采样就存在大量冗余，例如，仅穿过一个像素的Radon采样射线就不止一 条，而此处只需取其中一条进行计算，即求出了对应的Mojette投影； 因为系统方程Systenurad0的一行代表该投影角度下的一条投影射线，选择射线其实 就是在选择系统方程中的行矢量，选择这样的射线需要满足以下两个条件：1.该行矢量非 零；2.该行矢量中存在未求出的像素值，即该行矢量包含已求出的像素之外的新像素，从 而避免无用的重复计算； S4. 1导入由S3中生成的系统矩阵System_rad0，找出系统方程中第一个非零行，且该 行矢量中只有一个非零分量，即意味着该射线只穿过了一个像素，投射在探测器上，将探测 器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Moj^(I)利用下式 计算： Mojp；q(l) = RFT0 (l)/ap；q(l, 1) 其中，RFT0(I)表示Moj1^(I)的对应位置的Radon投影值，a M(l，l)表示投影值为 RFT0 (1)的射线穿过的对应像素的线段权重值； S4. 2寻找用来求出下一个Mojette投影的Radon投影射线；在S4. 1已求出的Mojette 投影Mojp，q(l)的基础上，若想求出下一个Mojette投影Mojp， q (2)，需要找到形如RFT0 (2)= Mojp，q(2) · · Ctp q(2, 1)的 Radon 投影，即，由且仅由已求得的 Mojette 投影Mojp，q(1)和待求Mojette投影Mo jp，q(2)构成的Radon投影，若所有的Radon投影都不 满足这一条件，则该待求的Mojette投影Mo jp，q(2)无法被求出； S4. 2. 1首先，找到待求Mojette投影Mojp，q(2)中穿过的一个像素点的索引坐标（i, j)， 并计算出其在列矢量化后的向量X中的列标I = (i-1) · N+j ; S4. 2. 2找到系统矩阵System_rad0的第1列，遍历System_rad 0 (:，I)，找出其中非 零值对应的行数范围，这相当于遍历所有的投影射线，并找出穿过（i，j)这一点，权重系 数值不为零的投影射线，设这些投影矢量在投影矩阵中的行标最小值为bin_up，最大值为 bin_down ;在这一范围内，遍历每一条投影射线在投影矩阵System_rad0中对应的横行， 即在 bin_up < bin < bin_down 范围内遍历 System_rad0 (bin,:)，若系统矩阵行 System_ rad0 (bin,:)中的非零系数处对应的像素点集，由且仅由Mojm(I)和Mojp，q(2)中穿过的像 素点构成，即1^1'0(2)=厘〇几(1(2).€%， (1(2,2)+厘〇夂(1(1).€%，(1(2，1)，则该投影射线即为 候选的更新射线； S4. 2. 3在满足这一要求的射线当中，选择α ^(2, 2)值最大的Radon投影射线作为下 一个用来更新Mojette投影的射线； S4. 3利用求得的Mojette投影来更新Radon投影，即在该Radon投影中减去已求得的 Mojette投影，再除以线段权值，就得到了新的Mojette投影： Mojp；q(2) = (RFT0 (2)-Mojp；q(l) * a p；q(2, I))/a p；q(2, 2)； 重复这一步骤，直到将所有的Mojette投影求出； 即以式（15)为标准，寻找用来更新Mojette投影的Radon投影：其中，j基于已经求出的Mojette投影，从式（15)中推导出各Mojette投影的转换公式，即如 式（16)所示：至此，已将投影角度Θ = tar^q/p下的Radon投影转化为投影矢量（p, q)下的Mojette 投影。
【专利摘要】基于固定分辨率条件下的离散Radon和Mojette投影转换方法，属于图像处理领域的计算机成像技术领域。其特征是将实际成像条件下的稀疏角度离散Radon投影转化为在离散域精确重建的Mojette投影。以计算机层析成像算法为基础，通过构造合理的成像系统来逼近Radon和Mojette投影场景，并在详细分析了Mojette投影和Radon投影的关系之后，给出了在固定分辨率前提下各投影角度下的Radon投影转化为对应的投影矢量下的Mojette投影的具体算法。本发明的效果和益处是在不同的稀疏投影角度下，不需要改变Radon投影图像的分辨率大小，克服了Mojette投影域分辨率随着投影矢量变化和Radon投影分辨率不变之间的转换障碍，搭建起这两种投影之间的转换桥梁。
【IPC分类】G06T3/20
【公开号】CN104899827
【申请号】CN201510274374
【发明人】孙怡, 李梦婕, 蒋敏
【申请人】大连理工大学
【公开日】2015年9月9日
【申请日】2015年5月26日
【技术领域】
[0001] 本发明属于图像处理领域里的计算机层析成像（ComputedTomography，CT)技术 领域，涉及在图像投影变换域的稀疏采样及高效复原问题。
【背景技术】
[0002] CT技术是一种基于计算机的三维重构技术，能够在不损坏物体表层和内部构造的 前提下，重现出物体不可见的内部切片结构。在CT成像算法上，按照成像原理的不同，分为 解析类和迭代类重建算法；按照重建精确性的不同，又分为近似重建和精确重建算法。在解 析算法中，奥地利数学学家Radon提出的变换模型利用低维的压缩投影数据来重构高维度 的物体结构，这一数学模型的提出对CT的图像重建技术起到了关键性的启蒙和指导作用。 在这一理论基础的指引下，滤波反投影等解析算法相继被提出，该算法在完备的投影角度 下，能精确的重建出图像断层，且能抑制在重建过程中产生的各种伪影，在CT成像领域具 有里程碑式的意义。
[0003] 但由于成像系统是数字离散系统，因而连续的投影模型无法适用于离散的图像， 故而离散Radon变换随机应运而生。离散Radon变换及其逆变换解决了从模拟域向数字域 转换的问题，但基于Radon逆变换的重建算法需要大量的采样样本和投影数量，才能重建 出较好的重建断层，这种重建需求必然会导致过高的辐射剂量和过长的重建时间，出于对 医学成像中病人身体健康和病灶动态变化会导致伪影的考虑，在短时间、低剂量和稀疏角 度下进行高质量成像就显得尤为重要，同时又充满挑战。
[0004] 在Radon变换的基础上，GuedonJP等提出了Mojette变换的概念，该变换是离 散Radon变换的一种特殊形式，Mojette变换通过改变不同投影矢量下的采样率，能够在最 大程度上避免像素点重复和冗余采样，以避免出现一部分像素点过采样，而另一部分像素 点欠采样而带来的重建效果上的不精确。并通过充分利用已重建出来的像素点的信息，从 而大大减小了重建断层所需的投影角度和投影射线条数。
[0005] 但由于Mojette投影系统要求探测器分辨率随着投影角度的变化而变化，且在投 影角度之间的步进角也不是一个固定的常值，这与现行的采样模式大相径庭，实际的投影 系统符合Radon变换的采样模式，探测器分辨率不随投影角度而变化，因此通过Radon变换 得到的投影值不可能直接参与到Mojette逆变换当中来，而是先要将Radon投影转换到相 关的Mojette投影上去，再利用Mojette逆变换来进行反重建。
[0006] 完成从Radon投影到Mojette投影的转换，其中需要解决几个实际的问题：1.投 影矢量和投影角度之间的匹配转换，由于Mojette投影系统中采用一对互质的整数来表达 投影方向，称之为投影矢量，而Radon投影系统中采用的是极坐标中的旋转角度，故而需要 通过简单的换算将Mojette投影矢量表达的方向转换成Radon投影中的旋转角；2.Mojette 变换中可变分辨率与Radon变换中固定分辨率之间的矛盾，由于Mojette变换中的探测器 分辨率或者说是投影射线条数为随着投影矢量而变化的函数，而Radon投影的探测器分辨 率是不随着投影角度而变化的，这就使得采样系统中得到的Radon投影需要通过本发明中 提出的算法进行转换，才能得到可以高效重建图像的Mojette投影。
【发明内容】
[0007] 本发明的目的是提供了一种基于固定分辨率采集的Radon投影转化为Mojette 投影的算法，该方法在详细分析了Mojette投影和Radon投影的关系之后，解决了固定分辨 率条件下的两投影之间的转化难题，给出了Radon投影转化为对应投影矢量下的Mojette 投影的充要条件，并基于此给出了具体的转化算法，使得实际成像条件下采集的Radon投 影能够转化为Mojette域的投影。
[0008] 本发明的技术方案的原理及步骤如下：
[0009] -、技术方案的原理是：
[0010] 1.建立投影空间坐标系：设投影空间坐标系用x-o-y表示，如图1所示，坐标原 点〇为待重建物体的几何中心，X轴正向与探测器像元索引号增加方向一致，y轴正向与 射线传播方向一致，x-o-y平面内的一点的空间坐标为（x，y)。光源与探测器在x-o-y平 面内逆时针旋转，旋转角度为0，设旋转后的空间坐标系用\-o_h表示，像素点的旋转坐 标为在连续域的二维图像空间内，将物体所在区域离散化为MXN的离散小块， 设ObjSize为待重建矩形区域的行方向上的边长，则每一个重建像素的物理尺寸大小为 ObjPixel=ObjSiZe/M，如图2所示。其中，（i，j)为标记离散小块行列位置的索引坐标， 即若利用二维矩阵存储该离散图像，则i代表在矩阵中的行数，j代表在矩阵中的列数，通 常将图像左上角第一个小块的位置记为索引坐标起始位置（1，1)。
[0011] 明确重建断层像素在x-〇-y平面内的空间坐标（X，y)和在旋转坐标系Xj-o-yj?平 面内旋转坐标0^，^)之间的关系，如式（1)所示。
[0013] 明确重建断层像素在x-0-y平面内的空间坐标（x，y)和离散域索引坐标（i，j)之 间的关系，通常，以待重建断层的左上角像素为起点，记该像素的索引坐标为（1，1)，索引坐 标为（1，1)的像素中心点的空间坐标（x，y)为：
[0014] ((l-(M+l)/2) ?ObjPixel，（l-(N+l)/2) ?ObjPixel)
[0015] 索引坐标为（i，j)的像素中心点的空间坐标（x，y)为：
[0016] ((i-(M+l)/2) ?ObjPixel, (j-(N+l)/2) ?ObjPixel)
[0017] 基于以上两种坐标的对应，则得到索引坐标为（i，j)的像素中心点在旋转坐标系 中的旋转坐标为：
[0018] xr= [(j-(N+l)/2) ?ObjPixel?sin( 0 ) + (i-(M+l)/2) ?ObjPixel?cos( 0 )]
[0019] yr= [-(i-(M+l)/2) ?ObjPixel?sin( 0 ) + (j-(N+l)/2) ?ObjPixel?cos( 0 )]
[0020] 2.离散Radon变换：对分辨率大小为MXN的离散图像f进行离散Radon变换，其 过程用式（2)来表述：
[0022] 其中，1^1\表示在投影角度0下的Radon投影，RFT0 〇〇表示在投影角度0下 打在探测器x,位置处的Radon投影值。f(i，j)代表待重建的图像切片上索引坐标为（i，j) 的一点的灰度值。若以（0，\)标记投影角度0下打在探测器上的物理偏移量为&的一 条投影射线，％，从^表示在投影射线（0，\)和象素点（i，j)之间的权重核函数，通常为 S函数或0次样条函数（即将投影贡献权值等比例的视为穿过该象素点的线段长值）。
[0023] 设在所有投影角度0下穿过像素中心点的射线的线段权值为1，而对穿过该像素 其它位置的的射线，则需要根据这些射线与中心射线之间的垂直距离IAX」，计算出其线 段权重值％，从，&，其具体的计算公式如式（3)所示。
[0025] 其中，ObjPixel代表每一个重建像素的物理尺寸。
[0026] 3.Mojette变换：Mojette变换中用一对互质的整数（p,q)来表达投影方向，一般 有pGZ，qeZ+，p代表了图像列方向上的整数位移，q代表了图像行方向上的整数位移， 投影矢量（P，q)表达的投影角度0 =tarTHq/p)，对分辨率大小为MXN的离散图像f进 行Mojette变换，其过程用式（4)来表述：
[0028] 其中，1\1〇几(1表示在投影矢量（p,q)下的Mojette投影，Moj^(bin)表示在投影矢 量（P，q)下打在探测器bin上Mojette投影值。f(i，j)代表待重建的图像切片上索引坐标 为（i，j)的一点的灰度值，P〇为一个探测器像元位置矫正量，当投影矢量P>〇时，P〇=〇; 当投影矢量P〈〇时，P〇= (N-1) ?p。
[0029] Mojette变换与Radon变换相比最大的不同在于，Radon变换中探测器分辨率为 固定的常值，穿过物体切片的Radon投影射线之间距离是固定的，因而不能保证每条Radon 投影射线都能穿过像素的中心点，穿过像素中心点的线段权值记为1，而不穿过像素中心点 的线段权值要根据具体线段长度而定，计算方法如式（3)所示；而Mojette变换中每条投影 射线都只穿过像素的中心点，因此不同投影矢量（P，q)下探测器分辨率B(M，N，p，q)不同， 探测器分辨率的值是由图像重建分辨率MXN和投影矢量（p，q)共同决定的，即投影矢量 (p,q)下的Mojette探测器分辨率为B(P,Q,p,q) = (Q-1) |p| + (P-1) |q|+1。这意味着在覆 盖相同直径范围内的切片时，不同投影角度下的Mojette投影射线之间的间距hi不同：
[0031] 此外，Radon变换在模拟域重建是精确的，在离散域重建是近似的。而Mojette变 换在离散域的重建是精确的。
[0032] 4.Mojette投影矩阵的构造方法：
[0033] 基于Mojette正变换原理，下面阐述一下Mojette投影矩阵的构造方法。
[0034] 1)初始化参数，设重建断层分辨率为MXN，投影矢量为（p，q)，投影矩阵为
'，其中，B= (N-l)|Pi| + (M-l)|qi|+l，V=M.N。待重建图像列矢 量化后的一维向量为X，投影矢量（p,q)下的Mojette投影为Mojj^，则有Systen^Mojp^*x=M〇jp；q；
[0035] 2)在每个投影矢量（p，q)下，遍历每一个像素点（i，j)，计算其打在探测器上的位 置bin:
[0037] bin=p? (i_l)+q? (j_l)-po+1
[0038] 3)将通过投影矢量（Pi，qi)建立起来的断层上的一点（i，j)与探测器像元bin之 间的投影映射关系，存储到线性投影矩阵SyStem_M〇jp，q中，即：
[0039] System_Mojp；q(bin, (i~l) ?N+j) = 1 ；
[0040] 4)当遍历完所有像素点后，投影矩阵System_Mojp;(1构造完成。
[0041] 投影矩阵Systen^Moj^if每一行代表着一条Mojette投影射线，每一列 代表着穿过一个像素点的所有投影射线的线段权值。投影方程中的一 点SyStem_ Mojp，q(bin，（i-1) *N+j)意味着列矢量化后的第（i-1) *N+j个像素在第bin条投影射线中 的贡献程度，在Mojette投影模型中，若第bin条投影射线穿过第（i-1) *N+j个像素的中 心点，则System-MoUbin，（i-1) .N+j) = 1，否则为 0。
[0042] 5.Radon投影矩阵的构造方法：
[0043] 在Radon投影当中，线性投影矩阵的构造方法是由探测器和射线源与待重建物体 之间的几何关系决定的。
[0044] 1)设重建断层分辨率为MXN，包含待重建断层的最小正方形区域的边长为 ObjSize，每一个重建像素的实际物理尺寸ObjPixel=ObjSize/M。设探测器分辨率 为DetRowNum，其中每一个探测器像元的物理尺寸为DetCXDSize，投影矩阵为System_ rad0GRJ=DetRowNum，V=M*N。待重建图像列矢量化后的一维向量为x，投影角度 Q下的Radon投影为RFT0，贝丨」有System_rad0 ?x=RFT0;
[0045] 2)在投影角度0下，遍历每一个像素点（i，j)，计算其打在探测器上的绝对偏移 量Xj?和相对探测器像元位置radbin:
[0046]
[0047] 其中，
^待重建图像的对角线像素点个数，radbin表示索引坐标 为（i，j)的像素点在投影角度9 打在第几个像元上。
[0048] 3)将断层上的一点（i，j)与探测器像元radbin之间的投影映射关系存储到投影 矩阵System_Rad0中，S卩有：
[0050] 其中，气/H.为小于等于1的线段权值系数，其计算方法如式⑶所示。
[0051] 4)当遍历完所有像素点后，投影矩阵System_Rad0构造完成。
[0052] 投影矩阵SyStem_Rad0中每一行代表着一条Radon投影射线，投影方程中的一点 System_rad0 (radbin, (i-1)，N+j)意味着列矢量化后的第（i-1)，N+j个像素在第radbin 条投影射线中的贡献程度，在Radon投影模型中，一般将投影贡献权值等比例的视为穿过 该象素点的线段长值。
[0053] 6.Radon系统和Mojette系统相互转化的可行性。
[0054] 通过对比System_Rad0和System_Moj^的构造方式得到，当Radon投影系统中的 投影角度为0 =tarT^q/p)，探测器像元物理尺寸为
时，Radon投 影射线之间的距离和Mojette投影射线之间的距离完全一致，Radon投影系统与Mojette投影系统的对比如图3所示，这时Radon投影射线的投影路径与Mojette投影射线路径一 致，但与Mojette投影不同之处在于，Radon投影中还包括了穿过像素区域但没有穿过中心 点的像素值的贡献，因此将Radon投影看作是由一组权重系数不同的Mojette投影线性累 加构成，图3中的投影关系用代数形式来表达，有：
[0056] 其中，apq(i，j)表示Radon投影中的第i个分量和Mojette投影中第j个分量之 间的依存关系。
[0057] 但在实际的Radon投影中，探测器像元尺寸的值往往是固定的，并不能在每个投 影角度下都满足
，即Radon投 影射线之间的距离大于Mojette投影射线之间的距离，此时的Radon投影系统与Mojette投影系统的对比如图4所示，这时Radon投影射线的投影路径与Mojette投影射线路径不 一致，Radon投影更为稀疏，图4中的投影关系用代数形式来表达，有：
[0059]
，Radon投影系统与Mojette投影系统的对比如图5所 示。这时Radon投影射线的投影路径与Mojette投影射线路径不一致，Radon投影更为紧 凑，虽然此时的Radon投影射线并不一定穿过每一组构成Mojette投影中像素点集的中心， 但在这种更紧致条件下的Radon投影能穿过构成全部Mojette投影中的像素点集。图5中 的投影关系用代数形式来表达，有：
[0061] 设两投影之间转换矩阵为Ap，qGRTXB，使得RFT0 =Ap，q ?Mojp，q，其 中，J=DetRowNum代表Radon投影射线条数，也代表了探测器像元数量，B= (N-l)|p| + (M_l)|q|+l代表了需要转换的Mojette投影总个数。该矩阵中的一个分量 aRq(i，j)表示Radon投影中的第i个分量和Mojette投影中第j个分量之间的关系， 即若ap，q(i，j)辛0,则说明RFT0⑴投影值中有Mojp，q(j)中像素和的贡献；该矩阵中 的一行表示该投影角度下Radon投影的第i个分量与所有Mojette投影之间的关系，即
，若aP，q(i，j'）# 〇,则说明RFT0⑴投影值中有 MoUj'）中像素和的贡献。
[0062] 由RFT0 =A?Mojp^可知，当转换矩阵A可逆时，显然Radon投影是能够直接 转换为Mojette投影的，但这并不意味着只有当转换矩阵八^可逆时，Radon投影才能转化 为Mojette投影。即lAp，」乒0是Radon投影转化为Mojette投影的充分条件，却不是转 换的必要条件。
[0063] 例如，在探测器分辨率较小的情况下，S卩
时，构成同一个 Mojette投影的像素点集内穿过了多条Radon射线，例如RFT0 (1)zMoUl) ?ap，q(l, 1)，RFTJ2) =Mojp，q(l) ?a'pq(l，l)。这使得转换矩阵Ap，q中有多行线性相关，但由于存在 用来转化为全部Mojette投影的完备Radon投影集，即存在满足下列条件的Radon投影： 需要穿过构成Mojette投影的像素点集，且至少需要穿过两个相邻Mojette投影路径上的 像素点的并集，第一条投影射线除外，如，RFT0(2) =Mojp，q(2) ?apq(2,2)+M〇jp，q(l) ?ap ，(1(2，1)，贝1」1^1'0(2)穿过了构成舭知忧6投影舭叉，(1(2)肩〇叉，(1(1)中的像素点集，且穿过了 2个彼此相邻的Mojette投影路径上的像素点的并集，当Mojp，q(l)由RFT0 (1)求出后，根据 RFT0 (2)的值求出Mojp，q (2)。故而此处依然能够通过对线性方程的化简，将Radon投影转化 为Mojette投影。即，当转换矩阵Ap，q经过线性变换化简至最简形式时，去掉其中的非零行， 若其剩余代数余子式的行列式不为零，则该投影模型下的Radon投影能够转换为Mojette 投影。
[0064] 用数学条件来更准确地描述的话，即转换矩阵Ap，q中线性不相关的行矢量构成的 满秩矩阵的主对角线元素非零，或者说，当转换矩阵需要满足rank(A)彡（M-l) |Pi| + (N-l) qi| +1时，Radon投影满足转换为Mojette投影的完备性要求。
[0065] 若探测器分辨率较大，即
时，则采集的Radon投影 较稀疏，根本不能满足上述的转换条件。故而，在投影矢量（P，q)的投影方向下，当
、时，对应投影角度0itar^q/p，下的Radon投影都满足转换完 备性的要求。
[0066] 二、技术方案的步骤是：
[0067] S1设置投影参数，即探测器分辨率和重建图像分辨率；设探测器分辨率为 BXB，重建图像分辨率为MXN，其中，重建图像分辨率和投影矢量（p，q)的选择需要满足 max{(N_l) |pi| + (M-l) |qi|+l}〈B〇
[0068] S2选择投影矢量（p,q)，计算出相应的投影角度0 =tar^q/p，并在选定的投影 角度下采集Radon投影。
[0069] S3在投影矢量（p,q)下建立相应的0阶B样条插值投影系统System_rad0，需要 明确两点，其一为穿过每个像素点的投影射线的投影位置，其二为穿过每个像素点的投影 射线的线段权值。
[0070] S3. 1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影 射线打在探测器上的有效范围。
[0071] 设重建像素的实际物理尺寸为ObjPixel，探测器像元的实际物理尺寸为 DetCCDSize，当DetCCDSiZe〈ObjPiXel时，说明在一个矩形区域形状的离散像素点范围内， 有不止一条射线穿过该像素点。遍历每一个像素（i，j)，其中有l〈i〈M，l〈j〈N，并计算出穿 过该像素点的射线打在探测器上的绝对物理偏移值，此处用Xr_up和Xr_d〇wn来标记这个 有效范围的最小值和最大值，并计算出与Xr_up相对应的探测器像元标号最小值bin_up 和与Xr_down相对应的探测器像元标号最大值bin_down，投影射线打在探测器上的有效范 围的计算公式如式（12)所示：
[0076] 其中，（(i-l-M/2) ?ObjPixel，（j-l-N/2) ?ObjPixel)为离散化后的最小单位块 的左上角那一临界点的空间坐标，（(i-M/2) .ObjPixel，（j-N/2) .ObjPixel)为单位块的 右下角那一临界点的空间坐标，（(i-M/2) ?ObjPixel，（j-l-N/2) ?ObjPixel)为单位块的 左下角那一临界点的空间坐标，（(i-l-M/2) ?ObjPixel，（j-N/2) ?ObjPixel)为单位块的 右上角那一临界点的空间坐标，
为待重建图像的对角线像素点个数，如图6 所示。
[0077] 穿过像素中心点的射线打在探测器上的具体位置的计算公式，如式（13)所示。
[0078] xr= [((i-l/2)-M/2) ?cos( 0 ) + ((j-l/2)-N/2) *sin(0)] ?ObjPixel
[0079] (13)
[0080] bin= (xr+D/2 ?ObjPixel)/DetCCDSize+1
[0081] 其中，
为待重建图像的对角线像素点个数。
[0082] S3. 2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长 度。设在所有投影角度下穿过像素中心点的射线的线段权值为1，而对穿过该像素其它位置 的的射线，则需要根据这些射线与中心射线之间的垂直距离IAXr|，计算出其线段权重值 apq，其具体的计算公式如式（14)所示。
[0084] 其中，ObjPixel?( |cos9 |-|sin9 | )/2为一个临界距离值，当没有穿过中心点的 射线与穿过中心点的射线之间的垂直距离 IAXr|小于等于这一个临界值时，穿过像素点 的线段长度是最长的，将其记为1;与此类似，〇1^1^狀1*(|(3 〇80| + |81110|)/2也是一个临 界距离值，它标记着当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离IAXr 大于等于这一个临界值时，穿过像素点的线段长度是该投影角度下最短，已不在该像素点 的范围之内，将其记为0。
[0085] 并将计算出的权重值存入到系统矩阵System_rad0中相应的位置bin处，即
，其中N代表图像的列数。System_rad0 中每一行代表一条具体的投影射线穿过图像中所有点的线段权值矢量，每一列代表所有投 影射线穿过图像中一点的线段权值矢量。
[0086] S4根据已求得的系统矩阵System_rad0，来进行投影转换。
[0087] 此处需要注意，由于一个像素里通常都有超过1条射线穿过，具体的射线条数取 决于P=ObjPixel/DetCXDSize，即像素物理大小和探测器实际大小的比值，故而对于 Mojette采样来说，Radon采样就存在大量冗余，例如，仅穿过一个像素的Radon采样射线就 不止一条，而此处只需取其中一条进行计算，即能够求出对应的Mojette投影。
[0088] 因为系统方程SyStem_rad0的一行代表该投影角度下的一条投影射线，选择射 线其实就是在选择系统方程中的行矢量，选择这样的射线需要满足以下两个条件：1.该行 矢量非零；2.该行矢量中存在未求出的像素值，即该行矢量包含已求出的像素之外的新像 素，从而避免无用的重复计算。
[0089] S4. 1导入由S3中生成的系统矩阵System_rad0，找出系统方程中第一个非零行， 且该行矢量中只有一个非零分量，即意味着该射线只穿过了一个像素，投射在探测器上，将 探测器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Mojp，q(l)利用 下式计算：
[0090] Mojp；q(l) =RFTe(l)/ap；Q(l, 1)
[0091] 其中，RFT0 (1)表示Mojp，q(l)的对应位置的Radon投影值，apq(l，1)表示投影值 为RFT0 (1)的射线穿过的对应像素的线段权重值。
[0092] S4. 2寻找用来求出下一个Mojette投影的Radon投影射线。在S4. 1已求出的 Mojette投影MoUl)的基础上，若想求出下一个Mojette投影Mojp^O，需要找到形如RFTJ2) =Mojp，q(2) ?ap，q(2,2)+M〇jp，q(l) ?ap，q(2，l)的Radon投影，即，由且仅由已求 得的Mojette投影MoUl)和待求Mojette投影MoU〗）构成的Radon投影，若所有的 Radon投影都不满足这一条件，则该待求的Mojette投影Mo(2)无法被求出。
[0093] S4. 2. 1首先，找到待求Mojette投影Mojp，q(2)中穿过的一个像素点的索引坐标 (i，j)，并计算出其在列矢量化后的向量x中的列标1 = (i-1) ?N+j。
[0094] S4. 2. 2找到系统矩阵System_rad0的第1列，遍历System_rad0 (:，1)，找出其中 非零值对应的行数范围，这相当于遍历所有的投影射线，并找出穿过（i，j)这一点，权重系 数值不为零的投影射线，设这些投影矢量在投影矩阵中的行标最小值为bin_up，最大值为 bin_down。在这一范围内，遍历每一条投影射线在投影矩阵System_rad0*对应的横行， 即在bin_up<bin<bin_down范围内遍历System_rad0 (bin,:)，若系统矩阵行System_ rad0 (bin,:)中的非零系数处对应的像素点集，由且仅由Mojp，q(l)和Mojp，q(2)中穿过的像 素点构成，即RFT0(2) =Mojp，q(2) ?ap，q(2,2)+M〇jp，q(l) ? % (1(2，1)，则该投影射线即为 候选的更新射线。
[0095] S4. 2. 3在满足这一要求的射线当中，选择apq(2, 2)值最大的Radon投影射线作 为下一个用来更新Mojette投影的射线。
[0096] S4. 3利用求得的Mojette投影来更新Radon投影，S卩在该Radon投影中减去已求 得的Mojette投影，再除以线段权值，就得到了新的Mojette投影以〇1，(1(2) = (1^1'0(2)-M〇jP，q(l) ?ap，q(2，l))/ap，q(2,2)〇
[0097] 重复这一步骤，直到将所有的Mojette投影求出。
[0098] 即以式（15)为标准，寻找用来更新Mojette投影的Radon投影：
[0101] 基于已经求出的Mojette投影，从式（15)中推导出各Mojette投影的转换公式， 即如式（16)所示：
[0103] 至此，已将投影角度0 =tarriq/p下的Radon投影转化为投影矢量（p，q)下的 Mojette投影。
[0104] 本发明的效果和益处是在不同的稀疏投影角度下，不需要改变Radon投影图像的 分辨率大小，克服了Mojette投影域分辨率随着投影矢量变化和Radon投影分辨率不变之 间的转换障碍，搭建起这两种投影之间的转换桥梁。
【附图说明】
[0105] 图1是待重建物体所处的连续域空间坐标图。图中，1为固定坐标系，2为旋转坐 标系。重建物体所在的二维固定坐标系用x-o-y表示，坐标原点〇为物体的几何中心。模 拟光源与探测器在x-o-y平面内逆时针旋转，其旋转坐标系用\-o_h表示。
[0106] 图2是离散域索引坐标图，即将待重建物体所处的区域划分成等间距的正方行像 素块，离散化后的重建区域中每一个像素点的行列位置通过索引坐标（i，j)来标记，i代表 行数，j代表列数。
[0107] 图3是当探测器分辨率
1寸，Radon投影与Mojette投影的 关系示意图，图中，（a)为Radon投影过程；(b)为Mojette投影过程。
[0108] 图4是当探测器分辨率
时，Radon投影与Mojette投影的 关系示意图，图中，（a)为Radon投影过程；(b)为Mojette投影过程。
[0109] 图5是当探测器分辨率
_时，Radon投影与Mojette投影 的关系示意图，图中，（a)为Radon投影过程；(b)为Mojette投影过程。
[0110] 图6是穿过一个像素点的全部Radon投影射线示意图，图中，（a)为投影角度小于 90°时，穿过该像素点的射线打在探测器上的有效范围。图中，（b)为投影角度大于90° 时，穿过该像素点的射线打在探测器上的有效范围。
【具体实施方式】
[0111] 下面结合技术方案和附图详细说明本发明的【具体实施方式】：
[0112] S1设置投影参数，设在投影矢量（1，2)下，对一个3X3的小块进行投影，设探测 器像元个数B= 10> (3-1) 111 + (3-1) | 2 | +1 = 7 ;
[0113] S2投影矢量为（1，2)时，相应的投影角度0 =tan-il/2 = 26. 5651° ;
[0114] S3在投影矢量（1，2)下建立相应的0阶B样条插值投影系统SyStem_rad26.6，此处 设DetCCDSize= 0? 4mm。
[0115] S3. 1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影 射线打在探测器上的有效范围。
[0116] 设重建像素的实际物理尺寸为ObjPixel，ObjPixel= 1mm;探测器像元的实际物 理尺寸为DetCCDSize，DetCCDSize= 0? 4mm〇
[0117] 遍历每一个像素{Pixel(i,j) |i,jGZ,l〈i〈3,l〈j〈3}，并计算出穿过每一个像素 点的射线打在探测器上的有效范围，根据式（6)计算得出，穿过（1，1)该点的探测器像元有 效范围为1~4,即在（1，2)的投影矢量下，穿过第一个像素小块的各射线分别打在第1~ 4个探测器像元上；穿过（1，2)该点的射线分别打在第4~6个探测器像元上；穿过（1，3) 该点的射线分别打在第6~8个探测器像元上；穿过（2, 1)该点的射线分别打在第2~5 个探测器像元上，依次类推。
[0118] S3. 2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长 度。根据式（7)计算得出各权重的值，并存储至Systenuradi#应的位置处。
[0119] 经过计算得到的系统矩阵如下所示：
[0121] S4根据已求得的系统矩阵SyStem_rad0，来进行投影转换。
[0122] S4. 1导入由S3中生成的系统矩阵System_rad0，找出系统方程中第一个非零行， 且该行矢量中只有一个非零分量，意味着该射线只穿过了一个像素，投射在探测器上，将探 测器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Mojp，q(l)利用下 式计算：
[0123] Mojp；q(l) =RFT0 (l)/〇. 2038
[0124] S4. 2寻找用来求出下一个Mojette投影的Radon投影射线。
[0125] S4. 2. 1首先，找到构成Mojp，q(2)的像素点的索引坐标位置，此处为（2, 1)，并计算 出其在列矢量化后的向量中的位置1 = (2-1) ? 3+1 = 4。
[0126] S4. 2. 2找到系统矩阵System_rad0的第4列，遍历System_rad0 (:，4)，找出其中 非零值所在的行数范围，设最小值为bin_up，最大值为bin_down，此处，bin_up= 2,bin_ down= 5,在2 <bin< 5范围内遍历System_rad0 (bin,:)，观察其投影路径上穿过的像 素点的位置。若该系统矩阵行中的非零系数对应位置处的像素点，由且仅由Mojette投影 射线Mojp，q(l)和Mojp，q(2)中穿过的像素点构成，此处，点（1，1)(在矩阵中第1列）构成 Mojp，q(l)，点（2, 1)(在矩阵中第4列）构成Mojp，q(2)，若有RFT0 (bin)zMoUl) ?a12(b in,D+Mojy⑵?auODin, 2)，贝该投影射线System_rad0 (bin,:)即为候选的更新射线。 此处投影射线System_rad0 (2,：）和System_rad0 (3,：）均满足要求 ，而System_rad0 (4,：） 和System_rad0 (5,：）中还包含了Mojette投影射线MOU3)中穿过的像素点，因此不能 作为求出Mojette投影Mojp^(2)的候选Radon射线。
[0127] S4. 2. 3在满足这一要求的射线当中，选择a12(bin，2)值最大的Radon投影射 线作为下一个用来更新Mojette投影的射线，即选择System_rad0 (3,：）这条射线得到的 Radon投影来进行求解，此时的a12 (3, 2) =0? 9927。
[0128] S4. 3利用求得的Mojette投影来更新Radon投影，S卩在该Radon投影中减去已求 得的Mojette投影，再除以线段权值，就得到了新的Mojette投影Mojp，q(2) = (RFT0 (3)_M ojp；q(l) ? 1)/0.9927〇
[0129] 重复这一步骤，直到将所有的Mojette投影求出。
[0130] 此处的转换矩阵\2如下所示：
[0132] 各Mojette投影的求解过程如下：
[0134] 至此，已将投影角度0 =tar^q/p下的Radon投影转化为投影矢量（p,q)下的Mojette投影。
【主权项】
1.基于固定分辨率条件下的离散Radon和Mojette投影转换方法，以计算机层析成像 算法为基础，通过构造合理的成像系统来逼近Radon和Mo jette投影场景，并在详细分析了 Mojette投影和Radon投影的关系之后，给出了在固定分辨率前提下各投影角度下的Radon 投影转化为对应的投影矢量下的Mojette投影的充要条件和具体的转化算法，其特征是， 步骤如下： Sl设置投影参数，即探测器分辨率和重建图像分辨率；设探测器分辨率为BXB， 重建图像分辨率为MXN，其中，重建图像分辨率和投影矢量（p，q)的选择需要满足 max {(N-I) I Pi I + (M-I) I qi I +1}〈B ; S2选择投影矢量（p，q)，计算出相应的投影角度Θ itar^q/p，并在选定的投影角度 下采集Radon投影； S3在投影矢量（p, q)下建立相应的0阶B样条插值投影系统System_rad0，需要明确 两点，其一为穿过每个像素点的投影射线的投影位置，其二为穿过每个像素点的投影射线 的线段权值； S3. 1遍历每个像素点，明确每个像素点由哪些投影射线穿过，并计算出这些投影射线 打在探测器上的有效范围； 设重建像素的实际物理尺寸为ObjPixel，探测器像元的实际物理尺寸为DetCXDSize， 当DetCXDSizKObjPixel时，说明在一个矩形区域形状的离散像素点范围内，有不止一条 射线穿过该像素点；遍历每一个像素 （i, j)，其中有l〈i〈M, l〈j〈N，并计算出穿过该像素点 的射线打在探测器上的绝对物理偏移值，此处用Xr_up和Xr_down来标记这个有效范围的 最小值和最大值，并计算出与Xr_up相对应的探测器像元标号最小值bin_up和与Xr_down 相对应的探测器像元标号最大值bin_down，投影射线打在探测器上的有效范围的计算公式 如式（12)所示：其中，（(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)所示： xr= [ ((i-1/2)-M/2) · cos ( Θ ) + ((j-l/2)-N/2) ?sin(0)] · ObjPixel (13) bin = (xr+D/2 · ObjPixel)/DetCCDSize+1 其中，为待重建图像的对角线像素点个数； S3. 2其次，需要求出每一条打在探测器中心处的投影射线穿过各像素的线段长度；设 在所有投影角度下穿过像素中心点的射线的线段权值为1，而对穿过该像素其它位置的的 射线，则需要根据这些射线与中心射线之间的垂直距离I ΔΧΓ|，计算出其线段权重值αρ(?， 其具体的计算公式如式（14)所示：其中，ObjPixel · (|cos Θ |-|sin Θ |)/2为一个临界距离值，当没有穿过中心点的射线 与穿过中心点的射线之间的垂直距离I ΔΧι·|小于等于这一个临界值时，穿过像素点的线 段长度是最长的，将其记为1;与此类似，ObjPixel · (Icos Θ | + | sin Θ |)/2也是一个临界距 离值，它标记着当没有穿过中心点的射线与穿过中心点的射线之间的垂直距离I ΔΧΓ|大 于等于这一个临界值时，穿过像素点的线段长度是该投影角度下最短，已不在该像素点的 范围之内，将其记为〇; 并将计算出的权重值存入到系统矩阵System_rad 0中相应的位置bin处，即，其中N代表图像的列数；System_rad θ 中每一行代表一条具体的投影射线穿过图像中所有点的线段权值矢量，每一列代表所有投 影射线穿过图像中一点的线段权值矢量； S4根据已求得的系统矩阵Systenurad0，来进行投影转换； 此处需要注意，由于一个像素里通常都有超过1条射线穿过，具体的射线条数取决于 P = ObjPixel/DetCCDSize，即像素物理大小和探测器实际大小的比值，故而对于Mojette 采样来说，Radon采样就存在大量冗余，例如，仅穿过一个像素的Radon采样射线就不止一 条，而此处只需取其中一条进行计算，即求出了对应的Mojette投影； 因为系统方程Systenurad0的一行代表该投影角度下的一条投影射线，选择射线其实 就是在选择系统方程中的行矢量，选择这样的射线需要满足以下两个条件：1.该行矢量非 零；2.该行矢量中存在未求出的像素值，即该行矢量包含已求出的像素之外的新像素，从 而避免无用的重复计算； S4. 1导入由S3中生成的系统矩阵System_rad0，找出系统方程中第一个非零行，且该 行矢量中只有一个非零分量，即意味着该射线只穿过了一个像素，投射在探测器上，将探测 器上的投影值视为该像素灰度值和射线穿过像素的线段权值的乘积，则Moj^(I)利用下式 计算： Mojp；q(l) = RFT0 (l)/ap；q(l, 1) 其中，RFT0(I)表示Moj1^(I)的对应位置的Radon投影值，a M(l，l)表示投影值为 RFT0 (1)的射线穿过的对应像素的线段权重值； S4. 2寻找用来求出下一个Mojette投影的Radon投影射线；在S4. 1已求出的Mojette 投影Mojp，q(l)的基础上，若想求出下一个Mojette投影Mojp， q (2)，需要找到形如RFT0 (2)= Mojp，q(2) · · Ctp q(2, 1)的 Radon 投影，即，由且仅由已求得的 Mojette 投影Mojp，q(1)和待求Mojette投影Mo jp，q(2)构成的Radon投影，若所有的Radon投影都不 满足这一条件，则该待求的Mojette投影Mo jp，q(2)无法被求出； S4. 2. 1首先，找到待求Mojette投影Mojp，q(2)中穿过的一个像素点的索引坐标（i, j)， 并计算出其在列矢量化后的向量X中的列标I = (i-1) · N+j ; S4. 2. 2找到系统矩阵System_rad0的第1列，遍历System_rad 0 (:，I)，找出其中非 零值对应的行数范围，这相当于遍历所有的投影射线，并找出穿过（i，j)这一点，权重系 数值不为零的投影射线，设这些投影矢量在投影矩阵中的行标最小值为bin_up，最大值为 bin_down ;在这一范围内，遍历每一条投影射线在投影矩阵System_rad0中对应的横行， 即在 bin_up < bin < bin_down 范围内遍历 System_rad0 (bin,:)，若系统矩阵行 System_ rad0 (bin,:)中的非零系数处对应的像素点集，由且仅由Mojm(I)和Mojp，q(2)中穿过的像 素点构成，即1^1'0(2)=厘〇几(1(2).€%， (1(2,2)+厘〇夂(1(1).€%，(1(2，1)，则该投影射线即为 候选的更新射线； S4. 2. 3在满足这一要求的射线当中，选择α ^(2, 2)值最大的Radon投影射线作为下 一个用来更新Mojette投影的射线； S4. 3利用求得的Mojette投影来更新Radon投影，即在该Radon投影中减去已求得的 Mojette投影，再除以线段权值，就得到了新的Mojette投影： Mojp；q(2) = (RFT0 (2)-Mojp；q(l) * a p；q(2, I))/a p；q(2, 2)； 重复这一步骤，直到将所有的Mojette投影求出； 即以式（15)为标准，寻找用来更新Mojette投影的Radon投影：其中，j基于已经求出的Mojette投影，从式（15)中推导出各Mojette投影的转换公式，即如 式（16)所示：至此，已将投影角度Θ = tar^q/p下的Radon投影转化为投影矢量（p, q)下的Mojette 投影。
【专利摘要】基于固定分辨率条件下的离散Radon和Mojette投影转换方法，属于图像处理领域的计算机成像技术领域。其特征是将实际成像条件下的稀疏角度离散Radon投影转化为在离散域精确重建的Mojette投影。以计算机层析成像算法为基础，通过构造合理的成像系统来逼近Radon和Mojette投影场景，并在详细分析了Mojette投影和Radon投影的关系之后，给出了在固定分辨率前提下各投影角度下的Radon投影转化为对应的投影矢量下的Mojette投影的具体算法。本发明的效果和益处是在不同的稀疏投影角度下，不需要改变Radon投影图像的分辨率大小，克服了Mojette投影域分辨率随着投影矢量变化和Radon投影分辨率不变之间的转换障碍，搭建起这两种投影之间的转换桥梁。
【IPC分类】G06T3/20
【公开号】CN104899827
【申请号】CN201510274374
【发明人】孙怡, 李梦婕, 蒋敏
【申请人】大连理工大学
【公开日】2015年9月9日
【申请日】2015年5月26日
最新回复(0)