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

地理与生物信息学院2011 / 2012 学年第 二 学期实 验 报 告 课程名称： 医学图像处理和成像技术 实验名称：CT反投影滤波重建算法设计 班 级 学 号: B090903* 学 生 姓 名: * 指 导 教 师: * 日 期： 2012 年 4 月1、 实验题目：CT反投影滤波重建算法设计二、实验内容：1.显示图像；2.获得仿真投影数据；3.基于获得的仿真投影数据重建图像。3、 实验要求：1. 画出Shepp-Logan头模型，显示尺寸为128128；2. 从头模型中获得投影数据，投影数据格式为180185；3. 基于获得的仿真投影数据重建图像，使用R-L卷积函数，重建 尺寸为128128。4、 实验过程：1. 显示图像：算法实现流程：I. S-L头模型由10个位置、大小、方向、密度各异的椭圆组成，象征一个脑断层图像。将模型中的椭圆参数写入一个p矩阵中，方便使用其中的数据，并设定所需参数。II. 使用循环语句给像素赋值：for i=1:10 for x. for y. 判断点(x, y)是否在第i个椭圆内; 如是,则将第i个椭圆折射指数赋给点(x, y); end end endIII. 显示仿真头模型，使用imshow(f,)函数显示出图像。实验代码：clear all; p=0 0 0.92 0.69 pi/2 1 0 -0.0184 0.874 0.6624 pi/2 2 0.22 0 0.31 0.11 72/180*pi 0 -0.22 0 0.41 0.16 108/180*pi 4 0 0.35 0.25 0.21 pi/2 5 0 0.1 0.046 0.046 0 6 0 -0.1 0.046 0.046 0 7 -0.08 -0.605 0.046 0.023 0 8 0 -0.605 0.023 0.023 0 8 0.06 -0.605 0.046 0.023 pi/2 8; N=256; x=linspace(-1,1,N); y=linspace(-1,1,N); f=zeros(N,N); for i=1:N for j=1:N for k=1:10 A=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); end end end end f=rot90(f); imshow(f,)运行结果：2. 获得仿真投影数据：算法实现流程：I. 取 00, 10, ., 1790， s-92, -91, ., 91，92，即表示在0 1790范围内角度每隔10度取样，每个角度下有185个探测器。II. 对于第i个椭圆求出对应和s的仿真投影数据：其中，(x0, y0)为中心坐标，A为长轴，B为短轴，a为旋转角度，为折射指数。III. 将10个椭圆求出的10个仿真投影数据全部累加起来即可得所要求的仿真投影数据： 。 实验代码：clear all;p=0 0 0.92 0.69 pi/2 1 0 -0.0184 0.874 0.6624 pi/2 2 0.22 0 0.31 0.11 72/180*pi 3 -0.22 0 0.41 0.16 108/180*pi 4 0 0.35 0.25 0.21 pi/2 5 0 0.1 0.046 0.046 0 6 0 -0.1 0.046 0.046 0 7 -0.08 -0.605 0.046 0.023 0 8 0 -0.605 0.023 0.023 0 8 0.06 -0.605 0.046 0.023 pi/2 8; g=zeros(180,185); M=185; s=linspace(-92/64,92/64,M); for i=0:179 for j=1:M for k=1:10 x0=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; end end end g=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. 再利用所得的重建的灰度值画图。实验代码：clear all;p=0 0 0.92 0.69 pi/2 1 0 -0.0184 0.874 0.6624 pi/2 2 0.22 0 0.31 0.11 72/180*pi 3 -0.22 0 0.41 0.16 108/180*pi 4 0 0.35 0.25 0.21 pi/2 5 0 0.1 0.046 0.046 0 6 0 -0.1 0.046 0.046 0 7 -0.08 -0.605 0.046 0.023 0 8 0 -0.605 0.023 0.023 0 8 0.06 -0.605 0.046 0.023 pi/2 8; g=zeros(180,185); M=185; s=linspace(-92/64,92/64,M); for i=0:179 for j=1:M for k=1:10 x0=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; end end endd=1.4531/92;h=zeros(1,370);for n=1:370 m=n-185; if m=0 h(n)=1/(4*d*d); elseif mod(m,2)=0 h(n)=0; else h(n)=-1/(m*m*pi*pi*d*d); endendfor i=1:180 g1(i,:)=conv(g(i,:),h);endf1=zeros(128,128);for i=1:128 for j=1:128 x=(i-64)*d; y=(j-64)*d; for a=1:180 o=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); end endendf1=real(f1);F1=rot90(f1);imshow(F1, ) 运行结果：实验小结：1.实验中使用S-L头模型好处之一，就是可以很简单的获得该模型的投影数据。模型由10个位置、大小、方向、密度各异的椭圆组成，所以需要一个p(10,6)的矩阵来存储数据。利用 来判断点(x,y)是否是在第i个椭圆内。2.实验内容二中，仿真投影数据的获得，则只需要用一个g(s,)的判定式就可以获得，但是要注意的就是由于投影数据格式为：180185，所以s的取值不仅仅只是简单的linespace (-1,1,M)，而是s=linspace(-92/64,92/64,M)。3.实验内容三中，首先用到了实验二中所获得的投影数据，对数据进行卷积反投影的重建。由g(s)和滤波器