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

摘要

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

主权项

基于固定分辨率条件下的离散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)·\cos \left(\theta \right)+((j-1)-N/2)·\sin \left(\theta \right)]·\mathrm{ObjPixel}\\ \mathrm{Xr}\_\mathrm{down}=[(\left(i\right)-M/2)·\cos \left(\theta \right)+(\left(j\right)-N/2)·\sin \left(\theta \right)]·\mathrm{ObjPixel}\end{array}\right.,\theta \le 90$$\left\{\begin{array}{c}\mathrm{Xr}\_\mathrm{up}=[(\left(i\right)-M/2)·\cos \left(\theta \right)+((j-1)-N/2)·\sin \left(\theta \right)]·\mathrm{ObjPixel}\\ \mathrm{Xr}\_\mathrm{down}=[((i-1)-M/2)·\cos \left(\theta \right)+(\left(j\right)-N/2)·\sin \left(\theta \right)]·\mathrm{ObjPixel}\end{array}\right.\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}·\left(\right|\cos \theta |-|\sin \theta \left|\right)}{2},\\ {\alpha }_{p,q}=1\\ \mathrm{if}\frac{\mathrm{ObjPixel}·\left(\right|\cos \theta |-|\sin \theta \left|\right)}{2}<\left|\mathrm{\Delta Xr}\right|\le \frac{\mathrm{ObjPixel}·\left(\right|\cos \theta |+|\sin \theta \left|\right)}{2}\\ {\alpha }_{p,q}=\frac{-2·\left|\mathrm{\Delta Xr}\right|+\left(\right|\cos \theta |+|\sin \theta \left|\right)·\mathrm{ObjPixel}}{2·\mathrm{ObjPixel}·\min \left(\right|\cos \theta ,\sin \theta \left|\right)}\\ \mathrm{if}\left|\mathrm{\Delta Xr}\right|>\frac{\mathrm{ObjPixel}·\left(\right|\cos \theta |+|\sin \theta \left|\right)}{2}\\ {\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}}_{\theta }(\mathrm{radbin},(i-1)·N+j)={\alpha }_{i,j,\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}}_{\theta }\left(1\right)={\mathrm{Moj}}_{p,q}\left(1\right)·{\alpha }_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{RFT}}_{\theta }\left(2\right)={\mathrm{Moj}}_{p,q}\left(2\right)·{\alpha }_{p,q}\left(\mathrm{2,2}\right)+{\mathrm{Moj}}_{p,q}\left(1\right)·{\alpha }_{p,q}\left(\mathrm{2,1}\right)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\theta }\left(k\right)={\mathrm{Moj}}_{p,q}\left(k\right)·{\alpha }_{p,q}(k,k)+{\Sigma }_{j=1}^{j\le k-1}{\mathrm{Moj}}_{p,q}(k-j)·{\alpha }_{p,q}(k-j,k-j)\\ .\\ .\\ .\\ {\mathrm{RFT}}_{\theta }\left(\mathrm{DetRN}\right)={\Sigma }_{k=1}^{k\le \mathrm{DetRN}}{\mathrm{Moj}}_{p,q}\left(k\right)·{\alpha }_{p,q}(k,k)\end{array}\right.---\left(15\right)$其中，${\mathrm{Moj}}_{p,q}\left(k\right)=\underset{k=0}{\overset{P-1}{\Sigma }}\underset{l=0}{\overset{Q-1}{\Sigma }}\mathrm{pixel}(i,j)\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}}_{\theta }\left(1\right)/{\alpha }_{p,q}\left(\mathrm{1,1}\right)\\ {\mathrm{Moj}}_{p,q}\left(2\right)=({\mathrm{RFT}}_{\theta }\left(2\right)-{\mathrm{Moj}}_{p,q}\left(1\right)·{\alpha }_{p,q}\left(\mathrm{2,1}\right))/{\alpha }_{p,q}\left(\mathrm{2,2}\right)\\ .\\ .\\ .\\ {\mathrm{Moj}}_{p,q}\left(k\right)=({\mathrm{RFT}}_{\theta }\left(k\right)-{\Sigma }_{j=1}^{j\le k-1}{\mathrm{Moj}}_{p,q}(k-j)·{\alpha }_{p,q}(k-j,k-j))/{\alpha }_{p,q}(k,k)\\ .\\ .\\ .\end{array}\right.---\left(16\right)$至此，已将投影角度θ＝tan^‑1q/p下的Radon投影转化为投影矢量(p,q)下的Mojette投影。