CT图像重建

（ 2009—2010学年 第 一 学期 ）

课程名称：医学成像系统与放射治疗装置 开课实验室： 3208 2008 年 12 月24 日 年级、专业、班 实验项目名称 教师评 学号 CT图像重建 姓名 指导教师 成绩 刘利军 语 教师签名： 年 月 日 一、实验目的与意义

医学成像技术是生物医学工程专业的一门重要的专业课程,课程主要涉及X光仪器，CT仪器，MRI仪器和核医学仪器的工作原理及成像方法。其中CT算法的出现又为后来数字化医学成像技术的发展提供了基础。该门课程为生物医学工程专业的专业基础课。

CT技术是医学成像系统中的一种重要手段。它通过特定的算法，利用计算机的高速运算功能，可以在短时间内快速呈现人体断层图像。让学生练习CT图像的重建有助于学生理解CT算法的内容，熟悉数字图像重建的过程。同时也能培养学生的团队精神和解决实际问题的能力。

二、实验算法原理

1、MATLAB处理数字图像的基本函数； 2、X-CT三维图像重建的基本算法。

CT图象重建有四种基本的算法:矩阵法,迭代法,傅立叶算法,反投影算法.我们采用的方法为卷积反投影. 卷积反投影有:平行光束投影的卷积反投影算法, 等角扇形光来投影的重建算法. 1）.平行光束投影的卷积反投影算法 从投影重建三维物体的图像，就是重建一个个横断面。这样三堆图像的重建就归结为二维图象的重建。二维图像的重建问题可以从数学上描述如下。

假定g(x,y)表示一个二维的未知函数，通过g(x,y) 的直线称为光钱(见图2．1)。沿光线g(x,y)的积分称作光线积分。沿相同方向的一组光线积分，就构成一个投影。图2．1中垂直于直线CC (与X轴夹角为)的光线所形成。

精选

'

图2.1 g(x,y)在方向的投影P(t)

的投影P(t)，称之为g(x,y)在方向的投影。光线积分和投影在数学上可以定义如下：

在图2.1中直线AB的方程为：

XcosYsint1 (2.1) 其中t1是AB到原点的距离，g(x,y)沿AB的积分为：

P(t1)g(x,y)dsg(x,y)(xcosysint1)dxdy (2.2)

AB对于给定的，g(x,y)在方向的投影P如果g(x,y)在各个方向的投影已知，g(x,y)(t)是t的函数。就可以唯一确定。下面就讨论卷积反投影重建算法。

假定投影方向，如图2.2，将坐标(x,y)旋转角（逆时针方向）形成坐标系(t,s)。g(x,y)在(t,s)坐标系中为g(t,s)。

精选

图2.2 傅立叶切片定理示意图

坐标系(t,s)与(x,y)之间的关系为： tcossinssincosxy 显然

Ptg(t,s)ds 令S(w)为P(t)的傅立叶变换则 S(w)P(t)exp(j2wt)dt g(t,s)exp(j2wt)dsdt 将上式变换到(x,y)坐标系中，注意到变换的可比行列式

tt tysscossinsincos1 xy从而得到： S(w)g(x,y)exp[j2(xcosysin)]dxdy g(x,y)exp[j2(uxvy)]dxdy 其中

ucosvsin 精选

2.3）2.4）2.5）2.6）2.7）2.8）（ （

（

（

（

（

若令g(x,y)的傅立叶变换为G(u,v)，由（2.8）可知

S()G(u,v)G(,) （2.9） 若g(x,y)的傅立叶变换为G(u,v)的极坐标表示。这说明g(x,y)在方向的投影P(t)

傅立叶变换S(w)等于G(u,v)在与u轴成角的直线上的值。这就是著名的傅立叶投影切片定理。可见在整个(u,v)平面G(u,v)可以利用各个方向的投影来得到，从而g(x,y)也可以通过求G(u,v)的傅立叶反交换的办法求得： g(x,y)变换到极坐标中 得到

g(x,y)经推导得 g(x,y)其中

txcosysin 若令

G(u,v)exp[j2(uxvy)]dudv （2.10）

ucos， 

vsin020G(,)exp[j2(xcosysin)]dd （2.11）

0S()exp(j2t)dd （2.12） Q(t)S()exp(j2t)d （2.13）

则

g(x,y)0Q(xcosysin)d (2.14)

(2.13)式右端是两频谱函数S(w)和H()的乘积的傅立叶反变换。S(w)是投影P(t) 傅立叶变换。若H()的傅立叶反变换为h(t)，则根据卷积定理有： Q(t)或

Q(t)P(t)h(t) 其中

h(t)P()h(t)d （2.15）

exp(j2t)d （2.16）

精选

当图像的频谱是有限带宽时，则上式变为

h(t)00exp(j2t)d （2.17）

由于图象及其频谱都是离散采样的， 假定图象采样间隔为， 则根据采样定理01/2。为了进行数学处理，只需知道h (t)在有限带宽上的离散采样点的值．这样我们有

1/(42)  h(n)0 （2.18）

1/n222其中n为正负整数。 （2.18）的离散形式为 Q(n)mP(m)h(nm) （2.19）

假定P(m)在m0,1,N1之外的值为0，则上式变为 Q(n)或

Q(n)m0P(m)h(nm) （2.20） P[nm]h(m) （2.21）

N1N1m(N1)其中n0,1,2,N1 从而可见为确定P(t)的N个采样点上的Q(n)的值，需要使用h(n)的2N— 1个点上的值，从n=一(N— 1)到(N — 1)。

为求得Q(n)，利用傅立叶变换计算卷积是比较快的方法，为清除循环卷积的周期交叠效应，实际上h(n)取2N个点，P(m)补0，使之有(2N—1)个元素，则P(t)在N个采样点上就避免了交叠，如果使用以2为基的FFT(快速傅立叶变换)算法， P(m)和h(n)都必需朴0至(2N一1)个元素，(2N一1)为大于等于2N—l的最小的2的整数幂。计算Q(n)的过程可以写为

Q(n)[FFT(FFT[P(n)0]FFT(h(n)0] （2.22）

其中FFT和IFFT分别表示快速傅立叶变换和反变换， 光滑窗是在滤波过程中加入的光滑因子，例如引用汉明窗 ，有时可以改进重建效果。对于各个方向的投影， 得到Q(n)之后就可以由(2.22)来计算

g(x,y)。重建步骤可以归纳为：

第一步：卷积，也称滤波，由(2．22)对每个方向计算Q(n)。 第二步。反投影，由(2．14)的近似形式

精选

g(x,y)MQ(xcosii1Miysini) （2.23）

ˆ(x,y)。M为投影个数i为投影方向角，他们均匀的分布在0~的范围内。 来计算g(x,y)的近似值g 当计算Qi(xcosiysini)时，txcosiysini，不一定在Q(n)的整离散点上，这就需要插值求得，预先将Q(n)插值加密，即最靠近的点，可以提高计算速度。

2).等角扇形光来投影的重建算法

几乎所有的快遗CT设备都是用的扇形光束。这里叙述的是等角度光束投影，如图2.3，测量投影数据的探测器等间距地分布在D1D2弧上，弧的半径为2D， D为光源到图像中心的距离。在下文中，f(r,)图象在极坐标中的表示。R()表示在方向角为的投影中位量角为的光线产生的投影数据。通过中心的光线其为0。L表示从光源到像素(r,)的距离。

图2.3 等扇形束投影重建算法中的变量

L(,,)D2r22Drsin() （2.24）

表示在方向角为的投影中通过像素(r,)的光线的位置角

(,,)tan'1PEcos()tan1 （2.25） SEDsin()图像f(r,)和扇形投影R()有下述关系

f(r,)2011''DcosrR()()h()d （2.26） 2'L2sin()精选

其中r和是投影中光线的最大位置角，从而可以得到这种重建算法的执行步骤： 第一步：投影的修改，假定投影的抽样间隔为，抽样数据R(n)通过下式进行修改， R'(n)R(n)D1cosn （2.27） 第二步：卷积（滤波）将第一步修改了的投影与响应函数g()进行卷积 g()'112()h() （2.28） 2sinQ()R(n)g(nm)d （2.29）

其离散形式为： Q(n)R(m)g(nm) （2.30）

'm1M这里也和上节一样可以加进一光滑窗，改进重建效应。 第三步：反投影，就是执行（2.30）的积分 f(r,)近似的有：

201Q()d （2.31）

L2(,,)M2 f(x,y)mL(,,)Q() （2.32）

i12'1其f(x,y)是图象f(x,y)的近似图像，xrcos，yrsin， M是投影数，假定投影均匀地分布在0~2内。为求得滤波后的投影Q，对象素(x,y)的贡献，首先由(2.24)和(2.25)计算出L和．所有

Q对(x,y)的贡献加起来再乘以2/M就得f(x,y)。

四、实验程序代码

原图象代码: p=phantom(256); imshow(p)

卷积反投影法代码：

H=[0 0 0.92 0.69 90*pi/180 1;...

0 -0.0184 0.874 0.6624 90*pi/180 -0.98;... 0.22 0 0.31 0.11 72*pi/180 -0.2;... -0.22 0 0.41 0.16 108*pi/180 -0.2;... 0 0.35 0.25 0.21 90*pi/180 0.1;... 0 0.1 0.046 0.046 0 0.2;... 0 -0.1 0.046 0.046 0 0.2;...

-0.08 -0.605 0.046 0.023 0 0.1;...

精选

0 -0.605 0.023 0.023 0 0.1;...

0.06 -0.605 0.046 0.023 90*pi/180 0.1]; angle=1; N=100; vstep=2/N; ax=N/2+1;

for k=1:(180/angle+1)

theta=(k-1)*angle*pi/180; for i=1:10

x0=H(i,1);y0=H(i,2);A=H(i,3);B=H(i,4);alpha=H(i,5);rho=H(i,6); R=0; forw=ax; back=ax;

MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(theta))^2); NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2; g(i,ax)=rho*2*A*B*MM/NN; for j=1:N/2 R=R+vstep; forw=forw+1;

MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(theta))^2); NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2; g(i,forw)=rho*2*A*B*MM/NN; R=-R;

back=back-1;

MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(theta))^2); NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2; g(i,back)=rho*2*A*B*MM/NN; R=-R; end end

radon0(k,:)=real(sum(g)); end

%%%%%%%%%%%%generate R-L function%%%%%%%%%%%%%%%% radon1=[zeros(k,N),radon0]; ax=N+1;

RL(ax)=1/(4*vstep^2); forw=ax+1; back=ax-1; for k=1:N/2 n=2*k-1;

RL(forw)=-1/(n*pi*vstep)^2; RL(back)=RL(forw); RL(forw+1)=0; RL(back-1)=0; forw=forw+2; back=back-2;

精选

end

for k=1:(180/angle+1)

radon2(k,:)=conv(radon1(k,:),RL); end

radonf=radon2(:,2*N:3*N);

%%%%%%%%%%%%generate S-L function%%%%%%%%%%%%%%%%%% radon1=[zeros(k,N),radon0]; for v=1:(2*N+1) n=v-N-1;

SL(v)=-2/(pi^2*vstep^2*(4*n^2-1)); end

for k=1:(180/angle+1)

radon2(k,:)=conv(radon1(k,:),SL); end

radonf=radon2(:,2*N:3*N); figure(1) subplot(321)

plot(1:(2*N+1),radon1(1,:)) title('投影函数（已补零）') subplot(323)

plot(1:(2*N+1),SL) title('S-L卷积函数') subplot(325)

plot(1:(4*N+1),radon2(1,:)) title('卷积结果')

Xradon1=fft(radon1(1,:)); subplot(322)

plot(1:(2*N+1),abs(Xradon1)) title('频谱') XRL=fft(SL); subplot(324)

plot(1:(2*N+1),abs(XRL)) Xradon2=fft(radon2(1,:)); subplot(326)

plot(1:(4*N+1),abs(Xradon2))

%%%%%%%%%%%%%%iradon%%%%%%%%%%%%%%%%%%%%%%%%%% for k=1:(180/angle+1)

theta=(k-1)*angle*pi/180;

C=N/2-(N-1)*(cos(theta)+sin(theta))/2; for i=1:N

R=(i-1)*cos(theta)+C; n0=floor(R);

if n0>0&&n0<(N+1) dot=R-n0;

I(i,1,k)=(1-dot)*radonf(k,n0)+dot*radonf(k,n0+1); else

精选

I(i,1,k)=nan; end

for j=2:N

R=R+sin(theta); n0=floor(R);

if n0>0&&n0<(N+1) dot=R-n0;

I(i,j,k)=(1-dot)*radonf(k,n0)+dot*radonf(k,n0+1); else

I(i,j,k)=nan; end end end end

Ifinal=sum(I,3);

Gfinal=mat2gray(Ifinal); Gfinal=imrotate(Gfinal,90); figure(2) imshow(Gfinal)

五、实验结果与分析

反投影一般步骤为：

原像 重建后图像 取投影

程序流程:

反投影重建

精选

提取图象数据 初始化 产生R-L函数 产生S-L函数 进行radon反变换 显示重建图象结果 结束

卷积反投影法结果：

精选

重建图象 原图象

六、心得体会

通过本次的实验,对CT图象重建的基本方法之一:卷积反投影,有了进一步的认识,在实验的过程中,采用的图象是经典的Shep-Logen头模型,得到的结果与原图象相比,有一定的差异,但影响不大.有待进一步的改进算法.

七、参考文献

[1] 《医用电子仪器及装备——医学成像系统及放射治疗装置类》，唐庆玉主编，清华大学电机系生物医学工程及仪器专业 [2] 《医学成像系统》，高上凯主编，清华大学出版社

[3] G．T.赫尔曼著(严洪编译)，由投影重建图象，科学出版社 [4] 董雏申、吴世法、王天童，卷积反投影法从x光图像重建三维轴对称图象，全国图基科学会议论文集，l989。

[5] G，N．Housfleld Computed Medical Imaging -Journal 0f Computed Assisted Tomography，4(5)，655—674．1980

[6] 赶荣椿等，《数字图像处理导论》，百北工业大学出版牡，西安，1995，p77—179 [7] 庄天戈，《CT原理与算法)，上海交通大学出版杜，1992，p31— 33

[8] Radon J，Uber die Bestimmumg von Funktionen much ihre integralwerte langs gewissser

Mannnmgfaltigkeiten Berichte Sachsische Akademie der Wissenschaften,Leipzig,Math．Phys．K1，1917.69：262~267

精选