《CT图像三维重建(附源码).doc》由会员分享，可在线阅读，更多相关《CT图像三维重建(附源码).doc（10页珍藏版）》请在三一文库上搜索。
1、程序流图:MATLAB 源码：clc;clear all;close all;% load mri% ph = squeeze(D);ph = pha ntom3d(128);prompt=E nter then ame=I nput nu mber：%载入mri数据，是matlab自带库 压缩载入的数据D,并赋值给phdefaulta nswer=1;numln put=in putdlg(prompt, name,1,defaulta nswer) %提示信息“输入1到27的片的数字” 弹出框名称默认数字弹出框，并得到用户的输入信息P= squeeze(ph(:,:,str2 num (c
2、ell2mat (numln put);% 从ph中得到相应的片信息P将用户的输入信息转换成数字，并imshow(P)D = 250;中心的像素距离。展示图片P将D赋值为250,是从扇束顶点到旋转dse nsorl = 2;%正实数指定扇束传感器的间距 2通过P, D等计算扇束的数据值正实数指定扇束传感器的间距1通过P, D等计算扇束的数据值F1 = fan beam(P,D,Fa nSen sorSpaci ng,dse nsorl); %dse nsor2 = 1;%F2 = fan beam(P,D,Fa nSen sorSpaci ng,dse nsor2); %dse nsorl =
3、 2;%正实数指定扇束传感器的间距 2dse nsorl = 2;%正实数指定扇束传感器的间距 2dse nsor3 = 0.25%正实数指定扇束传感器的间距0.25dse nsorl = 2;%正实数指定扇束传感器的间距 2F3, sen sor_pos3, fan _rot_a ngles3 = fan beam(P,D,.FanSensorSpacing,dsensor3);%通过P, D等计算扇束的数据值，并得到扇束传感器的位置 sensor_pos3和旋转角度fan_rot_angles3dse nsorl = 2;%正实数指定扇束传感器的间距 2dse nsorl = 2;%正实数
4、指定扇束传感器的间距 2figure,%imagesc(fa n_rot_a ngles3, sen sor_pos3, F3)%图片 colormap(hot);colorbar;创建窗口根据计算出的位置和角度展示F3的xlabel(Fa n Rotati on An gle (degrees) ylabel(Fa n Sen sor Positi on (degrees) output_size = max(size(P); output_sizeIfanl = ifan beam(F1,D, .Fa nSen sorSpac in g,dse nsor1,OutputSize,outpu
5、t_size);%设置色图为hot显示色栏定义x坐标轴定义y坐标轴得到P维数的最大值，并赋值给根据扇束投影数据 F1及D等数据重建图像figure, imshow(Ifa n1) title( 图一);disp(图一和原图的性噪比为：)；result=ps nr1(lfa n1,P);Ifan2 = ifan beam(F2,D, .Fa nSen sorSpac in g,dse nsor2,OutputSize,output_size);%创建窗口，并展示图片Ifanl根据扇束投影数据 F2及D等数据重建图像figure, imshow(Ifa n2)disp(图二和原图的性噪比为：)；r
6、esult=ps nr1(lfa n2,P);title( 图二);Ifan3 = ifan beam(F3,D, .Fa nSen sorSpac in g,dse nsor3,OutputSize,output_size);%创建窗口，并展示图片Ifan2根据扇束投影数据 F3及D等数据重建图像figure, imshow(Ifa n3) title( 图三);创建窗口，并展示图片Ifan3disp（图三和原图的性噪比为：）;result=ps nr1(lfa n3,P);function p,ellipse=pha ntom3d(vararg in) %PHANTOM3D Three-d
7、ime nsio nal an alogue of MATLAB Shepp-Logan pha ntom % P = PHANTOM3D(DEF,N) gen erates a 3D head pha ntom that ca n% be used to test 3-D recon structi on algorithms.% DEF is a stri ng that specifies the type of head pha ntom to gen erate.% Valid values are: %Shepp-Loga nA test image used widely by
8、researchers in%tomography% Modified Shepp-Loga n(default) A varia nt of the Shepp-Loga n pha ntom%in which the con trast is improved for better%visual percepti on.% % N is a scalar that specifies the grid size of P.% If you omit the argume nt, N defaults to 64.%P = PHANTOM3D(E,N) gen erates a user-d
9、efi ned pha ntom, where each row%of the matrix E specifies an ellipsoid in the image. E has ten colu mns,% with each colu mn containing a differe nt parameter for the ellipsoids:%Colu mn 1:Athe additive in te nsity value of the ellipsoid%Colu mn 2:athe len gth of the x semi-axis of the ellipsoid%Col
10、u mn 3:bthe len gth of the y semi-axis of the ellipsoid%Colu mn 4:cthe len gth of the z semi-axis of the ellipsoid%Colu mn 5:x0the x-coord in ate of the cen ter of the ellipsoid%Colu mn 6:y0the y-coord in ate of the cen ter of the ellipsoid%Colu mn 7:z0the z-coord in ate of the cen ter of the ellips
11、oid%Colu mn 8:phiphi Euler an gle (in degrees) (rotati on about z-axis)%Colu mn 9:thetatheta Euler an gle (in degrees) (rotati on about x-axis)%Colu mn 10:psipsi Euler an gle (in degrees) (rotati on about z-axis)%For purposes of gen erat ing the pha ntom, the doma ins for the x-, y-, and% z-axes spa
12、 n -1,1. Columns 2 through 7 must be specified in terms % of this ran ge.%P,E = PHANTOM3D(.) returns the matrix E used to gen erate the pha ntom.% Class Support% %All in puts must be of class double. All outputs are of class double.% Remarks% % For any give n voxel in the output image, the voxels va
13、lue is equal to the% sum of the additive inten sity values of all ellipsoids that the voxel is a% part of. If a voxel is not part of any ellipsoid, its value is 0.% The additive in te nsity value A for an ellipsoid can be positive or n egative;%if it is n egative, the ellipsoid will be darker tha n
14、the surro unding pixels.%Note that, depe nding on the values of A, some voxels may have values outside% the range 0,1.% Example% %ph = pha ntom3d(128);%figure, imshow(squeeze(ph(64,:,:)% Copyright 2005 Matthias Christian Schabel (matthias sta nfordalu mn i . org)%Uni versity of Utah Departme nt of R
15、adiology% Utah Cen ter for Adva need Imagi ng Research%729 Arapeen Drive%Salt Lake City, UT 84108-1218% This code is released un der the Gnu Public Lice nse (GPL). For more in formatio n,% see : http:/www.g nu.org/copyleft/gpl.html%Porti ons of this code are based on pha ntom.m, copyrighted by the M
16、athworks%ellipse ,n = parse_i nputs(vararg in :);p = zeros( n n n );rng =( (0: n-1)-( n-1)/2 ) / (n-1)/2);x,y,z = meshgrid(rng,rng,rng);coord = flatte n( x); flatte n( y); flatte n( z);p = flatte n( p);for k = 1:size(ellipse,1)A = ellipse(k,1);% Amplitude cha nge for this ellipsoidasq = ellipse(k,2)
17、A2;% aA2bsq=ellipse(k,3)A2;% bA2csq=ellipse(k,4)A2;% cA2x0 =ellipse(k,5);% x offsety0 =ellipse(k,6);% y offsetz0 =ellipse(k,7);% z offsetphi =ellipse(k,8)*pi/180;% first Euler an gle in radia nstheta = ellipse(k,9)*pi/180; % sec ond Euler an gle in radia nspsi =ellipse(k,10)*pi/180;% third Euler an
18、gle in radia nscphi = cos(phi); sphi = sin( phi); ctheta = cos(theta); stheta = si n( theta); cpsi = cos(psi); spsi = sin( psi);% Euler rotatio n matrixalpha = cpsi*cphi-ctheta*sphi*spsi cpsi*sphi+ctheta*cphi*spsi spsi*stheta; -spsi*cphi-ctheta*sphi*cpsi -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta; sthe
19、ta*sphi-stheta*cphictheta;% rotated ellipsoid coord in atescoordp = alpha*coord;idx = fin d(coordp(1,:)-x0)42./asq + (coordp(2,:)-y0)42./bsq + (coordp(3,:)-z0).A2./csq = 1); p(idx) = p(idx) + A;end p = reshape(p, n n n );return;function out = flatte n(i n)out = reshape(i n,1 prod(size(i n);return;fu
20、nction e,n = parse_ in puts(vararg in)% e is the m-by-10 array which defi nes ellipsoids% n is the size of the pha ntom brain imagen = 128;% The default sizee =；defaults = shepp-loga n, modified shepp-loga n, yu-ye-wa ng;for i=1: narginif ischar(vararg in i)% Look for a default pha ntomdef = lower(v
21、arargi ni);idx = strmatch(def, defaults);if isempty(idx)eid = spri ntf(lmages:%s: unknown Pha ntom,mfile name); msg = Unknown default pha ntom selected.; error(eid,%s,msg);endswitch defaultsidxcase shepp-loga ne = shepp_loga n;case modified shepp-loga ne = modified_shepp_loga n;case yu-ye-wa nge = y
22、u_ye_wa ng;endelseif nu mel(vararg in i)=1n = vararg in i;% a scalar is the image sizeelseif n dims(vararg in i)=2 & size(vararg in i,2)=10e = vararg in i;% user specified pha ntomelseeid = spri ntf(lmages:%s:i nvalidl nputArgs,mfile name);msg = In valid in put argume nts.; error(eid,%s,msg);endend%
23、 ellipse is not yet defi nedif isempty(e)e = modified_shepp_loga n;endreturn;% Default head pha ntoms:%function e = shepp_loga n e = modified_shepp_loga n;e(:,1) = 1 -.98 -.02 -.02 .01 .01 .01 .01 .01 .01;return;function e = modified_shepp_loga n%This head pha ntom is the same as the Shepp-Loga n ex
24、cept%the inten sities are cha nged to yield higher con trast in% the image. Taken from Toft, 199-200.% A a b c x0 y0 z0 phi theta psi% e = 1 .6900 .9201 .810000000-.8 .6624 .874.7800-.01840000-.2 .1100 .310.220.2200-18010-.2 .1600 .410.280-.220018010.1 .2100 .250.4100.35-.15000.1 .0460 .046.0500.1.2
25、5000.1 .0460 .046.0500-.1.25000.1 .0460 .023.050 -.08-.6050000.1 .0230 .023.0200-.6060000.1 .0230 .046.020.06-.6050000 ;return;function e = yu_ye_wa ng% Yu H, Ye Y, Wang G, Katsevich-Type Algorithms for Variable Radius SpiralCo ne-Beam CT% A a b c x0 y0 z0 phi theta psi% e = 1 .6900 .920 .900000000-
26、.8 .6624 .874 .880000000-.2 .4100 .160 .210 -.220-.2510800-.2 .3100 .110 .220.220-.257200.2 .2100 .250 .5000.35-.25000.2 .0460 .046 .0460.1-.25000.1 .0460 .023 .020 -.08-.65-.250001 .0460 .023 .020.06-.65-.2590002 .0560 .040 .100.06-.105.62590002 .0560 .056 .1000.100.625000 ;return;% func计算两幅图像的psnr
27、值function result=ps nr1(i n1,i n2) z=mse(i n1,i n2);result=10*log10(25592/z)% plot(result) fun ctio n z=mse(x,y) if n dims(x)=3 x=rgb2gray(x);endif n dims(y)=3 y=rgb2gray(y);endx=double(x); y=double(y);m1, n1=size(x);m2, n2=size(y); m=mi n( m1,m2); n=min(n1,n 2);z=0;for i=1:mfor j=1: n z=z+(x(i,j)-y(i,j)F2;endend z=z/(m* n);