地理与生物信息学院
2011 / 2012 学年第 二 学期
实 验 报 告
课程名称： 医学图像处理和成像技术 实验名称：CT反投影滤波重建算法设计
班 级 学 号: B******** 学 生 姓 名: **** 指 导 教 师: ****
日 期： 2012 年 4 月
一、实验题目：CT反投影滤波重建算法设计 二、实验内容：
1.显示图像；
2.获得仿真投影数据；
3.基于获得的仿真投影数据重建图像。
三、实验要求：
1. 画出Shepp-Logan头模型，显示尺寸为128×128；
2. 从头模型中获得投影数据，投影数据格式为180×185； 3. 基于获得的仿真投影数据重建图像，使用R-L卷积函数，重建 尺寸为128×128。
四、实验过程：
1. 显示图像： ①算法实现流程：
I. S-L头模型由10个位置、大小、方向、密度各异的椭圆组成，象征一个脑断层图像。将模型中的椭圆参数写入一个p矩阵中，方便使用其中的数据，并设定所需参数。
II. 使用循环语句给像素赋值： for i=1:10 for x….
for y…..
判断点(x, y)是否在第i个椭圆内;
2
如是,则将第i个椭圆折射指数赋给点(x, y); end end end
III. 显示仿真头模型，使用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,[])
3
③运行结果：
2. 获得仿真投影数据： ①算法实现流程：
I. 取θ∈ [00, 10, ..., 1790]， s∈[-92, -91, ..., 91，92]，即表示在[0 1790]范围内角度每隔10度取样，每个角度下有185个探测器。
II. 对于第i个椭圆求出对应θ和s的仿真投影数据：
g(s)i2ABA2cos2()B2sin2()(sx0cosy0sin)2A2cos2()B2sin2()其中，(x0, y0)为中心坐标，A为长轴，B为短轴，a为旋转角度，ρ为折射指数。
III. 将10个椭圆求出的10个仿真投影数据全部累加起来即可得所要求的仿真投影数据：g ( s) g i (s ) 。
i1②实验代码：
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/,92/,M); for i=0:179
10 4
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/,92/,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
5
-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/,92/,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
d=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); end end
for i=1:180
g1(i,:)=conv(g(i,:),h); end
f1=zeros(128,128); for i=1:128 for j=1:128
6
x=(i-)*d; y=(j-)*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 end end
f1=real(f1); F1=rot90(f1); imshow(F1,[ ])
③运行结果：
实验小结：
1.实验中使用S-L头模型好处之一，就是可以很简单的获得该模型的投影数据。模型由10个位置、大小、方向、密度各异的椭圆组x2y2成，所以需要一个p(10,6)的矩阵来存储数据。利用 2 2 1 来判
AB断点(x,y)是否是在第i个椭圆内。
2.实验内容二中，仿真投影数据的获得，则只需要用一个g(s,θ)的判定式就可以获得，但是要注意的就是由于投影数据格式为：180×185，所以s的取值不仅仅只是简单的linespace (-1,1,M)，而是s=linspace(-92/,92/,M)。
3.实验内容三中，首先用到了实验二中所获得的投影数据，对数据进行卷积反投影的重建。由gθ(s)和滤波器进行卷积运算获得g’θ(s)，进行插值，经过多次卷积运算，对数据进行累加就可以得到重建后的灰度值，进行画图了。
7
4.本次实验的关键在于函数表达式的正确输入以及对矩阵元素下表的熟练运用。在表达式的输入过程中，认识到了“*”和“.*”的区别，在程序运行调试时，了解到了“，”和“；“的区别。在函数的使用上，认识到了MATLAB的强大和Help的妙用，这些都对以后的实验打下了一定的基础。
8
