%%%%%%%%%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)));

————————————————

版权声明：本文为CSDN博主「fpga&amp;matlab」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。

原文链接：https://blog.csdn.net/ccsss22/article/details/115274040