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

地理与生物信息学院

2012 / 2013 学年第二学期

实验报告

课程名称：医学图像处理和成像技术

实验名称：CT反投影滤波重建算法设计

班级学号: B10090405

学生姓名: 陈洁

指导教师: 戴修斌

日期：2013 年 5 月

一、实验题目：CT反投影滤波重建算法设计

二、实验内容：

1.显示图像；

2.获得仿真投影数据；

3.基于获得的仿真投影数据重建图像。

三、实验要求：

1.Shepp-Logan头模型：

画出Shepp-Logan头模型，简称S-L模型，头模型尺寸设定为128×128；

2.仿真投影数据的获得：

从头模型中获得投影数据，投影数据格式为180×185，即[0，179°]范围内角度每隔1°取样，每个角度下有185个探测器；

3.卷积反投影重建算法的实现：

基于获得的仿真投影数据重建图像，使用R-L卷积函数，重建尺寸为128×128。

四、实验过程：

实验1. Shepp-Logan头模型

①算法实现流程：

I. S-L头模型由10个位置、大小、方向、密度各异的椭圆组成，

象征一个脑断层图像。

Shepp-Logan头模型中的椭圆参数：

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]；

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

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

)

(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 θθ