2022年CT图像重建 2.pdf

-1- 昆明理工大学信息工程与自动化学院学生实验报告（ 20092010 学年 第 一 学期 ）课程名称：医学成像系统与放射治疗装置开课实验室： 3208 2008 年 12 月 24 日年级、专业、班学号姓名成绩实验项目名称CT图像重建指导教师刘利军教师评语教师签名：年月日一、实验目的与意义医学成像技术是生物医学工程专业的一门重要的专业课程,课程主要涉及X 光仪器， CT 仪器， MRI仪器和核医学仪器的工作原理及成像方法。其中 CT算法的出现又为后来数字化医学成像技术的发展提供了基础。该门课程为生物医学工程专业的专业基础课。CT技术是医学成像系统中的一种重要手段。它通过特定的算法，利用计算机的高速运算功能，可以在短时间内快速呈现人体断层图像。让学生练习CT图像的重建有助于学生理解CT算法的内容，熟悉数字图像重建的过程。同时也能培养学生的团队精神和解决实际问题的能力。二、实验算法原理1、MATLAB 处理数字图像的基本函数；2、X-CT三维图像重建的基本算法。CT图象重建有四种基本的算法: 矩阵法 , 迭代法 , 傅立叶算法 , 反投影算法 . 我们采用的方法为卷积反投影 . 卷积反投影有: 平行光束投影的卷积反投影算法, 等角扇形光来投影的重建算法. 1）. 平行光束投影的卷积反投影算法从投影重建三维物体的图像，就是重建一个个横断面。这样三堆图像的重建就归结为二维图象的重建。二维图像的重建问题可以从数学上描述如下。假定),(yxg表示一个二维的未知函数，通过),(yxg的直线称为光钱( 见图 21) 。沿光线),(yxg的积分称作光线积分。沿相同方向的一组光线积分，就构成一个投影。图21中垂直于直线CC ( 与X轴夹角为) 的光线所形成。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 12 页 - - - - - - - - - -2- 图2.1 ),(yxg在方向的投影)(tP的投影)(tP，称之为),(yxg在方向的投影。光线积分和投影在数学上可以定义如下：在图 2.1 中直线 AB 的方程为：1sincostYX (2.1) 其中1t是AB 到原点的距离，),(yxg沿AB 的积分为：dxdytyxyxgdsyxgtPAB)sincos(),(),()(11 (2.2) 对于给定的，),(yxg在方向的投影)(tP是t的函数。 如果),(yxg在各个方向的投影已知，),(yxg就可以唯一确定。下面就讨论卷积反投影重建算法。假定投影方向，如图 2.2 ，将坐标),(yx旋转角（逆时针方向） 形成坐标系),( st。),(yxg在),(st坐标系中为),(stg。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 12 页 - - - - - - - - - -3- 图2.2 傅立叶切片定理示意图坐标系),( st与),(yx之间的关系为：yxstcossinsincos（2.3 ）显然dsstgtP),(（2.4 ）令)(wS为)(tP的傅立叶变换则dtwtjtPwS)2exp()()(dsdtwtjstg)2exp(),(（2.5 ）将上式变换到),(yx坐标系中，注意到变换的可比行列式1cossinsincosysytxstt（2.6 ）从而得到：dxdyyxjyxgwS)sincos(2exp),()(dxdyvyuxjyxg)(2exp),(（2.7 ）其中sincosvu（2.8 ）名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 12 页 - - - - - - - - - -4- 若令),(yxg的傅立叶变换为),(vuG，由（ 2.8 ）可知),(),()(GvuGS（2.9 ）若),(yxg的傅立叶变换为),(vuG的极坐标表示。这说明),(yxg在方向的投影)(tP傅立叶变换)(wS等于),(vuG在与 u轴成角的直线上的值。这就是著名的傅立叶投影切片定理。可见在整个),(vu平面),(vuG可以利用各个方向的投影来得到，从而),(yxg也可以通过求),(vuG的傅立叶反交换的办法求得：dudvvyuxjvuGyxg)(2exp),(),(（2.10 ）变换到极坐标中sincosvu，得到ddyxjGyxg)sincos(2exp),(),(200（2.11 ）经推导得0)2exp()(),(ddtjSyxg（2.12 ）其中sincosyxt若令dtjStQ)2exp()()(（2.13 ）则0)sincos(),(dyxQyxg (2.14) (2.13) 式右端是两频谱函数)(wS和)(H的乘积的傅立叶反变换。)(wS是投影)(tP傅立叶变换。若)(H的傅立叶反变换为)(th，则根据卷积定理有：dthPtQ)()()(（2.15 ）或)()()(thtPtQ其中dtjth)2exp()(（2.16 ）当图像的频谱是有限带宽时，则上式变为名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 12 页 - - - - - - - - - -5- dtjth)2exp()(00（2.17 ）由于图象及其频谱都是离散采样的，假定图象采样间隔为， 则根据采样定理2/10。为了进行数学处理，只需知道h (t)在有限带宽上的离散采样点的值这样我们有2222/10)4/(1)(nnh（2.18 ）其中 n为正负整数。（2.18 ）的离散形式为mmnhmPnQ)()()(（2.19 ）假定)(mP在1, 1 ,0Nm之外的值为 0，则上式变为10)()()(NmmnhmPnQ（2.20 ）或1)1()()(NNmmhmnPnQ（2.21 ）其中1,2 ,1 ,0Nn从而可见为确定)(tP的N个采样点上的)(nQ的值， 需要使用)(nh的2N1个点上的值，从n=一(N 1) 到(N 1) 。为求得)(nQ，利用傅立叶变换计算卷积是比较快的方法，为清除循环卷积的周期交叠效应，实际上)(nh取2N 个点，)(mP补0，使之有 (2N1) 个元素，则)(tP在N个采样点上就避免了交叠，如果使用以 2为基的 FFT(快速傅立叶变换) 算法，)(mP和)(nh都必需朴 0至(2N一1)个元素， (2N一1)为大于等于 2N l 的最小的 2的整数幂。计算)(nQ的过程可以写为0)(0)()(nhFFTnPFFTFFTnQ（2.22 ）其中 FFT和IFFT分别表示快速傅立叶变换和反变换，光滑窗是在滤波过程中加入的光滑因子，例如引用汉明窗，有时可以改进重建效果。对于各个方向的投影，得到)(nQ之后就可以由(2.22) 来计算),(yxg。重建步骤可以归纳为：第一步：卷积，也称滤波，由(222) 对每个方向计算)(nQ。第二步。反投影，由(2 14) 的近似形式名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 12 页 - - - - - - - - - -6- MiiiiyxQMyxg1)sincos(),(（2.23 ）来计算),(yxg的近似值),(?yxg。M 为投影个数i为投影方向角，他们均匀的分布在0的范围内。当计算)sincos(iiyxQi时，iiyxtsincos，不一定在)(nQ的整离散点上，这就需要插值求得，预先将)(nQ插值加密，即最靠近的点，可以提高计算速度。2). 等角扇形光来投影的重建算法几乎所有的快遗CT设备都是用的扇形光束。这里叙述的是等角度光束投影，如图 2.3 ，测量投影数据的探测器等间距地分布在1D2D弧上，弧的半径为2D， D为光源到图像中心的距离。在下文中，),(rf图象在极坐标中的表示。)(R表示在方向角为的投影中位量角为的光线产生的投影数据。通过中心的光线其为0。L表示从光源到像素),(r的距离。图2.3 等扇形束投影重建算法中的变量)sin(2),(22DrrDL（2.24 ）表示在方向角为的投影中通过像素),(r的光线的位置角)sin()cos(tantan),(11DESEP（2.25 ）图像),(rf和扇形投影)(R有下述关系dhrRDLrf)()sin(21)(cos1),(202（2.26 ）名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 12 页 - - - - - - - - - -7- 其中r和是投影中光线的最大位置角，从而可以得到这种重建算法的执行步骤：第一步：投影的修改，假定投影的抽样间隔为，抽样数据)(nR通过下式进行修改，nDnRnRcos)()(1（2.27 ）第二步：卷积（滤波）将第一步修改了的投影与响应函数)(g进行卷积)()sin1(21)(2hg（2.28 ）dmngnRQ)()()(（2.29 ）其离散形式为：)()()(1mngmRnQMm（2.30 ）这里也和上节一样可以加进一光滑窗，改进重建效应。第三步：反投影，就是执行（2.30 ）的积分dQLrf)(),(1),(202（2.31 ）近似的有：)(),(12),(12QLmyxfMi（2.32 ）其),(yxf是图象),(yxf的近似图像，cosrx，sinry， M是投影数，假定投影均匀地分布在02内。为求得滤波后的投影Q，对象素),(yx的贡献，首先由(2.24) 和(2.25) 计算出L和所有Q对),(yx的贡献加起来再乘以M/2就得),(yxf。四、实验程序代码原图象代码 : 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;. 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 12 页 - - - - - - - - - -8- 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(A2*cos(theta-alpha)2+B2*sin(theta-alpha)2-(R-x0*cos(theta)-y0*sin(theta)2); NN=A2*cos(theta-alpha)2+B2*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(A2*cos(theta-alpha)2+B2*sin(theta-alpha)2-(R-x0*cos(theta)-y0*sin(theta)2); NN=A2*cos(theta-alpha)2+B2*sin(theta-alpha)2; g(i,forw)=rho*2*A*B*MM/NN; R=-R; back=back-1; MM=sqrt(A2*cos(theta-alpha)2+B2*sin(theta-alpha)2-(R-x0*cos(theta)-y0*sin(theta)2); NN=A2*cos(theta-alpha)2+B2*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*vstep2); 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; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 12 页 - - - - - - - - - -9- 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/(pi2*vstep2*(4*n2-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 n00&n00&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) 五、实验结果与分析反投影一般步骤为：程序流程 : 原像取投影反投影重建重建后图像名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 12 页 - - - - - - - - - -11- 卷积反投影法结果：提取图象数据初始化产生R-L函数产生S-L函数进行radon反变换显示重建图象结果结束名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 12 页 - - - - - - - - - -12- 重建图象原图象六、心得体会通过本次的实验, 对 CT 图象重建的基本方法之一: 卷积反投影 , 有了进一步的认识, 在实验的过程中,采用的图象是经典的Shep-Logen 头模型 , 得到的结果与原图象相比, 有一定的差异, 但影响不大 . 有待进一步的改进算法. 七、参考文献1 医用电子仪器及装备医学成像系统及放射治疗装置类，唐庆玉主编， 清华大学电机系生物医学工程及仪器专业2 医学成像系统 ，高上凯主编，清华大学出版社3 G T.赫尔曼著 ( 严洪编译 )，由投影重建图象，科学出版社4 董雏申、吴世法、 王天童，卷积反投影法从x光图像重建三维轴对称图象，全国图基科学会议论文集，l989 。5 G ， NHousfleld Computed Medical Imaging -Journal 0f Computed Assisted Tomography，4(5) ，655 6741980 6 赶荣椿等，数字图像处理导论，百北工业大学出版牡，西安，1995，p77179 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,MathPhysK1，1917.69 ：262267名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 12 页，共 12 页 - - - - - - - - -