CT断层成像系列07——等距扇束前向投影和FBP重建(附Matlab代码)

时间地点_内蒙古呼和浩特20260215

等距扇束前向投影数据仿真

下图图1为等距扇束投影系统的几何结构说明图:

​ \(D_1D_2\)为探测器阵所在位置. \(S_0B\)为某一射线, 与探测器阵相较于点\(B\), 探测器阵中点为\(Q\). 为简化推导过程中的数学表达式，可设想将\(D_1D_2\)平移到正好穿过坐标原点的位置\(D_1'D_2'\)。它是\(D_1D_2\)的镜像，二者关系可由扇形的几何尺寸确定，故射线的相对位置也可由\(OA\)定出，线段\(OA\)的长度是探测器\(D_1D_2\)上的距离。因此，射线的投影函数可记为\(p^f_{\beta}(s)=p^f(s,B)\).

​ 射线\(S_0B\)可看作上章所述的平行射束中的一条射线，由\((t,\theta)\)唯一确定。\(t,\theta\)与\(\beta,s\)之间的关系为:

\[\theta=\beta+\gamma=\beta+\arctan\frac{s}{D}\\ t=s\cos\gamma=\frac{sD}{\sqrt{D^2+s^2}} \]

​ 与等角扇束一样, 将上面两个变换关系代入平行束前向投影S-L头公式就行了(博客06的\(p(t,\theta)\)这少写个\(\rho\)). 平行束S-L头前向投影公式为:

\[p(t,\theta)=\rho\frac{2AB\sqrt{r^2-d^2}}{r^2}\\ r^2=A^2\cos^2(\theta-\alpha)+B^2\sin^2(\theta-\alpha)\\ d=t-x_0\cos\theta-y_0\sin\theta \]

​ 代入后得到等距扇束的S-L头的前向投影公式:

\[p^f(s,B)=\rho\frac{2AB\sqrt{r^2-d^2}}{r^2}\\ r^2=A^2\cos^2(\beta+\arctan\frac{s}{D}-\alpha)+B^2\sin^2(\beta+\arctan\frac{s}{D} -\alpha)\\ d=\frac{sD}{\sqrt{D^2+s^2}}-x_0\cos(\beta+\arctan\frac{s}{D})-y_0\sin(\beta+\arctan\frac{s}{D}) \]

上代码:

%======主程序======%
clc
clear
close all
N=256;
N_d=256;
beta=0:1:359;
SOD=250;
t_max=sqrt(2)*0.5*N;
delta_dd=SOD*t_max/sqrt(SOD^2-t_max^2)/(N/2);
dd=delta_dd*(-N_d/2+0.5:N_d/2-0.5);
I=phantom(N);
P= medfuncFanBeamDistanceForwardProjection(N,beta,SOD,N_d,dd);
figure(1);
imshow(I,[0 1]);
figure(2);
imagesc(P),colormap(gray),colorbar;%======子程序1======%
%等距扇束前向投影
function P=medfuncFanBeamDistanceForwardProjection(N,beta,SOD,N_d,dd);
shep =[1 .69 .92 0 0 0-.8 .6624 .8740 0 -.0184 0-.2 .1100 .3100 .22 0 -18-.2 .1600 .4100 -.22 0 18.1 .2100 .2500 0 .35 0.1 .0460 .0460 0 .1 0.1 .0460 .0460 0 -.1 0.1 .0460 .0230 -.08 -.605 0.1 .0230 .0230 0 -.606 0.1 .0230 .0460 .06 -.605 0];
rho = shep(:,1).';
ae=0.5*N*shep(:,2).';
be = 0.5*N*shep(:,3).';
xe = 0.5*N*shep(:,4).';
ye = 0.5*N*shep(:,5).';
alpha = shep(:,6).';
alpha = alpha*pi/180;
beta=beta*pi/180;
beta_num=length(beta);
P=zeros(N_d,beta_num);
for k1=1:beta_numtheta=beta(k1)+atan(dd/SOD);P_beta=zeros(1,N_d);for k2=1:length(xe)rsq=(ae(k2)*cos(theta-alpha(k2))).^2+(be(k2)*sin(theta-alpha(k2))).^2;dsq=(SOD*dd./sqrt(SOD^2+dd.^2)-xe(k2)*cos(theta)-ye(k2)*sin(theta)).^2;temp=rsq-dsq;ind=temp>0;P_beta(ind)=P_beta(ind)+rho(k2)*(2*ae(k2)*be(k2)*sqrt(temp(ind)))./rsq(ind);endP(:,k1)=P_beta.';
end
end

运行结果如下:

等距扇束滤波反投影算法公式推导

先推公式. 思路与系列07的等角扇束精简版一致, 推导平行束重建公式与扇束重建公式的变换关系, 平行束FBP重建公式为:

\[f(x,y)=\int^{\pi}_{0}p(t,\theta)*h(t)|_{t=x\cos\theta+y\sin\theta}d\theta \]

还是平行束FBP重建公式, 但是重建图像的像素点\((x，y)\)换成极坐标\((r,\phi)\)表示, 且旋转角度从180变为360(因为扇束对称角度拍的不一样了, 180°不够):

\[f(r,\phi)=\frac{1}{2}\int^{2\pi}_{0}\int^{\infty}_{-\infty}p(t'',\theta)h(t-t'')|_{t=r\cos(\phi-\theta)}dt''d\theta \]

先上图03, 等距扇束相关变量的几何关系:

和系列06一样, 依旧从\(h(t-t'')\)的滤波窗函数为突破口找到补偿修正系数, 因此要算\(t\)和\(t''\), \(t''\)为:

\[t''=s''\cos\gamma=\frac{s''D}{\sqrt{D^2+s''^2}} \]

t要按右下角积分选择来(图中\(t\)写成\(t'\)了大概),所以有:

\[t=r\cos(\theta-\phi)=r\cos(\phi-\theta)=r\cos(\phi-\beta-\gamma'')\\=r\cos(\phi-\beta)\cos\gamma''-r\sin(\beta-\phi)\sin\gamma'' \]

\(r\sin(\beta-\phi)\)和\(r\cos(\beta-\phi)\)都能直接展开成\(x,y\)的形式, 但这样就不好化简了, 化简的思路是先用\(r\cos(\beta-\phi)\)和相似三角形的关系一通化简到只剩\(r\sin(\beta-\phi)\)再展成\(x,y\)的形式.

因为\(r\cos(\beta-\phi)=\overline{EP}\), 相似三角形有:

\[\frac{s_1}{\overline{EP}}=\frac{D}{D+r\sin(\beta-\phi)} \]

所以

\[r\cos(\beta-\phi)=s_1\frac{D+r\sin(\beta-\phi)}{D} \]

由几何图可知:

\[\cos\gamma''=\frac{D}{\sqrt{D^2+s''^2}}\\ \sin\gamma''=\frac{s''} {\sqrt{D^2+s''^2}} \]

将上面三个代入\(t\)的表达式可得:

\[t=s_1\frac{D+r\sin(\beta-\phi)}{D}\frac{D}{\sqrt{D^2+s''^2}}-r\sin(\beta-\phi)\frac{s''} {\sqrt{D^2+s''^2}} \]

再把\(t''\)带进去算\(t-t''\):

\[t-t''=s_1\frac{D+r\sin(\beta-\phi)}{\sqrt{D^2+s''^2}}-r\sin(\beta-\phi)\frac{s''} {\sqrt{D^2+s''^2}}-\frac{s''D}{\sqrt{D^2+s''^2}}\\ =(s-s'')\frac{UD}{\sqrt{D^2+s''^2}}\\ U=\frac{D+r\sin(\beta-\phi)}{D} \]

\(U\)里外多放个D可能是最后化简后更好看.

然后就可以求解\(h(t-t'')\)了, 来得到滤波函数的补偿参数:

\[h(t-t'')=h[(s_1-s'')\frac{UD}{\sqrt{D^2+s''^2}}]\\ =\int^{\infty}_{-\infty}|\omega|e^{i2\pi\omega(s_1-s'')\frac{UD}{\sqrt{D^2+s''^2}}}d\omega \]

引入变量

\[\omega'=\omega\frac{UD}{\sqrt{D^2+s''^2}} \]

后, \(h(t-t'')\)就可以写成:

\[h(t-t'')=h[(s_1-s'')\frac{UD}{\sqrt{D^2+s''^2}}]\\ =\int^{\infty}_{-\infty}\frac{|\omega'|}{\frac{UD}{\sqrt{D^2+s''^2}}}e^{i2\pi\omega'(s_1-s'')}\frac{d\omega'}{\frac{UD}{\sqrt{D^2+s''^2}}}\\ =\frac{D^2+s''^2}{U^2D^2}h(s_1-s'') \]

卷积核的补偿修正参数拿到了.

再看一眼我们的平行束重建公式, 我们要把它改造成等角扇束的摸样, \(h(t-t'')\)搞定已经完成了大半:

\[f(r,\phi)=\frac{1}{2}\int^{2\pi}_{0}\int^{\infty}_{-\infty}p(t'',\theta)h(t-t'')|_{t=r\cos(\phi-\theta)}dt''d\theta \]

然后我们再把公式里的\(dt''d\theta\)换成\(ds''d\beta\)就能找齐三个补正参数了, 使用雅可比行列式易得:

\[|detJ|=\begin{vmatrix} \frac{\partial t''}{\partial s''} \frac{\partial t''}{\partial \beta}\\ \frac{\partial \theta}{\partial s''} \frac{\partial \theta'}{\partial \beta} \end{vmatrix} =\frac{D^3}{\sqrt{(D^2+s''^2)^3}} \]

然而事实上一点也不易得, 这雅可比行列式给我解的汗出来了, 全靠豆包豆指导解围, 不感兴趣可以跳过这段雅可比行列式的推导. 先直接列出四个小项的结果:

\[\frac{\partial t''}{\partial s''} = \frac{D \cdot \sqrt{D^2+s''^2} - s'' D \cdot \frac{s''}{\sqrt{D^2+s''^2}}}{D^2+s''^2}\\ = \frac{D(D^2+s''^2) - D s''^2}{(D^2+s''^2)^{3/2}} = \frac{D^3}{(D^2+s''^2)^{3/2}}\\ \frac{\partial t''}{\partial \beta} = 0\\ \frac{\partial \theta}{\partial s''} = \frac{1}{1+\left(\frac{s''}{D}\right)^2} \cdot \frac{1}{D} = \frac{D}{D^2+s''^2}\\ \frac{\partial \theta}{\partial \beta} = 1 \]

计算雅可比行列式

\[|detJ| = \begin{vmatrix} \frac{D^3}{(D^2+s''^2)^{3/2}} & 0 \\ \frac{D}{D^2+s''^2} & 1 \end{vmatrix} = \frac{D^3}{(D^2+s''^2)^{3/2}} \cdot 1 - 0 \cdot \frac{D}{D^2+s''^2} = \frac{D^3}{(D^2+s''^2)^{3/2}} \]

最终答案：\(\boldsymbol{|detJ| = \frac{D^3}{(D^2+s''^2)^{3/2}}}\)

最复杂的是第一个小项$\frac{\partial t''}{\partial s''} $得再单独拿出来说, 用的商的求导法则. 商 = 分子 ÷ 分母 \(y=\frac{u}{v}\) 它的导数公式是：

\[y' = \frac{u'v - uv'}{v^2} \]

默认是对\(s''\)求导, 就不写出来了, 设：

\[u = D s'',\quad v = \sqrt{D^2 + s''^2} = (D^2 + s''^2)^{\frac12} \]

求导：

\[u' = D\\ v' = \frac12 (D^2 + s''^2)^{-\frac12} \cdot 2s'' = \frac{s''}{\sqrt{D^2 + s''^2}} \]

代入商法则：

\[ \frac{\partial t''}{\partial s''} = \frac{u'v - uv'}{v^2} = \frac{D\cdot\sqrt{D^2+s''^2} - D s''\cdot \frac{s''}{\sqrt{D^2+s''^2}}}{D^2+s''^2} \]

再化简, 分子通分：

\[D\cdot\sqrt{D^2+s''^2} - \frac{D s''^2}{\sqrt{D^2+s''^2}} = \frac{D(D^2+s''^2) - D s''^2}{\sqrt{D^2+s''^2}} = \frac{D^3}{\sqrt{D^2+s''^2}} \]

所以:

\[\frac{\partial t''}{\partial s''} = \frac{\;\displaystyle\frac{D^3}{\sqrt{D^2+s''^2}}\;}{D^2+s''^2} = \frac{D^3}{(D^2+s''^2)^{3/2}} \]

雅可比推导结束, 我们得到了投影补正系数.

投影值在代码里就是迭代的变量变了, 但在推公式里要特意说一嘴不然不能直接从\(p(t'',\theta)\)变\(p^f(s'',\beta)\):

\[p(t'',\theta)=p(\frac{s''D}{\sqrt{D^2+s''^2}},\beta+\arctan\frac{s}{D})=p^f(s'',\beta) \]

再看一次平行束重建公式:

\[f(r,\phi)=\frac{1}{2}\int^{2\pi}_{0}\int^{\infty}_{-\infty}p(t'',\theta)h(t-t'')|_{t=r\cos(\phi-\theta)}dt''d\theta \]

把我们求得的\(h(t-t'')\)和\(dt''d\theta\)代入, 得:

\[f(r,\phi)=\frac{1}{2}\int^{2\pi}_{0}\int^{\infty}_{-\infty}p^f(s'',\beta)\frac{D^2+s''^2}{U^2D^2}h(s_1-s'')|_{s=s_1}\frac{D^3}{\sqrt{(D^2+s''^2)^3}}ds''d\beta\\ =\int^{2\pi}_{0}\int^{\infty}_{-\infty}\frac{1}{U^2}[p^f(s'',\beta)\frac{D}{\sqrt{D^2+s''^2}}][\frac{1}{2}h(s_1-s'')]|_{s=s_1}ds''d\beta \]

这就是等距扇束公式完全体, \(s_1\)用公式8可得:

\[s_1=\frac{Dr\cos(\beta-\phi)}{D+r\sin(\beta-\phi)}=\frac{x\cos\beta+y\sin\beta}{U} \]

我们是把重建图像的像素点的极坐标形式下推导的, 最终公式还需换回笛卡尔坐标. 还有三大补正参数和真正的射线标号的计算我们和代码一起讲.

等距扇束滤波反投影算法代码实现

先上代码:

%======主程序======%
clc
clear
close all
N=256;%图像大小
N_d=256;%探测器个数
beta=0:359;%角度旋转
SOD=250;%焦距
t_max=sqrt(2)*0.5*N;
delta_dd=SOD*t_max/sqrt(SOD^2-t_max^2)/(N/2);
dd=delta_dd*(-N_d/2+0.5:N_d/2-0.5);
I=phantom(N);
P=medfuncFanBeamDistanceForwardProjection(N,beta,SOD,N_d,dd);
fh_RL=medfuncFanBeamRLFilter2(N_d,delta_dd);
rec_RL=medfuncFanBeamDistanceFBP(P,fh_RL,beta,SOD,N,N_d,delta_dd);
figure;
subplot(1,2,1),imshow(I,[]),title('250×250头模型原始图像');
subplot(1,2,2),imshow(rec_RL,[]),title('等间距FBP算法重建R_L函数');%======子程序1======%
%等距扇束前向投影
function P=medfuncFanBeamDistanceForwardProjection(N,beta,SOD,N_d,dd);
shep =[1 .69 .92 0 0 0-.8 .6624 .8740 0 -.0184 0-.2 .1100 .3100 .22 0 -18-.2 .1600 .4100 -.22 0 18.1 .2100 .2500 0 .35 0.1 .0460 .0460 0 .1 0.1 .0460 .0460 0 -.1 0.1 .0460 .0230 -.08 -.605 0.1 .0230 .0230 0 -.606 0.1 .0230 .0460 .06 -.605 0];
rho = shep(:,1).';
ae=0.5*N*shep(:,2).';
be = 0.5*N*shep(:,3).';
xe = 0.5*N*shep(:,4).';
ye = 0.5*N*shep(:,5).';
alpha = shep(:,6).';
alpha = alpha*pi/180;
beta=beta*pi/180;
beta_num=length(beta);
P=zeros(N_d,beta_num);
for k1=1:beta_numtheta=beta(k1)+atan(dd/SOD);P_beta=zeros(1,N_d);for k2=1:length(xe)rsq=(ae(k2)*cos(theta-alpha(k2))).^2+(be(k2)*sin(theta-alpha(k2))).^2;dsq=(SOD*dd./sqrt(SOD^2+dd.^2)-xe(k2)*cos(theta)-ye(k2)*sin(theta)).^2;temp=rsq-dsq;ind=temp>0;P_beta(ind)=P_beta(ind)+rho(k2)*(2*ae(k2)*be(k2)*sqrt(temp(ind)))./rsq(ind);endP(:,k1)=P_beta.';
end
end%======子程序2======%
%等距扇束RL滤波核
function fh_RL=medfuncFanBeamRLFilter2(N_d,delta_dd)
fh_RL=zeros(N_d,1);
for k1=1:N_dfh_RL(k1)=-1/(2*pi*pi*((k1-N_d/2-1)*delta_dd)^2);if mod(k1-N_d/2-1,2)==0fh_RL(k1)=0;end
end
fh_RL(N_d/2+1)=1/(8*delta_dd^2);
end%======子程序3======%
%等距扇束滤波反投影重建函数
function rec_RL=medfuncFanBeamDistanceFBP(P,fh_RL,beta,SOD,N,N_d,delta_dd)
dd=delta_dd*(-N/2+0.5:N/2-0.5);
beta=beta*pi/180;
beta_num=length(beta);
MX=N;MY=N;
roi=N*[-0.5 0.5 -0.5 0.5];
hx=(roi(2)-roi(1))/(MX-1);
xrange=roi(1)+hx*[0:MX-1];
hy=(roi(4)-roi(3))/(MY-1);
yrange=flipud((roi(3)+hy*[0:MY-1])');
x1=ones(MY,1)*xrange;
x2=yrange*ones(1,MX);
rec_RL=zeros(MY,MX);
for m=1:beta_numalphaj=beta(m);RF1=P(:,m).*(SOD./sqrt(SOD^2+dd.^2))';C_RL=conv(RF1,fh_RL,'same');aj=[cos(alphaj);sin(alphaj)];U=(SOD+x1.*aj(2)-x2.*aj(1))/SOD;t=real((x1.*aj(1)+x2.*aj(2))./U)/delta_dd;k=floor(t);u=t-k;k=max(1,k+N_d/2+1);k=min(k,N_d-1);P_RL=((1-u).*C_RL(k)+u.*C_RL(k+1));rec_RL=rec_RL+P_RL./U^2*2*pi/beta_num;
end
end

运行结果如下:

等距扇束与等角扇束一样, 除了三处修正和平行束FBP几乎一样.

第一步, 投影函数补正, 探测器等距间隔\(\Delta d\), 则有:

\[P_{\beta}(s)=P_{\beta}(n\Delta d)=p_{\beta}(n\Delta d)\frac{D}{\sqrt{D^2+(n\Delta d)^2}} \]

此处位于代码第87行:

RF1=P(:,m).*(SOD./sqrt(SOD^2+dd.^2))';

第二步, 滤波投影修正, 公式28里看到就乘个\(\frac{1}{2}\)就行, 平行束Ram-Lak滤波核为:

\[h(n\Delta\alpha)=\begin{cases} \frac{1}{4(\Delta d)^2},\qquad n=0\\ 0,\qquad\qquad n为偶数\\ -\frac{1}{(\pi n\Delta d)^2},\quad n=奇数\\\end{cases} \]

等距扇束滤波核:

\[h'(n\Delta\alpha)=\begin{cases} \frac{1}{8(\Delta d)^2},\qquad n=0\\ 0,\qquad\qquad n为偶数\\ -\frac{1}{2(\pi n\Delta d)^2},\quad n=奇数\\\end{cases} \]

对应代码61-67行:

for k1=1:N_dfh_RL(k1)=-1/(2*pi*pi*((k1-N_d/2-1)*delta_dd)^2);if mod(k1-N_d/2-1,2)==0fh_RL(k1)=0;end
end
fh_RL(N_d/2+1)=1/(8*delta_dd^2);

第三步, 加权反投影补正. 最后累加全角度的时候乘\(\frac{1}{U^2}\)就行, 真正运算时从极坐标变回笛卡尔, 此时的\(U\)为:

\[U=\frac{D+x\sin\beta_i-y\cos\beta_i}{D} \]

对应代码在行90和96:

U=(SOD+x1.*aj(2)-x2.*aj(1))/SOD;
rec_RL=rec_RL+P_RL./U^2*2*pi/beta_num;

还有最为核心的一步没提, 我们要计算每个像素点的射线编号, 然后才能进行线性插值:

\[n_0=\frac{s_1}{\Delta d}=\frac{\frac{x\cos\beta_i+y\sin\beta_i}{U}}{\Delta d} \]

对应代码在行91:

t=real((x1.*aj(1)+x2.*aj(2))./U)/delta_dd;

结束~~~