对滤波反投影重建算法的研究以phantom图为例，构建滤波器，重建图像

%%%%%%%%%projection.m%%%%%%%%%%%%

clear;
clear all;
M=256;
N=256;
I = phantom(M,N);
figure(1);
imshow(I);

nDetectors=512;
nViews=360;
projection=zeros(nViews,nDetectors);

length=20;
weight=20;
phyRatoDig=M/length;

stepBeam = 2 * pi/nViews;                %旋转步进的角度
focalDistance_phy=40;
detecDistance_phy=40;
focalDistance_dig=focalDistance_phy*phyRatoDig;
detecDistance_dig=detecDistance_phy*phyRatoDig;
diameter=sqrt(M*N+M*N);

detecLength=41.3;
detecLength_dig=detecLength*phyRatoDig;
detecResolution=512;
unitDis=detecLength/(detecResolution-1);
yDetector=([1:detecResolution]-1)*unitDis-detecLength/2;

gamma= atan(yDetector/(focalDistance_phy+detecDistance_phy));

sample=0.5;                                %沿射线进行步进的步长，即采样的间隔
sample_times=floor(sqrt((detecDistance_dig+focalDistance_dig)^2+(detecLength_dig/2)^2)/sample);
disp('projection')
for fanNum=1:nViews;
    if(rem(fanNum,10)==0)
        disp('Current Sample')
        fanNum
    end
    beta=(fanNum-1)*stepBeam;
    sourceX = focalDistance_dig * cos(pi/2 + beta);           %在极坐标情况下计算射线源的坐标（可包含负值）
    sourceY = focalDistance_dig * sin(pi/2 + beta);
           for  raynum=1: nDetectors
             deltaX=0;
             deltaY=0;
             value=0; 
             gamaTemp =gamma(1,raynum);
             full_angle=beta+gamaTemp;
                 if(full_angle<0)
                    full_angle=full_angle+2*pi;
                end
                if(full_angle>2*pi) 
                    full_angle=full_angle-2*pi;
                end
                if(full_angle >= 0 && full_angle < pi/2)
                
                    flagX = 1;
                    flagY = -1;
                else
                    if(full_angle >= pi/2 && full_angle < pi)
                   
                    flagX = 1;
                    flagY = 1;
                    
                    else
                        if(full_angle >= pi && full_angle < 3 * pi/2)
                
                     flagX = -1;
                    flagY = 1;
                    
                        else
                            if(full_angle >= 3 * pi/2 && full_angle < 2 * pi)
                
                           flagX = -1;
                           flagY = -1;
                            end
                        end
                    end
                end
                
                x_delta=abs(sin(full_angle));
                 y_delta=abs(cos(full_angle));
                 deltaX=flagX*x_delta;                                      %综合求得x和y方向的增长情况
                 deltaY=flagY*y_delta;
                 for  k=1:sample_times 
                      pixel = 0;                               %该射线路径上的一点投影值
                      x=sourceX+k*deltaX*sample+N/2;
                     y=sourceY+k*deltaY*sample+N/2;
                    if (x>=1&&x<=M&&y>=1&&y<=N)
                     
                         ix = floor(x);
                         x = x - ix;
                         iy = floor(y);
                         y = y - iy;
                        if((ix + 1 + (iy + 1) * N) >= (M * N))
                            pixel = I(iy,ix);
                        else
                            pixel=x*(I(iy,ix+1)-I(iy,ix))+y*(I(iy + 1,ix)-I(iy,ix))+(I(iy + 1,ix + 1)+I(iy,ix) -I(iy + 1,ix ) -I(iy + 1,ix + 1 ))*x*y+I(iy,ix);
                        end
                      
                    end
                      value=value+pixel;

             
        end
        projection(fanNum,raynum)=value;
        end        
end
showimge(projection,360,512,min(min(projection)),max(max(projection)));

D-022