CT图像重建.docx
- 文档编号:1117925
- 上传时间:2022-10-17
- 格式:DOCX
- 页数:15
- 大小:270.70KB
CT图像重建.docx
《CT图像重建.docx》由会员分享,可在线阅读,更多相关《CT图像重建.docx(15页珍藏版)》请在冰豆网上搜索。
CT图像重建
昆明理工大学信息工程与自动化学院学生实验报告
(2009—2010学年第一学期)
课程名称:
医学成像系统与放射治疗装置开课实验室:
32082008年12月24日
年级、专业、班
学号
姓名
成绩
实验项目名称
CT图像重建
指导教师
刘利军
教师评语
教师签名:
年月日
一、实验目的与意义
医学成像技术是生物医学工程专业的一门重要的专业课程,课程主要涉及X光仪器,CT仪器,MRI仪器和核医学仪器的工作原理及成像方法。
其中CT算法的出现又为后来数字化医学成像技术的发展提供了基础。
该门课程为生物医学工程专业的专业基础课。
CT技术是医学成像系统中的一种重要手段。
它通过特定的算法,利用计算机的高速运算功能,可以在短时间内快速呈现人体断层图像。
让学生练习CT图像的重建有助于学生理解CT算法的内容,熟悉数字图像重建的过程。
同时也能培养学生的团队精神和解决实际问题的能力。
二、实验算法原理
1、MATLAB处理数字图像的基本函数;
2、X-CT三维图像重建的基本算法。
CT图象重建有四种基本的算法:
矩阵法,迭代法,傅立叶算法,反投影算法.我们采用的方法为卷积反投影.卷积反投影有:
平行光束投影的卷积反投影算法,等角扇形光来投影的重建算法.
1).平行光束投影的卷积反投影算法
从投影重建三维物体的图像,就是重建一个个横断面。
这样三堆图像的重建就归结为二维图象的重建。
二维图像的重建问题可以从数学上描述如下。
假定表示一个二维的未知函数,通过的直线称为光钱(见图2.1)。
沿光线的积分称作光线积分。
沿相同方向的一组光线积分,就构成一个投影。
图2.1中垂直于直线(与X轴夹角为)的光线所形成。
图2.1在方向的投影
的投影,称之为在方向的投影。
光线积分和投影在数学上可以定义如下:
在图2.1中直线AB的方程为:
(2.1)
其中是AB到原点的距离,沿AB的积分为:
(2.2)
对于给定的,在方向的投影是t的函数。
如果在各个方向的投影已知,就可以唯一确定。
下面就讨论卷积反投影重建算法。
假定投影方向,如图2.2,将坐标旋转角(逆时针方向)形成坐标系。
在坐标系中为。
图2.2傅立叶切片定理示意图
坐标系与之间的关系为:
(2.3)
显然
(2.4)
令为的傅立叶变换则
(2.5)
将上式变换到坐标系中,注意到变换的可比行列式
(2.6)
从而得到:
(2.7)
其中
(2.8)
若令的傅立叶变换为,由(2.8)可知
(2.9)
若的傅立叶变换为的极坐标表示。
这说明在方向的投影
傅立叶变换等于在与u轴成角的直线上的值。
这就是著名的傅立叶投影切片定理。
可见在整个平面可以利用各个方向的投影来得到,从而也可以通过求的傅立叶反交换的办法求得:
(2.10)
变换到极坐标中
,
得到
(2.11)
经推导得
(2.12)
其中
若令
(2.13)
则
(2.14)
(2.13)式右端是两频谱函数和的乘积的傅立叶反变换。
是投影
傅立叶变换。
若的傅立叶反变换为,则根据卷积定理有:
(2.15)
或
其中
(2.16)
当图像的频谱是有限带宽时,则上式变为
(2.17)
由于图象及其频谱都是离散采样的,假定图象采样间隔为,则根据采样定理。
为了进行数学处理,只需知道h(t)在有限带宽上的离散采样点的值.这样我们有
(2.18)
其中n为正负整数。
(2.18)的离散形式为
(2.19)
假定在之外的值为0,则上式变为
(2.20)
或
(2.21)
其中从而可见为确定的N个采样点上的的值,需要使用的2N—1个点上的值,从n=一(N—1)到(N—1)。
为求得,利用傅立叶变换计算卷积是比较快的方法,为清除循环卷积的周期交叠效应,实际上取2N个点,补0,使之有(2N—1)个元素,则在N个采样点上就避免了交叠,如果使用以2为基的FFT(快速傅立叶变换)算法,和都必需朴0至(2N一1)个元素,(2N一1)为大于等于2N—l的最小的2的整数幂。
计算的过程可以写为
(2.22)
其中FFT和IFFT分别表示快速傅立叶变换和反变换,光滑窗是在滤波过程中加入的光滑因子,例如引用汉明窗,有时可以改进重建效果。
对于各个方向的投影,得到之后就可以由(2.22)来计算。
重建步骤可以归纳为:
第一步:
卷积,也称滤波,由(2.22)对每个方向计算。
第二步。
反投影,由(2.14)的近似形式
(2.23)
来计算的近似值。
M为投影个数为投影方向角,他们均匀的分布在0~的范围内。
当计算时,,不一定在的整离散点上,这就需要插值求得,预先将插值加密,即最靠近的点,可以提高计算速度。
2).等角扇形光来投影的重建算法
几乎所有的快遗CT设备都是用的扇形光束。
这里叙述的是等角度光束投影,如图2.3,测量投影数据的探测器等间距地分布在弧上,弧的半径为2D,D为光源到图像中心的距离。
在下文中,图象在极坐标中的表示。
表示在方向角为的投影中位量角为的光线产生的投影数据。
通过中心的光线其为0。
L表示从光源到像素的距离。
图2.3等扇形束投影重建算法中的变量
(2.24)
表示在方向角为的投影中通过像素的光线的位置角
(2.25)
图像和扇形投影有下述关系
(2.26)
其中和是投影中光线的最大位置角,从而可以得到这种重建算法的执行步骤:
第一步:
投影的修改,假定投影的抽样间隔为,抽样数据通过下式进行修改,
(2.27)
第二步:
卷积(滤波)将第一步修改了的投影与响应函数进行卷积
(2.28)
(2.29)
其离散形式为:
(2.30)
这里也和上节一样可以加进一光滑窗,改进重建效应。
第三步:
反投影,就是执行(2.30)的积分
(2.31)
近似的有:
(2.32)
其是图象的近似图像,,,M是投影数,假定投影均匀地分布在0~内。
为求得滤波后的投影,对象素的贡献,首先由(2.24)和(2.25)计算出和.所有对的贡献加起来再乘以就得。
四、实验程序代码
原图象代码:
p=phantom(256);
imshow(p)
卷积反投影法代码:
H=[000.920.6990*pi/1801;...
0-0.01840.8740.662490*pi/180-0.98;...
0.2200.310.1172*pi/180-0.2;...
-0.2200.410.16108*pi/180-0.2;...
00.350.250.2190*pi/1800.1;...
00.10.0460.04600.2;...
0-0.10.0460.04600.2;...
-0.08-0.6050.0460.02300.1;...
0-0.6050.0230.02300.1;...
0.06-0.6050.0460.02390*pi/1800.1];
angle=1;
N=100;
vstep=2/N;
ax=N/2+1;
fork=1:
(180/angle+1)
theta=(k-1)*angle*pi/180;
fori=1:
10
x0=H(i,1);y0=H(i,2);A=H(i,3);B=H(i,4);alpha=H(i,5);rho=H(i,6);
R=0;
forw=ax;
back=ax;
MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(theta))^2);
NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2;
g(i,ax)=rho*2*A*B*MM/NN;
forj=1:
N/2
R=R+vstep;
forw=forw+1;
MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(theta))^2);
NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2;
g(i,forw)=rho*2*A*B*MM/NN;
R=-R;
back=back-1;
MM=sqrt(A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2-(R-x0*cos(theta)-y0*sin(theta))^2);
NN=A^2*cos(theta-alpha)^2+B^2*sin(theta-alpha)^2;
g(i,back)=rho*2*A*B*MM/NN;
R=-R;
end
end
radon0(k,:
)=real(sum(g));
end
%%%%%%%%%%%%generateR-Lfunction%%%%%%%%%%%%%%%%
radon1=[zeros(k,N),radon0];
ax=N+1;
RL(ax)=1/(4*vstep^2);
forw=ax+1;
back=ax-1;
fork=1:
N/2
n=2*k-1;
RL(forw)=-1/(n*pi*vstep)^2;
RL(back)=RL(forw);
RL(forw+1)=0;
RL(back-1)=0;
forw=forw+2;
back=back-2;
end
fork=1:
(180/angle+1)
radon2(k,:
)=conv(radon1(k,:
),RL);
end
radonf=radon2(:
2*N:
3*N);
%%%%%%%%%%%%generateS-Lfunction%%%%%%%%%%%%%%%%%%
radon1=[zeros(k,N),radon0];
forv=1:
(2*N+1)
n=v-N-1;
SL(v)=-2/(pi^2*vstep^2*(4*n^2-1));
end
fork=1:
(180/angle+1)
radon2(k,:
)=conv(radon1(k,:
),SL);
end
radonf=radon2(:
2*N:
3*N);
figure
(1)
subplot(321)
plot(1:
(2*N+1),radon1(1,:
))
title('投影函数(已补零)')
subplot(323)
plot(1:
(2*N+1),SL)
title('S-L卷积函数')
subplot(325)
plot(1:
(4*N+1),radon2(1,:
))
title('卷积结果')
Xradon1=fft(radon1(1,:
));
subplot(322)
plot(1:
(2*N+1),abs(Xradon1))
title('频谱')
XRL=fft(SL);
subplot(324)
plot(1:
(2*N+1),abs(XRL))
Xradon2=fft(radon2(1,:
));
subplot(3
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- CT 图像 重建