CT实验一Word文件下载.docx
《CT实验一Word文件下载.docx》由会员分享，可在线阅读，更多相关《CT实验一Word文件下载.docx（26页珍藏版）》请在冰点文库上搜索。
f3=fftshift（fab）;
%对所得的傅立叶变换进行频谱转移
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:
imshow（C,[]）;
figure,imshow（I,[]）;
得出结果如图2-1所示。
图2-1a弦图biradon重建图像
2）首先，我们分别对前面所创建的shepp-logan进行平行束反投影重建以及滤波反投影重建。
（1）反投影重建
Step1.生成shepp-logan体模图像，并得到它的投影值。
I=phantom（256）;
%formtheimage
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:
[X,Y]=meshgrid（x,y）;
xpr=X-（sideSize+1）/2;
ypr=Y-（sideSize+1）/2;
%loopovereachprojection
fori=1:
m
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;
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）;
ifMY>
hy=（roi（4）-roi（3））/（MY-1）;
yrange=flipud（（roi（3）+hy*[0:
MY-1]）'
yrange=roi（3）;
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）;
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:
来将文本格式matlab并以图像形式输出。
并用D=radon（I,0:
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,...
'
0.5,...
FanCoverage'
'
minimal'
...
FanRotationIncrement'
1）;
phReconstructed=ifanbeam（F,200,...
Filter'
Shepp-Logan'
OutputSize'
256,...
imshow（phReconstructed）;
图2-7a认值重建图像b改变参数后重建图像
由上图可以得到，通过修改参数，可以使图像重建结果更加清晰。
升级会员 