CT仿真实验报告

西安交通大学实验报告
共 7 页课程医学成像实验
系别生物医学工程实验日期 2012 年 12 月 XX日专业班级医电 01 班组别交报告日期 2013 年 01 月 02日姓名学号报告退发 (订正、重做)
同组者教师审批签字
实验名称 CT重建原理——投影数据采集实验1
一、实验目的以及要求
实验目的：利用CTSim模拟软件生成投影数据，为滤波反投影重
建实验做准备。

实验基本要求：用CTSim程序完成实验模拟，分析评价结果。

二、实验内容
1、利用CTSim模拟软件生成椭圆的平行束投影数据；
2、利用CTSim模拟软件生成Shepp-Logan图的平行束投影数据；
3、对生成的投影数据进行初步评价。

三、实验步骤
A、完成CTSim模拟软件生成椭圆的平行束投影数据；
1、点击软件ctsim，打开软件界面，点击File，选择creat phantom，
选择Herman Head，得到椭圆的灰度图像，如图：
图1 软件界面图2 选择界面
图3 椭圆的原始数据
2、在选择椭圆窗口的情况下，点击Process，选择rasterize，点击OK，将图像进行光栅化，如图：
图4 光栅化参数图5 光栅化后的图像
3、在选择unnamed3窗口，选择View，选择Auto Scale Parameters，
将Standard Deviation Factor参数改为，点击OK，得到处理后图像，如图：
图6 光栅化参数图7 参数优化后图像
4、回到herman窗口，点击Process，选择Projection Paramaters，参数默认即可，点击ok，开始采集数据，如图：
图8 参数选择界面
5、得到投影参数后，在选择unnamed4窗口界面下，选择Analyze，
在选择Plot Histogram，得到平行束投影分析数据，如图：
图9 平行束投影后的数据图10 投影数据分析图
6、在选择平行束投影后数据窗口情况下，选择Reconstruct，选择Filtered Backprojection Parameters，选择默认参数即可，点击OK，得到重建数据，如图：
图11 参数选择界面图12 重建后的图像
7、可选择不同的View参数对重建后的图像进行参数的优化调整，得到最优的观察效果。

B、完成CTSim模拟软件生成Shepp-Logan图的平行束投影数据；
1、由于生成Shepp-Logan图的平行束投影数据和前面A中的步骤基本类似，只是选择的模拟图像不同而已，故在此不再一一列出步骤。

四、实验结果及分析
A、椭圆平行束投影
图13 椭圆平行束投影及其重建图像
图13 采用Projection方式图14 时域的采样参数时的数据采集过程
图15 椭圆投影数据的分析直方图
B、Shepp-Logan图的平行束投影
图16 Shepp-Logan图平行束投影及其重建图像
图17 数据采集过程图18 时域的采样参数
图19 Shepp-Logan图投影数据的分析直方图结果分析：
由图15和图19的投影数据的分析可知，它们的统计直方分布图还是有明显区别的，这是因为Shepp-Logan的图像更复杂，厚度更大，因此，X射线通过后的衰减程度也是不一样的。

在椭圆的统计直方图数据中，只有在处，有一个峰值，而在Shepp-Logan中，在到之间有两个峰值。

同时，从重建中的图像中可以看出，与原来的图像相比较，还是存在一定的噪声干扰的，在后期的数据处理中应该加入一定滤波技术，使图像更加清楚地展示。

西安交通大学实验报告
成绩
共 17 页课程医学成像实验
系别生物医学工程实验日期 2012 年 12 月 XX日专业班级医电 01 班组别交报告日期 2013 年 01 月 02日姓名谭礼茂学号报告退发 (订正、重做)
同组者教师审批签字
实验名称 CT重建原理——投影数据采集实验2
一、实验目的
对仿真数据与实际图像数据实现滤波反投影重建，加深对CT成像过程及原理的理解。

二、实验要求
用CTSim程序完成滤波反投影重建实验，分析评价结果。

三、实验内容
1、对人脑体模仿真数据与实际的CT图像数据进行滤波反投影重
建实验，比较直接反投影与R——L和S——L滤波方法的重建结果；
2、分析不同视角条件下各种重建的结果；
3、分析噪声对各种滤波反投影重建方法的影响，比较各种滤波反
投影重建方法的抗噪能力；
4、计算重建误差。

四、实验步骤
在第一次实验时，已经获得了图像的投影方式，并在实验报告中详细列出了操作的步骤，故在此不再重复。

1、在仿真软件中，有很多可以选择的选项，比如在扫描数据时，
选择Process，选择Projection Parameters，在这个对话框里面，Geometry选项是用来选择投影的方式的，一共有三种选择的方式，如图：
图1 投影方式的选择
2、同样在图像重建时，可以选择不同的滤波方式，重建的图像进行滤波，选择Reconstruct，选择Filted Backprojection，在对话框中有Filter选项，其中包含四种滤波方式可供选择，如图所示：
图2 滤波方式的选择
3、重建误差的计量可以使用软件自带的相减，在Image选项下面，点击subtract，可以获得重建后的图像与原始图像之间的差值，如图：
图3 图像之间的相减
五、实验结果
A、椭圆数据
1、首先，利用模拟软件进行椭圆的不同投影方式下的数据采集：
图4 采用软件自带的三种投影方式：平行束投影、等角投影和等线投影，得到的K空间的图像
2、三种不同的投影方式在同一种滤波方式：Bandlimit滤波，下
得到的重建图像以及重建误差，如图：
图5 投影数据（第一行）、重建图像（第二行）、重建误差（第三行）
3、三种不同的投影方式在同一种滤波方式：hamming滤波，下得到的重建图像以及重建误差，如图：
图6 投影数据（第一行）、重建图像（第二行）、重建误差（第三行）
4、三种不同的投影方式在同一种滤波方式：hanning滤波，下得到的重建图像以及重建误差，如图：
图7 重建图像（第一行）、重建误差（第二行）
5、三种不同的投影方式在同一种滤波方式：Cosine滤波，下得到的重建图像以及重建误差，如图：
图8 重建图像（第一行）、重建误差（第二行）
B、Shepp-Logan图
1、首先，利用模拟软件进行的不同投影方式下的数据采集：
图9图4 采用软件自带的三种投影方式：平行束投影、等角投影和等线投影，得到的K空间的图像
2、三种不同的投影方式在同一种滤波方式：Bandlimit滤波，下得到的重建图像以及重建误差，如图：
图10 重建图像（第一行）、重建误差（第二行）
到的重建图像以及重建误差，如图：
图11 重建图像（第一行）、重建误差（第二行）
4、三种不同的投影方式在同一种滤波方式：hanning滤波，下得到的重建图像以及重建误差，如图：
图12 重建图像（第一行）、重建误差（第二行）
的重建图像以及重建误差，如图：
图13 重建图像（第一行）、重建误差（第二行）
C 、在MATLAB 平台下，调用phantom 函数，获得标准图像，然后进行投影，并在无滤波、Ram-Lak 滤波重建、Shepp-Logan 滤波重建，（代码见附录程序1）得到的结果如下图所示：
原始图
像
Ram-Lak 滤波重建
50100150200250
无滤波重建
50100150200250
Shepp-Logan 滤波重建
50100150200250
图14 原始图像以及三种滤波情况下的重建图像
D、利用MATLAB平台，选用实际中CT脑部图像，然后进行投影，并在无滤波、Ram-Lak滤波重建、Shepp-Logan滤波重建，（代码见附录程序2）得到的结果如下图所示：
原始图像
Ram-Lak滤波重建
100
200
300
400
500
无滤波重建
100
200
300
400
500
Shepp-Logan滤波重建
100
200
300
400
500
图15 原始图像以及三种滤波情况下的重建图像
E、分析不同的视角条件下，投影并进行重建的图像情况，本次采用两种情况下的视角
a、从0度到180度，间隔为10度，共18个视角；
b、从0度到180度，间隔为2度，共90个视角。

同样分别采用三种滤波方式进行重建（代码见附录程序3）。

结果如下图所示：
原始图像
无滤波重建
50100150200250
Ram-Lak 滤波重建
50100150200250
Shepp-Logan 滤波重建
50100150200250
图16 a 情况下的重建图像
原始图像
无滤波重建
50100150200250
Ram-Lak 滤波重建
50100150200250
Shepp-Logan 滤波重建
50100150200250
图17 b 情况下的重建图像
F 、噪声测试，在实际测试中添加高斯白噪声（均值为0，方差为），然后进行投影之后，滤波重建，这三种情况下的滤波效果均不理想（代码见附录程序4），然后，再选择添加Poisson 噪声，然后投影并在三种滤波情况下进行重建图像，并采用模拟图像和实际的CT 图像分别进行测试（代码见附录程序5），结果如下图所示：
original
无滤波重建
50
100
150
200
250
50100150200250
Ram-Lak 滤波重建
50
100
150
200
250
50100150200250
Shepp-Logan 滤波重建
50
100
150
200
250
50100150200250
图18 模拟图像的高斯白噪声测试
original
无滤波重建
100
200
300
400
500
100200300400500
Ram-Lak 滤波重建
100
200
300
400
500
100200300400500
Shepp-Logan 滤波重建
100
200
300
400
500
100200300400500
图19 实际CT 图像的高斯白噪声测试
original
无滤波重建
50
100150200250
50100150200250
Ram-Lak 滤波重建
50
100
150
200
250
50100150200250
Shepp-Logan 滤波重建
50
100
150200250
50100150200250
图20 模拟图像的Poisson 噪声测试
original
无滤波重建
100
200
300
400
500
100200300400500
Ram-Lak 滤波重建
100
200
300
400
500
100200300400500
Shepp-Logan 滤波重建
100
200
300
400
500
100200300400500
图21 实际CT 图像的Poisson 噪声测试
六、结果分析
1、在没有噪声干扰或者噪声干扰很小的情况下，三种情况：无滤波、Ram-Lak 滤波重建、Shepp-Logan 滤波下的重建，无滤波重建得到的图像显示细节并没有另外两种优，后两种滤波重建得到的图像并没有很大的差异，如图14和图15所示。

2、不同的视角条件下，投影后重建的图像是有明显差异的，视角的间隔越小，最后重建得到的图像越优，如图16和图17所示。

3、在有噪声的干扰的情况下，首先测试高斯白噪声，结果表明，三种滤波方式均未能从噪声中很好地提取图像信息，如图18和图19所示，而在测试Poisson噪声情况下，无滤波重建显然不及后两种滤波投影优，而R-L和S-L均能达到理想的重建效果，但是相对来说，S-L的重建效果更有些，如图20和图21所示。

程序1
image= phantom('Modified Shepp-Logan',256); subplot(2,2,1);
imshow(image);
title('原始图像');
[R,xp]= radon(image,0:179);
theta = 0:180;
I1 = iradon(R,0:179,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray)
title('无滤波重建');
I2 = iradon(R,0:179,'linear','Ram-Lak'); subplot(2,2,3);
imagesc(I2);colormap(gray)
title('Ram-Lak滤波重建');
I3 = iradon(R,0:179,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray)
title('Shepp-Logan滤波重建');
程序2
image = imread('');
subplot(2,2,1);
imshow(image);
title('原始图像');
[R,xp]= radon(image,0:179);
theta = 0:180;
I1 = iradon(R,0:179,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray)
title('无滤波重建');
I2 = iradon(R,0:179,'linear','Ram-Lak'); subplot(2,2,3);
imagesc(I2);colormap(gray)
title('Ram-Lak滤波重建');
I3 = iradon(R,0:179,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray)
title('Shepp-Logan滤波重建');
figure(1)
a=0:10:180;
image = phantom('Modified Shepp-Logan',256); subplot(2,2,1);
imshow(image);
title('原始图像');
[R,xp] = radon(image,a);
theta = 0:18;
I1 = iradon(R,a,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray)
title('无滤波重建');
I2 = iradon(R,a,'linear','Ram-Lak');
subplot(2,2,3);
imagesc(I2);colormap(gray)
title('Ram-Lak滤波重建');
I3 = iradon(R,a,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray)
title('Shepp-Logan滤波重建');
figure(2)
b=0:2:180;
image = phantom('Modified Shepp-Logan',256); subplot(2,2,1);
imshow(image);
title('原始图像');
[R,xp]= radon(image,b);
theta = 0:36;
I1 = iradon(R,b,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray)
title('无滤波重建');
I2 = iradon(R,b,'linear','Ram-Lak');
subplot(2,2,3);
imagesc(I2);colormap(gray)
title('Ram-Lak滤波重建');
I3 = iradon(R,b,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray)
title('Shepp-Logan滤波重建');
figure(1);
image1= phantom('Modified Shepp-Logan',256); subplot(2,2,1);
imshow(image1);
title('original');
[R,xp]= radon(image1,0:179);
R=1e12*imnoise(1e-13*R,'gaussian',0,;
theta = 0:180;
I1 = iradon(R,0:179,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray)
title('无滤波重建');
I2 = iradon(R,0:179,'linear','Ram-Lak'); subplot(2,2,3);
imagesc(I2);colormap(gray)
title('Ram-Lak滤波重建');
I3 = iradon(R,0:179,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray)
title('Shepp-Logan滤波重建');
figure(2);
image2 = imread('');
subplot(2,2,1);
imshow(image2); title('original');
[R,xp]= radon(image2,0:179);
R=1e12*imnoise(1e-13*R,'gaussian',0,;
theta = 0:180;
I1 = iradon(R,0:179,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray);
title('无滤波重建');
I2 = iradon(R,0:179,'linear','Ram-Lak'); subplot(2,2,3);
imagesc(I2);colormap(gray);
title('Ram-Lak滤波重建');
I3 = iradon(R,0:179,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray);
title('Shepp-Logan滤波重建');
figure(1);
image1= phantom('Modified Shepp-Logan',256); subplot(2,2,1);
imshow(image1);
title('original');
[R,xp]= radon(image1,0:179);
R=1e12*imnoise(1e-13*R,'Poisson');
theta = 0:180;
I1 = iradon(R,0:179,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray)
title('无滤波重建');
I2 = iradon(R,0:179,'linear','Ram-Lak'); subplot(2,2,3);
imagesc(I2);colormap(gray)
title('Ram-Lak滤波重建');
I3 = iradon(R,0:179,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray)
title('Shepp-Logan滤波重建');
figure(2);
image2 = imread('');
subplot(2,2,1);
imshow(image2); title('original');
[R,xp]= radon(image2,0:179);
R=1e12*imnoise(1e-13*R,'Poisson');
theta = 0:180;
I1 = iradon(R,0:179,'linear','none');
subplot(2,2,2);
imagesc(I1);colormap(gray);
title('无滤波重建');
I2 = iradon(R,0:179,'linear','Ram-Lak'); subplot(2,2,3);
imagesc(I2);colormap(gray);
title('Ram-Lak滤波重建');
I3 = iradon(R,0:179,'linear','Shepp-Logan'); subplot(2,2,4);
imagesc(I3);colormap(gray);
title('Shepp-Logan滤波重建');。