CT图像重建

-1-昆明理工大学信息工程与自动化学院学生实验报告（2009—2010学年第一学期）课程名称：医学成像系统与放射治疗装置开课实验室：32082008年12月24日年级、专业、班学号姓名成绩实验项目名称CT图像重建指导教师刘利军教师评语教师签名：年月日一、实验目的与意义医学成像技术是生物医学工程专业的一门重要的专业课程,课程主要涉及X光仪器，CT仪器，MRI仪器和核医学仪器的工作原理及成像方法。其中CT算法的出现又为后来数字化医学成像技术的发展提供了基础。该门课程为生物医学工程专业的专业基础课。CT技术是医学成像系统中的一种重要手段。它通过特定的算法，利用计算机的高速运算功能，可以在短时间内快速呈现人体断层图像。让学生练习CT图像的重建有助于学生理解CT算法的内容，熟悉数字图像重建的过程。同时也能培养学生的团队精神和解决实际问题的能力。二、实验算法原理1、MATLAB处理数字图像的基本函数；2、X-CT三维图像重建的基本算法。CT图象重建有四种基本的算法:矩阵法,迭代法,傅立叶算法,反投影算法.我们采用的方法为卷积反投影.卷积反投影有:平行光束投影的卷积反投影算法,等角扇形光来投影的重建算法.1）.平行光束投影的卷积反投影算法从投影重建三维物体的图像，就是重建一个个横断面。这样三堆图像的重建就归结为二维图象的重建。二维图像的重建问题可以从数学上描述如下。假定),(yxg表示一个二维的未知函数，通过),(yxg的直线称为光钱(见图2．1)。沿光线),(yxg的积分称作光线积分。沿相同方向的一组光线积分，就构成一个投影。图2．1中垂直于直线'CC(与X轴夹角为)的光线所形成。-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。-3-图2.2傅立叶切片定理示意图坐标系),(st与),(yx之间的关系为：yxstcossinsincos（2.3）显然dsstgtP),(（2.4）令)(wS为)(tP的傅立叶变换则dtwtjtPwS)2exp()()(dsdtwtjstg)2exp(),(（2.5）将上式变换到),(yx坐标系中，注意到变换的可比行列式1cossinsincosysytxstt（2.6）从而得到：dxdyyxjyxgwS)]sincos(2exp[),()(dxdyvyuxjyxg)](2exp[),(（2.7）其中sincosvu（2.8）-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）当图像的频谱是有限带宽时，则上式变为-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,0Nm之外的值为0，则上式变为10)()()(NmmnhmPnQ（2.20）或1)1()(][)(NNmmhmnPnQ（2.21）其中1,2,1,0Nn从而可见为确定)(tP的N个采样点上的)(nQ的值，需要使用)(nh的2N—1个点上的值，从n=一(N—1)到(N—1)。为求得)(nQ，利用傅立叶变换计算卷积是比较快的方法，为清除循环卷积的周期交叠效应，实际上)(nh取2N个点，)(mP补0，使之有(2N—1)个元素，则)(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。重建步骤可以归纳为：第一步：卷积，也称滤波，由(2．22)对每个方向计算)(nQ。第二步。反投影，由(2．14)的近似形式-6-MiiiiyxQMyxg1)sincos(),(（2.23）来计算),(yxg的近似值),(ˆyxg。M为投影个数i为投影方向角，他们均匀的分布在0~的范围内。当计算)sincos(iiyxQi时，iiyxtsincos，不一定在)(nQ的整离散点上，这就需要插值求得，预先将)(nQ插值加密，即最靠近的点，可以提高计算速度。2).等角扇形光来投影的重建算法几乎所有的快遗CT设备都是用的扇形光束。这里叙述的是等角度光束投影，如图2.3，测量投影数据的探测器等间距地分布在1D2D弧上，弧的半径为2D，D为光源到图像中心的距离。在下文中，),(rf图象在极坐标中的表示。)(R表示在方向角为的投影中位量角为的光线产生的投影数据。通过中心的光线其为0。L表示从光源到像素),(r的距离。图2.3等扇形束投影重建算法中的变量)sin(2),,(22DrrDL（2.24）表示在方向角为的投影中通过像素),(r的光线的位置角)sin()cos(tantan),,(11'DESEP（2.25）图像),(rf和扇形投影)(R有下述关系dhrRDLrf)())sin((21)(cos1),('''202（2.26）-7-其中r和是投影中光线的最大位置角，从而可以得到这种重建算法的执行步骤：第一步：投影的修改，假定投影的抽样间隔为，抽样数据)(nR通过下式进行修改，nDnRnRcos)()(1'（2.27）第二步：卷积（滤波）将第一步修改了的投影与响应函数)(g进行卷积)()sin1(21)(2'hg（2.28）dmngnRQ)()()(（2.29）其离散形式为：)()()(1'mngmRnQMm（2.30）这里也和上节一样可以加进一光滑窗，改进重建效应。第三步：反投影，就是执行（2.30）的积分dQLrf)(),,(1),(202（2.31）近似的有：)(),,(12),(12QLmyxfMi（2.32）其),(yxf是图象),(yxf的近似图像，cosrx，sinry，M是投影数，假定投影均匀地分布在0~2内。为求得滤波后的投影Q，对象素),(yx的贡献，首先由(2.24)和(2.25)计算出L和'．所有Q对),(yx的贡献加起来再乘以M/2就得),(yxf。四、实验程序代码原图象代码:p=phantom(256);imshow(p)卷积反投影法代码：H=[000.920.6990*pi/1801;...0-0.01840.8740.662490*pi/180-0.98;...0.2200.310.1172*pi/180-0.2;...-0.2200.410.16108*pi/180-0.2;...00.350.250.2190*pi/1800.1;...00.10.0460.04600.2;...0-0.10.0460.04600.2;...-0.08-0.6050.0460.02300.1;...-8-0-0.6050.0230.02300.1;...0.06-0.6050.0460.02390*pi/1800.1];angle=1;N=100;vstep=2/N;ax=N/2+1;fork=1:(180/angle+1)theta=(k-1)*angle*pi/180;fori=1:10x0=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;forj=1:N/2R=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)=rh