CT反投影滤波重建算法设计(汇编)

地理与生物信息学院

2011 / 2012 学年第二学期

实验报告

课程名称：医学图像处理和成像技术实验名称：CT反投影滤波重建算法设计

班级学号: B090903**

学生姓名: ****

指导教师: ****

日期：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个椭圆内;

如是,则将第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,[])

③运行结果：

2. 获得仿真投影数据：

①算法实现流程：

I. 取θ∈ [00, 10, ..., 1790]， s ∈[-92, -91, ..., 91，92]，即表示在[0 1790]范围内角度每隔10度取样，每个角度下有185个探测器。

II. 对于第i 个椭圆求出对应θ和s 的仿真投影数据：

其中，(x 0, y 0)为中心坐标，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

)

(sin )(cos )sin cos ()(sin )(cos 2)(22222

002222αθαθθθαθαθρθ-+-----+-=B A y x s B A AB s g i ∑==10

1)()(i i s g s g θθ

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