两种扇形束CT投影生成方法



在学习CT重建算法时最先遇到的问题是如何建立自己的投影数据对算法进行测试，本文总结列两种常用的扇形束CT投影生成方法。

首先是插值法

1、插值法假定原始断层图像 f（x,y） 是连续的，非离散的。

2、对于f（x,y）中不在整数点的位置则用插值法补齐，使之连续。

3、用L（w,s）表示光源与y轴夹角为w时，第s个探测器的投影射线。

4、对 L （w,s）以0.5（此值越小，精度越高）为间隔进行取样，当取样点不在整点时，用插值法确定其数值

5、对L（w,s）上所有取样点进行加和。

射线法

1、射线法吧断层图像f（x,y）看成是一个一个离散的像素块。

2、同样有射线 L（w,s）表示光源与y轴夹角为w时，第s个探测器的投影射线。

3、以L 经过各个像素块的长度为权值，对其灰度值进行积分。

下面是两种方法的matlab代码

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

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

nDetectors=512;%探测器成像矩阵512*512me
nViews=360;%投影个数
projection=zeros(nViews,nDetectors);


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

stepBeam = 2 * pi/nViews;                %旋转步进的角度
focalDistance_phy=40;%射线源到旋转中心的实际距离me
detecDistance_phy=40;
focalDistance_dig=focalDistance_phy*phyRatoDig;%像素距离me 512
detecDistance_dig=detecDistance_phy*phyRatoDig;%512



detecLength=41.3;%检测器大小me
detecLength_dig=detecLength*phyRatoDig;%528.64
detecResolution=512;%检测器分辨率
unitDis=detecLength/(detecResolution-1);%单位像素的大小me 0.0808
yDetector=([1:detecResolution]-1)*unitDis-detecLength/2;%s 射线与探测器交点坐标[-1/2:1/512:1/2]*41.3

gamma= atan(yDetector/(focalDistance_phy+detecDistance_phy));%探测器角度r me



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;%中心射线与y轴夹角
    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)));%%%%%%%%%projection.m%%%%%%%%%%%%

clear;
clear all;
M=256;
N=256;
I = phantom(M,N);
figure(1);
imshow(I);
tempx = -N/2:N/2-1;
tempx = repmat(tempx,M,1);

tempy = -M/2:M/2-1;
tempy = repmat(tempy',1,N);

nDetectors = 512;%探测器成像矩阵512*512me
nViews=360;%投影个数
projections = zeros(nViews,nDetectors);


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

stepBeam = 2 * pi/nViews;                %旋转步进的角度
focalDistance_phy=40;%射线源到旋转中心的实际距离me
detecDistance_phy=40;
focalDistance_dig=focalDistance_phy*phyRatoDig;%像素距离me 512
detecDistance_dig=detecDistance_phy*phyRatoDig;%512


detecLength=41.3;%检测器大小me
detecLength_dig=detecLength*phyRatoDig;%528.64
detecResolution=512;%检测器分辨率
unitDis=detecLength/(detecResolution-1);%单位像素的大小me 0.0808
yDetector=([1:detecResolution]-1)*unitDis-detecLength/2;%s 射线与探测器交点坐标[-1/2:1/512:1/2]*41.3

gamma= atan(yDetector/(focalDistance_phy+detecDistance_phy));%探测器角度r me



sample=0.5;                                %沿射线进行步进的步长，即采样的间隔
sample_times=floor(sqrt((detecDistance_dig+focalDistance_dig)^2+(detecLength_dig/2)^2)/sample);
disp('projection');


for fanNum=1:nViews;
    beta=(fanNum-1)*stepBeam;%中心射线与y轴夹角
    sourceX = focalDistance_dig * cos(beta + pi/2);           %在极坐标情况下计算射线源的坐标（可包含负值）
    sourceY = focalDistance_dig * sin(beta + pi/2);
    source = [sourceX sourceY];
   	for  raynum=1: nDetectors
         value=0; 
         gamaTemp =gamma(1,raynum);
         full_angle=beta+gamaTemp;%投影角
         oned = [cos(full_angle - pi/2) sin(full_angle - pi/2)]; 
         shot = source - sum(source.*oned)*oned;
         oneshot = shot/norm(shot);
         d = norm(shot) - tempx*oneshot(1) + tempy*oneshot(2);
         ind = find(abs(d)<0.5);
         t = d(ind);
         nd = zeros(size(d));
         nd(ind) = 2*(0.25-t.*t);
         dd = I.*nd;
         projections(fanNum,raynum) = sum(dd(:));
    end
end