CT实验一.docx

1、CT实验一实验课程：现代医学影像设备 生物医学工程学院实验技术管理中心实验题目：1. 傅立叶切片定理证明实验目的：加深对傅立叶切片定理的理解以及强化matlab编程能力。实验原理：傅立叶切片定理实验内容：用matlab证明在x线束平行于x轴时傅立叶切片定理的正确性实验步骤1.运行matlab程序，在命令窗口输入以下exit moxing输入以下代码：I=phantom(512); /生成shepp-logan体模H=imrotate(I,270); /figure,imshow(I);figure,imshow(H);显示图像如下：图1 shepp-ladon体模 512X512图2体模旋转1

2、80度 矩阵512X5122、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); %显示图像图33)直接把原旋转180角度下的体模进行一维傅立叶变换，图如下：f4=

3、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平行线束反投影重建及滤波

4、反投影重建实验目的：更深入的掌握反投影以及滤波反投影重建方法及其matlab实现实验原理：滤波反投影的重建算法:原理：“断层平面中某一点的密度值可看作这一平面内所有经过该点的射线投影之和（的平均值）” 先对1D投影进行滤波,然后再进行BP实验内容：1. 实验步骤1) 通过 help radon 和 help iradon 了解radon与iradon函数的用法; 然后编写生成图像的m文件；程序：P = phantom(256);C = radon(P,0:179);I = iradon(C, 0:179);imshow(C,);figure,imshow(I,);得出结果如图2-1所示。 图2

5、-1 a 弦图 b iradon重建图像2) 首先，我们分别对前面所创建的shepp-logan进行平行束反投影重建以及滤波反投影重建。（1） 反投影重建 Step1.生成shepp-logan体模图像，并得到它的投影值。I=phantom(256); %form the imagefigure,imshow(I);IMG=double(I);THETA=0:179;% This MATLAB code takes an image matrix and vector of angles and then finds the 1D projection (Radon transform) at

6、 each of the angles. It returns a matrix whose columns are the projections at each angle.pad the image with zeros so we dont lose anything when we rotate.iLength, iWidth = size(IMG);iDiag = sqrt(iLength2 + iWidth2);LengthPad = ceil(iDiag - iLength) + 2;WidthPad = ceil(iDiag - iWidth) + 2;padIMG = ze

7、ros(iLength+LengthPad, iWidth+WidthPad);padIMG(ceil(LengthPad/2):(ceil(LengthPad/2)+iLength-1), . ceil(WidthPad/2):(ceil(WidthPad/2)+iWidth-1) = IMG;% loop over the number of angles, rotate 90-theta (because we can easily sum if we look at stuff from the top), and then add up. Dont perform any inter

8、polation on the rotating.n = length(THETA);PR = zeros(size(padIMG,2), n);for i = 1:n tic tmpimg = imrotate(padIMG, 90-THETA(i), bilinear, crop); PR(:,i) = (sum(tmpimg); THETA(i) tocendfigure,imshow(PR,);得到的图像如图2-1所示 图2-2 a原始体模 b所得的投影图Step2.对所得的图2-2b进行反投影重建。程序： % back projectionn = size(PR,1);sideSiz

9、e = n; % convert THETA to radiansth = (pi/180)*THETA;% set up the imagem = length(THETA); BPI = zeros(sideSize,sideSize);% find the middle index of the projectionsmidindex = (n+1)/2;% set up x and y matricesx = 1:sideSize;y = 1:sideSize;X,Y = meshgrid(x,y);xpr = X - (sideSize+1)/2;ypr = Y - (sideSiz

10、e+1)/2; % loop over each projection for i = 1:m tic disp(On angle , num2str(THETA(i); filtIndex = round(midindex + xpr*sin(th(i) - ypr*cos(th(i); BPIa = zeros(sideSize,sideSize); spota = find(filtIndex 0) & (filtIndex 0;RF(ind) = RF(ind)+rho(mu)*(2*u(mu)*v(mu)*sqrt(test(ind)/a;end % mu-loopslkernel.

11、m：function u = slkernel(t)u = zeros(size(t);i1 = 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*pi3);window3.m:function pic1 = window3(mi,ma,roi,pic);%function pic1 = window3(mi,ma,roi,pic);% displays image pic with coordinates given by roi% roi

12、 = xmin xmax ymin ymaxx = 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 1 hx = (roi(2)-roi(1)/(MX-1); xrange = roi(1) + hx*0:MX-1;else hx = 0; xrange = roi(1);endif MY 1 hy = (roi(4)-roi(3

13、)/(MY-1); yrange = flipud(roi(3) + hy*0:MY-1);else hx = 0; yrange = roi(3);endcenter = (roi(1)+roi(2), (roi(3)+roi(4)/2;x1 = ones(MY,1)*xrange; %x-coordinate matrixx2 = yrange*ones(1,MX); %y-coordinate matrixif circle = 1 re = min(roi(2)-roi(1),roi(4)-roi(3)/2; chi = (x1-center(1).2 + (x2-center(2).

14、2 = re2); %char. fct of roi;else chi = isfinite(x1);endx1 = 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)/(rps2); %compute discrete convolution kernel.for j = 1:p j phi = (pi*(j-1)/p); theta = cos(phi);sin(phi); RF = Rad(theta,ph

15、i,s,xe,ye,ae,be,alpha,rho); %compute line integrals% convolution C = conv(RF,wb); Q = h*C(2*q+1:4*q); Q(2*q+1)=0; % interpolation and backprojection 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=

16、Pchi+ Pupdate;end % j-loopP(chi) = Pchi*(2*pi/p);pmin = min(min(P); pmax = max(max(P);window3(pmin,pmax,roi,P); % view the computed image所得的实验结果如图2-4所示 图2-4滤波反投影的结果实验结果：由上图可知，反投影重建和滤波反投影重建都可以重建图像，但是反投影的结果会使图像模糊，滤波反投影的结果则有较好的效果。实验讨论：1) 根据实验指导图2-3和图像2-4，试说明为何反投影重建会出现图2-3的模糊效果，并将反投影重建和滤波反投影重建加以比较。讨论：反投

17、影法是利用穿过某些像素的所有射线的投影值反过来估算该像素的吸收系数值。这样就把从各个方向上得到的投影看作为这个方向上的像素具有同等的贡献，邻近结构对反投影图像中所有点的密度都要产生影响。以下面的图为例：图中的黑点表示一个射线密集的物质，此时X线管和探测器面对面地沿着被检体扫描，这些信号射线的反投影，产生一个星状图案，焦点处相当于铁钉的高密度处，其周围则有模糊阴影存在，此种阴影的浓淡程度随着与高密度物体的距离远近成反比例减少，从而造成重建后图像的模糊。2）. 所给的文件夹projdata.zip中会给出几组实际临床图像，以.txt形式给出，以临床腹部图像为例，首先对文件进行更改：把图像前的数字去

18、除，利用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 的弦图滤波反投影重建图像结果：第三

19、幅图像：dsitest: distest 的原图像 distest 的弦图滤波反投影重建图像结果：第四幅图像：lung: lung 的原图像 lung 的弦图滤波反投影重建图像结果：实验三：扇束重建实验实验题目：3. 扇束重建实验(选作)实验目的：验证扇束重建的原理及其matalb函数。实验原理：扇形线束滤波反投影重建算法实验内容：理解matlab中fanbeam和ifanbeam函数的参数，通过修改参数，对比图像重建的结果并分析原因。实验步骤：扇形束滤波反投影重建Step1.打开matlab 在command window 中键入 help fanbeam 与 help ifanbeam ，

20、了解扇形线束重建函数以及相关参数。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-5 a 平行束重建结果 b 扇形线束重建结果图2-6扇形重建投影图Step3.自定义matlab中ifanbeam函数的参数，我们采用等角度扇束几何结构，即保持i

21、fanbeam的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-7 a认值重建图像 b改变参数后重建图像实验结果：由上图可以得到，通过修改参数，可以使图像重建结果更加清晰。实验讨论：