CT反投影滤波重建算法设计

地理与生物信息学院2011/2012学年第二学期实验报告课程名称：医学图像处理和成像技术实验名称：CT反投影滤波重建算法设计班级学号:B090903**学生姓名:****指导教师:****日期：2012年4月2一、实验题目：CT反投影滤波重建算法设计二、实验内容：1.显示图像；2.获得仿真投影数据；3.基于获得的仿真投影数据重建图像。三、实验要求：1.画出Shepp-Logan头模型，显示尺寸为128×128；2.从头模型中获得投影数据，投影数据格式为180×185；3.基于获得的仿真投影数据重建图像，使用R-L卷积函数，重建尺寸为128×128。四、实验过程：1.显示图像：①算法实现流程：I.S-L头模型由10个位置、大小、方向、密度各异的椭圆组成，象征一个脑断层图像。将模型中的椭圆参数写入一个p矩阵中，方便使用其中的数据，并设定所需参数。II.使用循环语句给像素赋值：fori=1:10forx….fory…..判断点(x,y)是否在第i个椭圆内;3如是,则将第i个椭圆折射指数赋给点(x,y);endendendIII.显示仿真头模型，使用imshow(f,[])函数显示出图像。②实验代码：clearall;p=[000.920.69pi/210-0.01840.8740.6624pi/220.2200.310.1172/180*pi0-0.2200.410.16108/180*pi400.350.250.21pi/2500.10.0460.046060-0.10.0460.04607-0.08-0.6050.0460.023080-0.6050.0230.023080.06-0.6050.0460.023pi/28];N=256;x=linspace(-1,1,N);y=linspace(-1,1,N);f=zeros(N,N);fori=1:Nforj=1:Nfork=1:10A=p(k,3);B=p(k,4);x0=p(k,1);y0=p(k,2);x1=(x(i)-x0)*cos(p(k,5))+(y(j)-y0)*sin(p(k,5));y1=-(x(i)-x0)*sin(p(k,5))+(y(j)-y0)*cos(p(k,5));if((x1*x1)/(A*A)+(y1*y1)/(B*B)=1)%判断条件f(i,j)=p(k,6);endendendendf=rot90(f);imshow(f,[])4③运行结果：2.获得仿真投影数据：①算法实现流程：I.取θ∈[00,10,...,1790]，s∈[-92,-91,...,91，92]，即表示在[01790]范围内角度每隔10度取样，每个角度下有185个探测器。II.对于第i个椭圆求出对应θ和s的仿真投影数据：其中，(x0,y0)为中心坐标，A为长轴，B为短轴，a为旋转角度，ρ为折射指数。III.将10个椭圆求出的10个仿真投影数据全部累加起来即可得所要求的仿真投影数据：。②实验代码：clearall;p=[000.920.69pi/210-0.01840.8740.6624pi/220.2200.310.1172/180*pi3-0.2200.410.16108/180*pi400.350.250.21pi/2500.10.0460.046060-0.10.0460.04607-0.08-0.6050.0460.023080-0.6050.0230.023080.06-0.6050.0460.023pi/28];g=zeros(180,185);M=185;s=linspace(-92/64,92/64,M);fori=0:179)(sin)(cos)sincos()(sin)(cos2)(22222002222BAyxsBAABsgi101)()(iisgsg5forj=1:Mfork=1:10x0=p(k,1);y0=p(k,2);A=p(k,3);B=p(k,4);a=p(k,5);%角度b=p(k,6);%折射指数gg=b*2*A*B*sqrt(A*A*(cos(i/180*pi-a))^2+B*B*(sin(i/180*pi-a))^2-(s(j)-x0*cos(i/180*pi)-y0*sin(i/180*pi))^2)/(A*A*(cos(i/180*pi-a))^2+B*B*(sin(i/180*pi-a))^2);g(i+1,j)=g(i+1,j)+gg;endendendg=real(g);g=rot90(g);iptsetpref('ImshowBorder','tight');imshow(g,[])③运行结果：3.基于获得的仿真投影数据重建图像：①算法实现流程：I.进行坐标区间的转换，使用linspace(-92/64,92/64,185)指令。II.选择合适的滤波器，将投影是数据gθ(s)和滤波器进行卷积运算获得g’θ(s)。即，获得标号为(i,j)的像素点对应的灰度值f(i,j)。III.对于不同角度重复上一步骤，直至遍历所有角度。将获得的投影线数据累加起来就得到了重建的点(x0,y0)的灰度值。IV.再利用所得的重建的灰度值画图。②实验代码：clearall;p=[000.920.69pi/210-0.01840.8740.6624pi/220.2200.310.1172/180*pi36-0.2200.410.16108/180*pi400.350.250.21pi/2500.10.0460.046060-0.10.0460.04607-0.08-0.6050.0460.023080-0.6050.0230.023080.06-0.6050.0460.023pi/28];g=zeros(180,185);M=185;s=linspace(-92/64,92/64,M);fori=0:179forj=1:Mfork=1:10x0=p(k,1);y0=p(k,2);A=p(k,3);B=p(k,4);a=p(k,5);%角度b=p(k,6);%折射指数gg=b*2*A*B*sqrt(A*A*(cos(i/180*pi-a))^2+B*B*(sin(i/180*pi-a))^2-(s(j)-x0*cos(i/180*pi)-y0*sin(i/180*pi))^2)/(A*A*(cos(i/180*pi-a))^2+B*B*(sin(i/180*pi-a))^2);g(i+1,j)=g(i+1,j)+gg;endendendd=1.4531/92;h=zeros(1,370);forn=1:370m=n-185;ifm==0h(n)=1/(4*d*d);elseifmod(m,2)==0h(n)=0;elseh(n)=-1/(m*m*pi*pi*d*d);endendfori=1:180g1(i,:)=conv(g(i,:),h);endf1=zeros(128,128);fori=1:128forj=1:1287x=(i-64)*d;y=(j-64)*d;fora=1:180o=a/180*pi;sx=x*cos(o)+y*sin(o);nx=sx/d;ni=floor(nx);k=ni+277;f1(i,j)=f1(i,j)+g1(a,k)*(nx-ni)+g1(a,k+1)*(ni+1-nx);endendendf1=real(f1);F1=rot90(f1);imshow(F1,[])③运行结果：实验小结：1.实验中使用S-L头模型好处之一，就是可以很简单的获得该模型的投影数据。模型由10个位置、大小、方向、密度各异的椭圆组成，所以需要一个p(10,6)的矩阵来存储数据。利用来判断点(x,y)是否是在第i个椭圆内。2.实验内容二中，仿真投影数据的获得，则只需要用一个g(s,θ)的判定式就可以获得，但是要注意的就是由于投影数据格式为：180×185，所以s的取值不仅仅只是简单的linespace(-1,1,M)，而是s=linspace(-92/64,92/64,M)。3.实验内容三中，首先用到了实验二中所获得的投影数据，对数据进行卷积反投影的重建。由gθ(s)和滤波器进行卷积运算获得g’θ(s)，进行插值，经过多次卷积运算，对数据进行累加就可以得到重建后的灰度值，进行画图了。12222ByAx84.本次实验的关键在于函数表达式的正确输入以及对矩阵元素下表的熟练运用。在表达式的输入过程中，认识到了“*”和“.*”的区别，在程序运行调试时，了解到了“，”和“；“的区别。在函数的使用上，认识到了MATLAB的强大和Help的妙用，这些都对以后的实验打下了一定的基础。