CT图像重建.docx

1、CT图像重建昆明理工大学信息工程与自动化学院学生实验报告(2009 2010学年第一学期)课程名称：医学成像系统与放射治疗装置 开课实验室：3208 2008 年12月24日年级、专业、班学号姓名成绩实验项目名称CT图像重建指导教师刘利军教师 评 语教师签名：年 月 日、实验目的与意义医学成像技术是生物医学工程专业的一门重要的专业课程 ，课程主要涉及 X光仪器，CT仪器，MRI仪器和核医学仪器的工作原理及成像方法。 其中CT算法的出现又为后来数字化医学成像技术的发展提供了基础。该门课程为生物医学工程专业的专业基础课。CT技术是医学成像系统中的一种重要手段。它通过特定的算法，利用计算机的高速运算

2、功能，可以 在短时间内快速呈现人体断层图像。让学生练习 CT图像的重建有助于学生理解 CT算法的内容，熟悉数字图像重建的过程。同时也能培养学生的团队精神和解决实际问题的能力。二、实验算法原理1、 MATLA处理数字图像的基本函数；2、 X-CT三维图像重建的基本算法。CT图象重建有四种基本的算法 ：矩阵法，迭代法，傅立叶算法，反投影算法.我们采用的方法为卷积反 投影卷积反投影有：平行光束投影的卷积反投影算法 ，等角扇形光来投影的重建算法 1).平行光束投影的卷积反投影算法从投影重建三维物体的图像， 就是重建一个个横断面。这样三堆图像的重建就归结为二维图象的重建。 二维图像的重建问题可以从数学上

3、描述如下。假定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的方程为：X cos Ysin ti (2.1)其中ti是AE到原点的距离，g(x,y)沿AB的积分为：P (ti) ABg(x, y)ds g(x, y) (xcos ysin tjdxdy (22)对于给定

4、的 ，g(x, y)在 方向的投影P (t)是t的函数。如果g(x, y)在各个方向的投影已知， g(x, y)就可以唯一确定。下面就讨论卷积反投影重建算法。假定投影方向 ，如图2.2，将坐标(x,y)旋转 角(逆时针方向)形成坐标系(t,s)。 g(x,y)在(t,s)坐标系中为g(t, s)。图2.2 傅立叶切片定理示意图坐标系(t,s)与(x,y)之间的关系为：t cos sin xs sin cos y显然P t g(t,s)ds令S (w)为P (t)的傅立叶变换则(2.3)(2.4)S (w) P (t)exp( j2 wt)dtg (t, s)exp( j 2 wt)dsdt将上

5、式变换到(x, y)坐标系中，注意到变换的可比行列式(2.5)IIcos sinsin cos(2.6)从而得到：S (w)g(x, y) exo j2 (xcos ysin ) |dxdyg(x, y)exp j2 (ux vy)dxdy(2.7)其中cos sin(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轴成 角的直线上的值。这就是著名的傅立叶投影切片定理。可见在整个

6、(u,v)平面G(u,v)可以利用各个方向的投影来得到， 从而g(x, y)也可以通过求 G(u,v)的傅立叶反交换的办法求得:II变换到极坐标中cos sin得到经推导得其中若令其中当图像的频谱是有限带宽时，则上式变为(2.17 )0h(t) J |e)p(j2 t)d由于图象及其频谱都是离散采样的， 假定图象采样间隔为 ，则根据采样定理 0 1/2。为了进行数学处理，只需知道h (t)在有限带宽上的离散采样点的值这样我们有1/(4 2)h(n ) 0 (2.18)2 2 21/n2 2 2其中n为正负整数。(2.18 )的离散形式为Q (n ) P (m )h (n m) (2.19)m假

7、定P (m )在m 0,1, N 1之外的值为0,则上式变为其中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

8、 )的过程可以写为Q (n ) FFT (FFTP (n )0 FFT (h(n )0 (2.22 )其中FFT和IFFT分别表示快速傅立叶变换和反变换， 光滑窗是在滤波过程中加入的光滑因子，例如引用汉明窗，有时可以改进重建效果。对于各个 方向的投影， 得到Q (n )之后就可以由(2.22)来计算 g(x, y)。重建步骤可以归纳为: 第二步。反投影，由(2 . 14)的近似形式Mg(x, y) Qi(xcos i ysin J (2.23)M i i来计算g(x, y)的近似值0(x,y)。M为投影个数j为投影方向角，他们均匀的分布在 0的范围内。当计算Q i(xcos i ysin J时

9、，t xcos i ysin i，不一定在Q (n )的整离散点上，这就需要插值求得，预先将 Q (n )插值加密，即最靠近的点，可以提高计算速度。2).等角扇形光来投影的重建算法几乎所有的快遗CT设备都是用的扇形光束。 这里叙述的是等角度光束投影， 如图2.3，测量投影数据的探测器等间距地分布在 D1 D2弧上，弧的半径为2D, D为光源到图像中心的距离。在下文中， f(r,)图象在极坐标中的表示。 R ()表示在方向角为 的投影中位量角为 的光线产生的投影数据。通过中心的光线其 为0。L表示从光源到像素(r,)的距离。图2.3 等扇形束投影重建算法中的变量L( , , ) , D2 r2

10、2Dr sin( ) (2.24)表示在方向角为 的投影中通过像素(r,)的光线的位置角1 PE 丄 1 cos( )(,)tan tan ( 2.25)SE D sin( )图像f(r,)和扇形投影R ()有下述关系其中r和 是投影中光线的最大位置角，从而可以得到这种重建算法的执行步骤:第一步：投影的修改，假定投影的抽样间隔为 ，抽样数据R (n )通过下式进行修改，R (n ) R (n Dosn ( 2.27)第二步：卷积(滤波)将第一步修改了的投影与响应函数 g()进行卷积 1 1 2g ( ) -( ) h( ) (2.28)2 sinQ ( ) R (n )g(n m) d (2.

11、29)其离散形式为：MQ (n ) R (m )g(n m) (2.30)m 1这里也和上节一样可以加进一光滑窗，改进重建效应。第三步：反投影，就是执行(2.30 )的积分f(r,)2 1(2.31 )2 Q ( )d0 L2(,) )近似的有：f (x,y)2 M 1()(2.32 )2 Q ()m i 1 L (,)其f (x, y)是图象f(x,y)的近似图像，x rcos ， y rsin , M是投影数，假定投影均匀地分布在 02内。为求得滤波后的投影 Q，对象素(x,y)的贡献，首先由(2.24)和(2.25)计算出L和所有 Q对(x, y)的贡献加起来再乘以 2 / M就得f (

12、x, y)。四、实验程序代码原图象代码：p=pha ntom(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

13、 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(AA2*cos(theta-alpha)A2+BA2*si n(theta-alpha)A2-(

14、R-xO*cos(theta)-yO*si n(theta)2);NN=AA2*cos(theta-alpha)A2+BA2*si n(theta-alpha)A2; g(i,ax)=rho*2*A*B*MM/NN;for j=1:N/2R=R+vstep;forw=forw+1;MM=sqrt(AA2*cos(theta-alpha)A2+BA2*sin(theta-alpha)A2-(R-x0*cos(theta)-y0*sin(theta)A2);NN=AA2*cos(theta-alpha)A2+BA2*sin(theta-alpha)A2; g(i,forw)=rho*2*A*B*M

15、M/NN;R=-R;back=back-1;MM=sqrt(AA2*cos(theta-alpha)A2+BA2*sin(theta-alpha)A2-(R-x0*cos(theta)-y0*sin(theta)A2);NN=AA2*cos(theta-alpha)A2+BA2*sin(theta-alpha)A2; g(i,back)=rho*2*A*B*MM/NN;R=-R;endendradon0(k,:)=real(sum(g);end%generate R-L function%radon1=zeros(k,N),radon0;ax=N+1;RL(ax)=1/(4*vstepA2);

16、forw=ax+1;back=ax-1;for k=1:N/2n=2*k-1;RL(forw)=-1/(n*pi*vstep)A2;RL(back)=RL(forw);RL(forw+1)=0;RL(back-1)=0;forw=forw+2;back=back-2;endfor k=1:(180/angle+1)radon2(k,:)=conv(radon1(k,:),RL);endradonf=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/

17、(p2*vstepA2*(4* 门人2-1);endfor k=1:(180/angle+1)radon2(k,:)=conv(radon1(k,:),SL);endradonf=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,:);subpl

18、ot(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:NR=(i-1)*cos(theta)+C; nO=floor(R);if n00&n 0

19、0&n 0(N+1)dot=R-n0;I(i,j,k)=(1-dot)*rado nf(k,n 0)+dot*rado nf(k, n0+1); elseI(i,j,k)=nan;endendendendIfin al=sum(I,3);Gfin al=mat2gray(lfi nal);Gfin al=imrotate(Gfi nal,90);figure(2)imshow(Gfi nal)五、实验结果与分析反投影一般步骤为：原像 重建后图像取投影 反投影重建程序流程:卷积反投影法结果:重建图象 原图象六、 心得体会通过本次的实验，对CT图象重建的基本方法之一：卷积反投影，有了进一步的认识，

20、在实验的过程中， 采用的图象是经典的 Shep-Logen头模型，得到的结果与原图象相比，有一定的差异，但影响不大有待进 一步的改进算法七、 参考文献1医用电子仪器及装备一一医学成像系统及放射治疗装置类 ，唐庆玉主编，清华大学电机系生物医学工程及仪器专业2医学成像系统，高上凯主编，清华大学出版社3G . T.赫尔曼著(严洪编译)，由投影重建图象，科学出版社4董雏申、吴世法、王天童，卷积反投影法从x光图像重建三维轴对称图象， 全国图基科学会议论文集，1989。 G , N. Housfleld Computed Medical Imaging -Journal 0f Computed Assis

21、ted Tomography , 4(5),655 67419806赶荣椿等，数字图像处理导论，百北工业大学出版牡，西安， 1995，p77 1797庄天戈，CT原理与算法），上海交通大学出版杜， 1992, p31 338Radon J ， Uber die Bestimmumg von Funktionen much ihre integralwerte langs gewissserMannnmgfaltigkeiten Berichte Sachsische Akademie derWissenschaften,Leipzig,Math Phys K1 , 1917.69 ：262267（注：专业文档是经验性极强的领域,无法思考和涵盖全面,素材和资料部分来自网络, 供参考。可复制、编制,期待你的好评与关注）