CT实验一.docx
- 文档编号:5389705
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:26
- 大小:487.78KB
CT实验一.docx
《CT实验一.docx》由会员分享,可在线阅读,更多相关《CT实验一.docx(26页珍藏版)》请在冰点文库上搜索。
CT实验一
实验课程:
现代医学影像设备
生物医学工程学院实验技术管理中心
实验题目:
1.傅立叶切片定理证明
实验目的:
加深对傅立叶切片定理的理解以及强化matlab编程能力。
实验原理:
傅立叶切片定理
实验内容:
用matlab证明在x线束平行于x轴时傅立叶切片定理的正确性
实验步骤
1.运行matlab程序,在命令窗口输入以下exitmoxing
输入以下代码:
I=phantom(512);//生成shepp-logan体模
H=imrotate(I,270);//
figure,imshow(I);
figure,imshow(H);
显示图像如下:
图1shepp-ladon体模512X512
图2体模旋转180度矩阵512X512
2、X线以与Y轴平行的方向进行投影,根据傅立叶切片定理知道
=180度,
1)获取在
=180度角度下的体模的投影数据
projection=sum(H,1);
figure,plot(1:
512,projection);
图如下:
2)再把此角度下的投影数据进行傅立叶变换,图如下
f=fft(projection);%进行傅立叶变换
fab=f.*conj(f);%由于变换是复数,取模
f3=fftshift(fab);%对所得的傅立叶变换进行频谱转移
figure,plot(1:
512,f3);%显示图像
图3
3)直接把原旋转180角度下的体模进行一维傅立叶变换,图如下:
f4=fft2(H);%直接对原图傅立叶变换
fab1=f4.*conj(f4);%取模
f5=fftshift(fab1);
figure,plot(f5);
图4
实验结论
由图3和图4可知,在θ=180度时,傅立叶切片定理是成立的。
实验讨论
把程序中H=imrotate(I,270);或者其他角度并且把H和I互换,可知在θ=180度时,傅立叶切片定理是成立的。
由于坐标系是任意选择的,当把坐标系旋转一个角度时,以上结论仍然是正确的。
也即是:
目标物质在任意角度的投影的傅里叶变换,等于它的二维傅里叶变换在同方向的直线。
评阅成绩
2.CT平行线束反投影重建及滤波反投影重
实验题目:
2.CT平行线束反投影重建及滤波反投影重建
实验目的:
更深入的掌握反投影以及滤波反投影重建方法及其matlab实现
实验原理:
滤波反投影的重建算法:
原理:
“断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和(的平均值)”先对1D投影进行滤波,然后再进行BP
实验内容:
1.实验步骤
1)通过helpradon和helpiradon了解radon与iradon函数的用法;
然后编写生成图像的m文件;程序:
P=phantom(256);
C=radon(P,0:
179);
I=iradon(C,0:
179);
imshow(C,[]);figure,imshow(I,[]);
得出结果如图2-1所示。
图2-1a弦图biradon重建图像
2)首先,我们分别对前面所创建的shepp-logan进行平行束反投影重建以及滤波反投影重建。
(1)反投影重建
Step1.生成shepp-logan体模图像,并得到它的投影值。
I=phantom(256);%formtheimage
figure,imshow(I);
IMG=double(I);
THETA=0:
179;
%ThisMATLABcodetakesanimagematrixandvectorofanglesandthenfindsthe1Dprojection(Radontransform)ateachoftheangles.Itreturnsamatrixwhosecolumnsaretheprojectionsateachangle.padtheimagewithzerossowedon'tloseanythingwhenwerotate.
[iLength,iWidth]=size(IMG);
iDiag=sqrt(iLength^2+iWidth^2);
LengthPad=ceil(iDiag-iLength)+2;
WidthPad=ceil(iDiag-iWidth)+2;
padIMG=zeros(iLength+LengthPad,iWidth+WidthPad);
padIMG(ceil(LengthPad/2):
(ceil(LengthPad/2)+iLength-1),...
ceil(WidthPad/2):
(ceil(WidthPad/2)+iWidth-1))=IMG;
%loopoverthenumberofangles,rotate90-theta(becausewecaneasilysumifwelookatstufffromthetop),andthenaddup.Don'tperformanyinterpolationontherotating.
n=length(THETA);
PR=zeros(size(padIMG,2),n);
fori=1:
n
tic
tmpimg=imrotate(padIMG,90-THETA(i),'bilinear','crop');
PR(:
i)=(sum(tmpimg))';
THETA(i)
toc
end
figure,imshow(PR,[]);
得到的图像如图2-1所示
图2-2a原始体模b所得的投影图
Step2.对所得的图2-2b进行反投影重建。
程序:
%backprojection
n=size(PR,1);
sideSize=n;
%convertTHETAtoradians
th=(pi/180)*THETA;
%setuptheimage
m=length(THETA);
BPI=zeros(sideSize,sideSize);
%findthemiddleindexoftheprojections
midindex=(n+1)/2;
%setupxandymatrices
x=1:
sideSize;
y=1:
sideSize;
[X,Y]=meshgrid(x,y);
xpr=X-(sideSize+1)/2;
ypr=Y-(sideSize+1)/2;
%loopovereachprojection
fori=1:
m
tic
disp(['Onangle',num2str(THETA(i))]);
filtIndex=round(midindex+xpr*sin(th(i))-ypr*cos(th(i)));
BPIa=zeros(sideSize,sideSize);
spota=find((filtIndex>0)&(filtIndex<=n));
newfiltIndex=filtIndex(spota);
BPIa(spota)=PR(newfiltIndex(:
),i);
BPI=BPI+BPIa;
toc
end
BPI=BPI./m;
figure,imshow(BPI,[]);
得到的结果如图2-3
图2-3反投影重建结果
(2)滤波反投影重建
滤波反投影就是在反投影之前对投影数据进行前置滤波,实验的代码如下:
先生成三个m文件:
Rad.m:
function[RF]=Rad(theta,phi,s,x,y,u,v,alpha,rho)
%ThisfunctioncomputestheRadontransformofellipses
%centeredat(x,y)withmajoraxisu,minoraxisv,
%rotatedthroughanglealpha,withweightrho.
RF=zeros(size(s));
formu=1:
max(size(x))
a=(u(mu)*cos(phi-alpha(mu)))^2+(v(mu)*sin(phi-alpha(mu)))^2;
test=a-(s-[x(mu);y(mu)]'*theta).^2;
ind=test>0;
RF(ind)=RF(ind)+rho(mu)*(2*u(mu)*v(mu)*sqrt(test(ind)))/a;
end%mu-loop
slkernel.m:
functionu=slkernel(t)
u=zeros(size(t));
i1=abs(abs(t)-pi/2)<=1.e-6;
u(i1)=ones(size(u(i1)))/pi;
t1=t(abs(abs(t)-pi/2)>1.e-6);
v=(pi/2-t1.*sin(t1))./((pi/2)^2-t1.^2);
u(abs(abs(t)-pi/2)>1.e-6)=v;
u=u/(2*pi^3);
window3.m:
functionpic1=window3(mi,ma,roi,pic);
%functionpic1=window3(mi,ma,roi,pic);
%displaysimagepicwithcoordinatesgivenbyroi
%roi=[xminxmaxyminymax]
x=[roi
(1),roi
(2)];y=[roi(3),roi(4)];
colors=128;co=colors-1;
pic1=pic-mi*ones(size(pic));
pic1=(co/(ma-mi))*pic1;
P=(pic1>=0);
pic1=pic1.*P;
P=(pic1<=co);
pic1=pic1.*P+co*(ones(size(pic1))-P);
colormap(gray(colors));
image(x,fliplr(y),flipud(pic1));
axis('square');
进行滤波反投影,源程序如下:
%Parallel-beamfilteredbackprojectionalgorithm
%forthestandardlattice
%Lastrevision:
August29,2001
%specifyinputparametershere
p=200;%numberofviewanglesbetween0andpi
q=64;%q=1/d,d=detectorspacing
MX=128;MY=128;%matrixdimensions
roi=[-11-11];%roi=[xminxmaxyminymax]
%regionofinterestwhere
%reconstructioniscomputed
circle=1;%Ifcircle=1imagecomputedonlyinside
%circleinscribedinroi.
%Specifyparametersofellipsesformathematicalphantom.
%xe=vectorofx-coordinatesofcenters
%ye=vectorofy-coordinatesofcenters
%ae=vectoroffirsthalfaxes
%be=vectorofsecondhalfaxes
%alpha=vectorofrotationangles(degrees)
%rho=vectorofdensities
xe=[000.22-0.22000-0.0800.060.5538];
ye=[0-0.0184000.350.1-0.1-0.605-0.605-0.605-0.3858];
ae=[0.690.66240.110.160.210.0460.0460.0460.0230.0230.0333];
be=[0.920.8740.310.410.250.0460.0460.0230.0230.0460.206];
alpha=[00-1818000000-18];
rho=[1-0.98-0.02-0.020.010.010.010.010.010.010.03];
%endofinputsection
b=pi*q;rps=1/b;
alpha=alpha*pi/180;
ifMX>1
hx=(roi
(2)-roi
(1))/(MX-1);
xrange=roi
(1)+hx*[0:
MX-1];
else
hx=0;xrange=roi
(1);
end
ifMY>1
hy=(roi(4)-roi(3))/(MY-1);
yrange=flipud((roi(3)+hy*[0:
MY-1])');
else
hx=0;yrange=roi(3);
end
center=[(roi
(1)+roi
(2)),(roi(3)+roi(4))]/2;
x1=ones(MY,1)*xrange;%x-coordinatematrix
x2=yrange*ones(1,MX);%y-coordinatematrix
ifcircle==1
re=min([roi
(2)-roi
(1),roi(4)-roi(3)])/2;
chi=((x1-center
(1)).^2+(x2-center
(2)).^2<=re^2);%char.fctofroi;
else
chi=isfinite(x1);
end
x1=x1(chi);x2=x2(chi);
P=zeros(MY,MX);Pchi=P(chi);
h=1/q;
s=h*[-q:
q-1];
bs=[-2*q:
2*q-1]/(q*rps);
wb=slkernel(bs)/(rps^2);%computediscreteconvolutionkernel.
forj=1:
p
j
phi=(pi*(j-1)/p);
theta=[cos(phi);sin(phi)];
RF=Rad(theta,phi,s,xe,ye,ae,be,alpha,rho);%computelineintegrals
%convolution
C=conv(RF,wb);
Q=h*C(2*q+1:
4*q);Q(2*q+1)=0;
%interpolationandbackprojection
Q=[real(Q)';0];
t=theta
(1)*x1+theta
(2)*x2;
k1=floor(t/h);
u=(t/h-k1);
k=max(1,k1+q+1);k=min(k,2*q);
Pupdate=((1-u).*Q(k)+u.*Q(k+1));
Pchi=Pchi+Pupdate;
end%j-loop
P(chi)=Pchi*(2*pi/p);
pmin=min(min(P));
pmax=max(max(P));
window3(pmin,pmax,roi,P);%viewthecomputedimage
所得的实验结果如图2-4所示
图2-4滤波反投影的结果
实验结果:
由上图可知,反投影重建和滤波反投影重建都可以重建图像,但是反投影的结果会使图像模糊,滤波反投影的结果则有较好的效果。
实验讨论:
1)根据实验指导图2-3和图像2-4,试说明为何反投影重建会出现图2-3的模糊效果,并将反投影重建和滤波反投影重建加以比较。
讨论:
反投影法是利用穿过某些像素的所有射线的投影值反过来估算该像素的吸收系数值。
这样就把从各个方向上得到的投影看作为这个方向上的像素具有同等的贡献,邻近结构对反投影图像中所有点的密度都要产生影响。
以下面的图为例:
图中的黑点表示一个射线密集的物质,此时X线管和探测器面对面地沿着被检体扫描,这些信号射线的反投影,产生一个星状图案,焦点处相当于铁钉的高密度处,其周围则有模糊阴影存在,此种阴影的浓淡程度随着与高密度物体的距离远近成反比例减少,从而造成重建后图像的模糊。
2).所给的文件夹projdata.zip中会给出几组实际临床图像,以.txt形式给出,以临床腹部图像为例,
首先对文件进行更改:
把图像前的数字去除,利用matlab语句:
%事先将文本放入matlab当前目录
I=load('abdom.txt');%载入文本
I=double(I);%数据类型转换
I=iradon(I,0:
179);
figure,imshow(I,[]);
来将文本格式matlab并以图像形式输出。
并用D=radon(I,0:
179);figure,imshow(D);得到图像的弦图。
理解所给平行束的滤波反投影算法,并根据算法重建四副临床图像。
第一幅图像:
abdom:
abodm的原图像
abodm的弦图
滤波反投影重建图像结果:
第二幅图像:
brain:
brain的原图像
brain的弦图
滤波反投影重建图像结果:
第三幅图像:
dsitest:
distest的原图像distest的弦图
滤波反投影重建图像结果:
第四幅图像:
lung:
lung的原图像
lung的弦图
滤波反投影重建图像结果:
实验三:
扇束重建实验
实验题目:
3.扇束重建实验(选作)
实验目的:
验证扇束重建的原理及其matalb函数。
实验原理:
扇形线束滤波反投影重建算法
实验内容:
理解matlab中fanbeam和ifanbeam函数的参数,通过修改参数,对比图像重建的结果并分析原因。
实验步骤:
扇形束滤波反投影重建
Step1.打开matlab在commandwindow中键入helpfanbeam与helpifanbeam,了解扇形线束重建函数以及相关参数。
Step2.在新的m-file中输入对shepp-logan进行扇束滤波反投影重建。
程序:
ph=phantom(256);
d=200;
F=fanbeam(ph,d);
I=ifanbeam(F,d,'FanSensorSpacing',0.25);
figure,imshow(F,[]);
figure,imshow(I);
figure,imshow(ph);
实验结果如图:
图2-5a平行束重建结果b扇形线束重建结果
图2-6扇形重建投影图
Step3.自定义matlab中ifanbeam函数的参数,我们采用等角度扇束几何结构,即保持ifanbeam的FanSensorGeometry默认值arc。
在插值方式和滤波器不变的前提下,我们通过自行扇束探测器间距离值,探测器重建角度以及探测器角度增量来对比重建图像如图2-7.
修改参数的代码如下:
ph=phantom(256);
P=radon(ph);
[F,obeta,otheta]=para2fan(P,200,...
'FanSensorSpacing',0.5,...
'FanCoverage','minimal',...
'FanRotationIncrement',1);
phReconstructed=ifanbeam(F,200,...
'FanSensorSpacing',0.5,...
'Filter','Shepp-Logan',...
'OutputSize',256,...
'FanCoverage','minimal',...
'FanRotationIncrement',1);
imshow(phReconstructed);
图2-7a认值重建图像b改变参数后重建图像
实验结果:
由上图可以得到,通过修改参数,可以使图像重建结果更加清晰。
实验讨论:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- CT 实验