S-L头模型程序

x=linspace(-1,1,128);y=linspace(-1,1,128);[x1,y1]=meshgrid(x,y);x0=[000.22-0.22000-0.0800.06];y0=[0-0.0184000.350.1-0.1-0.605-0.605-0.605];a=[0.920.8740.310.410.250.0460.0460.0460.0230.046];b=[0.690.66240.110.160.210.0460.0460.0230.0230.023];alpha=[90907210890000090]/180*pi;p=[2.0-0.98-0.2-0.20.10.10.10.10.10.1];%%%%%%%%%模型数据产生F=0;fori=1:10X=(x1-x0(i)).*cos(alpha(i))+(y1-y0(i)).*sin(alpha(i));Y=-(x1-x0(i)).*sin(alpha(i))+(y1-y0(i)).*cos(alpha(i));R0=(X/a(i)).^2+(Y/b(i)).^2;f=[R0=1].*p(i);F=F+f;end%处理显示数据%forj=1:64%矩阵上下翻转%t=F(j,1:128);%F1(j,1:128)=F(129-j,1:128);%F1(129-j,1:128)=t;%%endF1=flipud(F);%矩阵上下翻转F2=mat2gray(F1)%矩阵转换为灰度图像figure(1);imshow(F2);%显示仿真头模型title('S-L头模型的灰度显示图像');%投影数据产生N=128;%%%%%%%%%%%%%vstep=2.0/N;ax=N/2+1;fortheta=0:180forl=1:10R=0;forw=ax;back=ax;MM=sqrt(a(l)^2.*cos(theta*pi/180-alpha(l))^2+b(l)^2.*sin(theta*pi/180-alpha(l))^2-(R-x0(l).*cos(theta*pi/180)-y0(l).*sin(theta*pi/180))^2);NN=a(l)^2.*cos(theta*pi/180-alpha(l))^2+b(l)^2.*sin(theta*pi/180-alpha(l))^2;g(l,ax)=2*p(l)*a(l)*b(l)*MM/NN;form=1:N/2R=R+vstep;forw=forw+1;MM=sqrt(a(l)^2.*cos(theta*pi/180-alpha(l))^2+b(l)^2.*sin(theta*pi/180-alpha(l))^2-(R-x0(l).*cos(theta*pi/180)-y0(l).*sin(theta*pi/180))^2);NN=a(l)^2.*cos(theta*pi/180-alpha(l))^2+b(l)^2.*sin(theta*pi/180-alpha(l))^2;g(l,forw)=2*p(l).*a(l).*b(l).*MM./NN;R=-R;back=back-1;MM=sqrt(a(l)^2.*cos(theta*pi/180-alpha(l))^2+b(l)^2.*sin(theta*pi/180-alpha(l))^2-(R-x0(l).*cos(theta*pi/180)-y0(l).*sin(theta*pi/180))^2);NN=a(l)^2.*cos(theta*pi/180-alpha(l))^2+b(l)^2.*sin(theta*pi/180-alpha(l))^2;g(l,back)=2*p(l).*a(l).*b(l).*MM./NN;R=-R;endendradon0(1:129,theta+1)=real(sum(g));%%%%%%%%%%%%%%%%endfigure(2);imshow(radon0,[]);%显示投影数据title('显示投影数据');%卷积反投影M=128;daltax=2.0/M;daltay=2.0/M;dalta=2.0/M;T=(F)';%转置%fori=1:M%forj=1:M%T(i,j)=t(-1.0+i*daltax-daltax/2,-1.0+j*daltay-daltay/2);%end%end%求h（n）ax=M+1;h(ax)=1/(4*dalta^2);forw=ax+1;back=ax-1;forn=1:2:127h(forw)=-1/(n*pi*dalta)^2;h(back)=h(forw);h(forw+1)=0;h(back-1)=0;forw=forw+2;back=back-2;end%求g（n）%radon1=[zeros(181,M),radon0']';%上补零%figure(3);%imshow(radon1,[]);radon1=[radon0',zeros(181,M)]';%或下补零%求卷积fortheta=0:180Q1(:,theta+1)=conv(radon1(:,theta+1),h);end%Q2=Q1(2*M:3*M,:);%上补零%figure(5);%imshow(Q2,[]);Q2=Q1(M:2*M,:);%或下补零figure(3);imshow(Q2,[]);%显示最终卷积运算结果title('显示最终卷积运算结果');%反投影%fori=1:128%forj=1:128%T(i,j)=0.0;%fortheta=0:180%C=M/2-(M-1)*(cos(theta*pi/180)+sin(theta*pi/180))/2;%R0=(i-1)*cos(theta*pi/180)+(j-1)*sin(theta*pi/180)+C%n0=floor(R0);%ifn00&&n0(M+1)%d=R-n0;%T(i,j)=T(i,j)+(1-d)*Q2(n0,theta+1)+d*Q2(n0+1,theta+1);%else%T(i,j)=T(i,j);%end%end%end%endfortheta=0:180C=M/2-(M-1)*(cos(theta*pi/180)+sin(theta*pi/180))/2;fori=1:MR=(i-1)*cos(theta*pi/180)+C;n0=floor(R);ifn00&&n0(N+1)dot=R-n0;T(i,1,theta+1)=(1-dot)*Q2(n0,theta+1)+dot*Q2(n0+1,theta+1);elseT(i,1,theta+1)=0;endforj=2:MR=R+sin(theta*pi/180);n0=floor(R);ifn00&&n0(N+1)dot=R-n0;T(i,j,theta+1)=(1-dot)*Q2(n0,theta+1)+dot*Q2(n0+1,theta+1);elseT(i,j,theta+1)=0;endendendendT2=sum(T,3);Gfinal=mat2gray(T2);Gfinal=imrotate(Gfinal,90);figure(4);imshow(Gfinal);title('显示卷积反投影重建的结果');