CT图像重建的MATLAB代码

clear % 原始图像参数 pic_wide=20;%原始图像长宽 % 椭圆参数 xc=2;yc=-1;a=4;b=5;alpha=0;color=0.8; xc1=-4;yc1=5;a1=3;b1=2;alpha1=pi/4;color1=0.5; % (xc,yc) is the position of center % a is long radio % b is short radio % alpha is the position of the angle % color is the color of the ellipse, 0=empty,1=black % 探测器参数 point_size=64;%探测点数目 angle_size=360;%角度变化次数（从0度到359度，步进1度，共360个角度值) gamma_max=pi/6;%探测点最大张角 delta_gamma=gamma_max/point_size;%各探测点之间的张角 d=40;%射线源到中心的距离 % 反投影图像参数 pic_wide_new=64;% 重建图像的长宽 %画出原始图像 subplot(3,2,1) hold on %确定原始图像长宽并画出黑色背景 plot(-pic_wide/2,-pic_wide/2) plot(pic_wide/2,-pic_wide/2) plot(-pic_wide/2,pic_wide/2) plot(pic_wide/2,pic_wide/2) bak1=[-pic_wide/2,pic_wide/2,pic_wide/2,-pic_wide/2]; bak2=[-pic_wide/2,-pic_wide/2,pic_wide/2,pic_wide/2]; fill(bak1,bak2,[0,0,0]) line([-pic_wide/2,pic_wide/2],[0,0]) line([0,0],[-pic_wide/2,pic_wide/2]) axis square my_ellipse(xc,yc,a,b,alpha,color) my_ellipse(xc1,yc1,a1,b1,alpha1,color1) title('原始图像') % 计算投影 projection=zeros(point_size,angle_size); for k=1:point_size for kk=1:angle_size gamma=delta_gamma*(k-point_size/2-1/2); dd=d*sin(gamma);%计算射线到中心的距离 projection(k,kk)=my_projection(xc,yc,a,b,alpha,color,dd,(kk-1)/180*pi); projection(k,kk)=projection(k,kk)+my_projection(xc1,yc1,a1,b1,alpha1,color1,dd,(kk-1)/180*pi); end end % 中间探测器处投影随角度的变化 subplot(3,2,2) axis auto hold on plot(0,0) plot(0:angle_size-1,projection(point_size/2,:),'g.-') title('中间探测器处投影随角度的变化(横向曲线)') % 初始位置处各探测器的投影 subplot(3,2,3) axis auto hold on plot(1:point_size,projection(:,1),'r.-') title('初始位置处各探测器的投影(纵向曲线)') % 投影二维图像 subplot(3,2,4) axis square hold on projection_pic=zeros(point_size,angle_size,3); projection_pic(:,:,1)=projection(:,:)/2/pic_wide*255; projection_pic(:,:,2)=projection(:,:)/2/pic_wide*255; projection_pic(:,:,3)=projection(:,:)/2/pic_wide*255; projection_pic=uint8(projection_pic); imshow(projection_pic) title('投影数据二维图像') imwrite(projection_pic,'投影数据.bmp','bmp'); % 滤波函数 subplot(3,2,5) axis auto hold on g=-point_size:point_size; h=-point_size:point_size; for kk=-point_size:point_size g(kk+point_size+1)=g_myfunction(kk,delta_gamma);%用于反投影 h(kk+point_size+1)=h_myfunction(kk,delta_gamma);%用于显示 end plot(-point_size:point_size,g,'b.-') plot(-point_size:point_size,h,'r.-') title('RL滤波函数') % 计算反投影 pic_data=zeros(pic_wide_new,pic_wide_new);%储存像素的重建结果 for k=1:pic_wide_new for kk=1:pic_wide_new sum=0;%积分初始值为0 %第一步，计算图像点坐标 x=(k-1/2-pic_wide_new/2)/pic_wide_new*pic_wide; y=(kk-1/2-pic_wide_new/2)/pic_wide_new*pic_wide; %第二步，转换成极坐标 [r,thet]=trans(x,y); for kkk=1:angle_size %角度从0度到359度 beta=(kkk-1)/angle_size*2*pi; %下面开始对固定的角度beta进行卷积计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %第三步，计算gamma_0，即理想的射线方向 gamma_0=asin(r*cos(thet-beta)/sqrt(d^2+r^2+2*d*r*sin(beta-thet))); %第四步，计算最接近的采样点（不同采样点对应不同的方向） num=gamma_0/delta_gamma+point_size/2+1/2;% gamma=delta_gamma*(k-point_size/2-1/2)的逆运算 num1=floor(num);%最接近的采样点（较小的编号） if num1<=0 num1=1; else if num1>=point_size num1=point_size-1; else end end num2=num1+1;%最接近的采样点（较大的编号） gamma1=(num1-point_size/2-1/2)*delta_gamma; gamma2=(num2-point_size/2-1/2)*delta_gamma; %第五步，计算卷积(对gamma积分) conv1=my_conv(num1,delta_gamma,gamma1,point_size,beta,d,angle_size,projection(:,kkk),g);%参数传送只有beta角对应的投影值 conv2=my_conv(num2,delta_gamma,gamma2,point_size,beta,d,angle_size,projection(:,kkk),g);%参数传送只有beta角对应的投影值 %第六步，对2个卷积结果插值 conv=conv1*(num2-num)+conv2*(num-num1); %结束对固定的角度beta进行的卷积计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %第七步，beta角循环360度,对beta进行积分 L2=d^2+r^2+2*d*r*sin(beta-thet); sum=sum+conv/L2; end %得到积分结果 pic_data(k,kk)=sum*(2*pi/angle_size); end end % 反投影成像 subplot(3,2,6) axis square hold on new_pic=zeros(pic_wide_new,pic_wide_new,3); %下面进行坐标变换，以适应图像的坐标规则 pic_data=pic_data'; pic_data=flipud(pic_data); %变换完成 new_pic(:,:,1)=pic_data(:,:)*255; new_pic(:,:,2)=pic_data(:,:)*255; new_pic(:,:,3)=pic_data(:,:)*255; new_pic=uint8(new_pic); imshow(new_pic) title('反投影重建图像') imwrite(new_pic,'重建图像.bmp','bmp');