医学成像系统实验报告姓名:__陆想想_______学号:___030840110____指导老师:___陶玲____南京航空航天大学2011年11月实验一仿真头模型一、实验目的1、掌握MATLAB制图;2、绘制S-L头模型。二、实验内容1、用MATLAB绘制S-L头模型,该模型由十个位置、大小、方向、密度各异的椭圆组成。三、编程实现程序源代码:h=zeros(512,512);%ax0=256;y0=256;a=236;b=177;w=90/180*pi;fori=1:512forj=1:512c=b*b*((j-x0)*cos(w)+(i-y0)*sin(w))^2+a*a*(-(j-x0)*sin(w)-(i-y0)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=1.0;endendEnd%bx1=256;y1=261;a=224;b=170;w=90/180*pi;fori=1:512forj=1:512c=b*b*((j-x1)*cos(w)+(i-y1)*sin(w))^2+a*a*(-(j-x1)*sin(w)-(i-y1)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.7;endendend%cx2=312;y2=256;a=79;b=28;w=72/180*pi;fori=1:512forj=1:512c=b*b*((j-x2)*cos(w)+(i-y2)*sin(w))^2+a*a*(-(j-x2)*sin(w)-(i-y2)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.5;endendend%dx3=200;y3=256;a=105;b=41;w=108/180*pi;fori=1:512forj=1:512c=b*b*((j-x3)*cos(w)+(i-y3)*sin(w))^2+a*a*(-(j-x3)*sin(w)-(i-y3)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.5;endendend%ex4=256;y4=166;a=64;b=54;w=90/180*pi;fori=1:512forj=1:512c=b*b*((j-x4)*cos(w)+(i-y4)*sin(w))^2+a*a*(-(j-x4)*sin(w)-(i-y4)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.8;endendend%重叠部分x3=200;y3=256;a1=105;b1=41;w1=108/180*pi;fori=1:512forj=1:512c1=b1*b1*((j-x3)*cos(w1)+(i-y3)*sin(w1))^2+a1*a1*(-(j-x3)*sin(w1)-(i-y3)*cos(w1))^2;c=b*b*((j-x4)*cos(w)+(i-y4)*sin(w))^2+a*a*(-(j-x4)*sin(w)-(i-y4)*cos(w))^2;d=a*a*b*b;d1=a1*a1*b1*b1;if(c=d&c1=d1)h(i,j)=0.6;endendend%fx5=256;y5=230;a=12;b=12;w=0/180*pi;fori=1:512forj=1:512c=b*b*((j-x5)*cos(w)+(i-y5)*sin(w))^2+a*a*(-(j-x5)*sin(w)-(i-y5)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.8;endendend%重叠部分x4=256;y4=166;a1=64;b1=54;w1=90/180*pi;fori=1:512forj=1:512c1=b1*b1*((j-x4)*cos(w1)+(i-y4)*sin(w1))^2+a1*a1*(-(j-x4)*sin(w1)-(i-y4)*cos(w1))^2;c=b*b*((j-x5)*cos(w)+(i-y5)*sin(w))^2+a*a*(-(j-x5)*sin(w)-(i-y5)*cos(w))^2;d=a*a*b*b;d1=a1*a1*b1*b1;if(c=d&c1=d1)h(i,j)=0.9;endendend%gx6=256;y6=282;a=12;b=12;w=0/180*pi;fori=1:512forj=1:512c=b*b*((j-x6)*cos(w)+(i-y6)*sin(w))^2+a*a*(-(j-x6)*sin(w)-(i-y6)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.9;endendend%重叠部分x3=200;y3=256;a1=105;b1=41;w1=108/180*pi;fori=1:512forj=1:512c1=b1*b1*((j-x3)*cos(w1)+(i-y3)*sin(w1))^2+a1*a1*(-(j-x3)*sin(w1)-(i-y3)*cos(w1))^2;c=b*b*((j-x6)*cos(w)+(i-y6)*sin(w))^2+a*a*(-(j-x6)*sin(w)-(i-y6)*cos(w))^2;d=a*a*b*b;d1=a1*a1*b1*b1;if(c=d&c1=d1)h(i,j)=0.7;endendend%hx7=236;y7=411;a=12;b=6;w=0/180*pi;fori=1:512forj=1:512c=b*b*((j-x7)*cos(w)+(i-y7)*sin(w))^2+a*a*(-(j-x7)*sin(w)-(i-y7)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.8;endendend%ix8=256;y8=411;a=6;b=6;w=0/180*pi;fori=1:512forj=1:512c=b*b*((j-x8)*cos(w)+(i-y8)*sin(w))^2+a*a*(-(j-x8)*sin(w)-(i-y8)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.8;endendend%jx9=271;y9=411;a=12;b=6;w=90/180*pi;fori=1:512forj=1:512c=b*b*((j-x9)*cos(w)+(i-y9)*sin(w))^2+a*a*(-(j-x9)*sin(w)-(i-y9)*cos(w))^2;d=a*a*b*b;if(c=d)h(i,j)=0.8;endendend四、实验结果实验二仿真投影的实现一、实验目的将头模型进行投影得头模型投影图。二、实验内容要求每隔5度进行进行一次投影。三、编程实现程序源代码%存放0~90°投影数据g1=zeros(18,512);%存放90~180°投影数据g2=zeros(18,512);%存放0~180°投影数据g=zeros(36,512);%存放0~90°时,像素点在投影轴各个点的个数c1=zeros(18,512);%存放90~180°时,像素点在投影轴各个点的个数c2=zeros(18,512);%反投影矩阵h1=zeros(512,512);%反投影矩阵上得到反投影像素点在的个数c=zeros(512,512);%0~90°投影forw1=1:18w=(w1-1)*5/180*pi;fori=1:512forj=1:512if(w1=10)r1=(j-256.5)*cos(w)^2+(256.5-i)*sin(w)*cos(w)+256.500000000001;r=round(r1);elser1=(j-256.5)*cos(w)*sin(w)+(256.5-i)*sin(w)^2+256.500000000001;r=round(r1);endif(r0&r513)g1(w1,r)=g1(w1,r)+h(i,j);c1(w1,r)=c1(w1,r)+1;endendendend%90~180°投影forw1=1:18w=(w1-1)*5/180*pi;fori=1:512forj=1:512if(w1=10)r1=(-1)*(j-256.5)*cos(w)*sin(w)+(256.5-i)*cos(w)^2+256.500000000001;r=round(r1);elser1=(-1)*(j-256.5)*sin(w)^2+(256.5-i)*cos(w)*sin(w)+256.500000000001;r=round(r1);endif(r0&r513)g2(w1,r)=g2(w1,r)+h(i,j);c2(w1,r)=c1(w1,r)+1;endendendend%投影数据取平均值fori=1:18forj=1:512g1(i,j)=g1(i,j)/(c1(i,j)+1);g2(i,j)=g2(i,j)/(c2(i,j)+1);endendfori=1:18forj=1:512g(i,j)=g1(i,j);endendfori=19:36forj=1:512g(i,j)=g2(i-18,j);endend四、实验结果实验三仿真反投影的实现一、实验目的将所得的头模型投影图进行反投影得头模型反投影图。二、编程实现程序源代码:%0~90°数据反投影forw1=1:18w=(w1-1)*5/180*pi;fori=1:512forj=1:512if(w1=10)r1=(j-256.5)*cos(w)^2+(256.5-i)*sin(w)*cos(w)+256.500000000001;r=round(r1);elser1=(j-256.5)*cos(w)*sin(w)+(256.5-i)*sin(w)^2+256.500000000001;r=round(r1);endif(r0&r513)h1(i,j)=h1(i,j)+g1(w1,r);c(i,j)=c(i,j)+1;endendendend%90~180°数据反投影forw1=1:18w=(w1-1)*5/180*pi;fori=1:512forj=1:512if(w1=10)r1=(-1)*(j-256.5)*cos(w)*sin(w)+(256.5-i)*cos(w)^2+256.500000000001;r=round(r1);elser1=(-1)*(j-256.5)*sin(w)^2+(256.5-i)*cos(w)*sin(w)+256.500000000001;r=round(r1);endif(r0&r513)h1(i,j)=h1(i,j)+g2(w1,r);c(i,j)=c(i,j)+1;endendendend%反投影矩阵上像素点取平均值fori=1:512forj=1:512h1(i,j)=h1(i,j)/(c(i,j));endend三、实验结果