1.算法描述
CT重建算法大致分為解析重建算法和迭代重建算法烦秩，隨著CT技術的發(fā)展，重建算法也變得多種多樣郎仆，各有各的有特點只祠。本文使用目前應用最廣泛的重建算法——濾波反投影算法（FBP）作為模型的基礎算法。FBP算法是在傅立葉變換理論基礎之上的一種空域處理技術扰肌。它的特點是在反投影前將每一個采集投影角度下的投影進行卷積處理抛寝，從而改善點擴散函數(shù)引起的形狀偽影，重建的圖像質量較好曙旭。
上圖應可以清晰的描述傅立葉中心切片定理的過程：對投影的一維傅立葉變換等效于對原圖像進行二維的傅立葉變換盗舰。傅立葉切片定理的意義在于，通過投影上執(zhí)行傅立葉變換桂躏，可以從每個投影中得到二維傅立葉變換钻趋。從而投影圖像重建的問題，可以按以下方法進行求解：采集不同時間下足夠多的投影（一般為180次采集）剂习，求解各個投影的一維傅立葉變換蛮位，將上述切片匯集成圖像的二維傅立葉變換，再利用傅立葉反變換求得重建圖像鳞绕。
投影重建的過程是失仁，先把投影由線陣探測器上獲得的投影數(shù)據(jù)進行一次一維傅立葉變換，再與濾波器函數(shù)進行卷積運算们何，得到各個方向卷積濾波后的投影數(shù)據(jù)萄焦；然后把它們沿各個方向進行反投影，即按其原路徑平均分配到每一矩陣單元上冤竹，進行重疊后得到每一矩陣單元的CT值拂封；再經(jīng)過適當處理后得到被掃描物體的斷層圖像
算法步驟如下：
1. 將原始投影進行一次一維傅立葉變換
2. 設計合適的濾波器，在φ_i的角度下將得到原始投影p(x_r,φ_i)進行卷積濾波贴见，得到濾波后的投影烘苹。
3. 將濾波后的投影進行反投影躲株，得到滿足x_r=r cos?（(θ - φ_i)）方向上的原圖像的密度片部。
4. 將所有反投影進行疊加，得到重建后的投影。
2.仿真效果預覽
matlab2013B仿真結果如下：
3.MATLAB部分代碼預覽
coordinateSource=[];
projMatrix=[];
detector=[];
proj=load(char('projection.mat'));
phyRatoDig=proj.phyRatoDig;
projMatrix=proj.projection;
yDetector=proj.yDetector;
nDetectors=proj.nDetectors;
figure(2)
showimge(projMatrix,360,512,0,max(max(projMatrix)));
D_dig=proj.focalDistance_dig;
sourceToDetector_dig=proj.focalDistance_dig+proj.detecDistance_dig;
s=[];
s=D_dig/sourceToDetector_dig*yDetector(1,:)*phyRatoDig;
Detector=yDetector(1,:)*phyRatoDig;
%
pe=[];
M=D_dig./sqrt(D_dig.^2+s.^2);
nViews=proj.nViews;
%
for i=1:nViews
pe(i,:)=projMatrix(i,:).*M;
end
figure(3);
showimge(pe,360,512,0,max(max(pe)));
disp('Filtering')
filternum=128;
filter_ramp=zeros(filternum,1);
for j=1:filternum ??% 16 point ramp filter
i=j-1-filternum/2;
if(i==0)
filter_ramp(j,1)=1/(8.0);
elseif (mod(i,2)==0)
filter_ramp(j,1)=0;
elseif (mod(i,2)==1)
filter_ramp(j,1)=-0.5/(i*i*pi*pi);
end
end
m=1;
figure(4);
plot(filter_ramp);
pfilter=[];
length_conv=filternum+nDetectors-1;
pPro=zeros(length_conv,1);
temp_pro=zeros(nDetectors,1);
h_filter=filternum/2;
ii=length_conv-h_filter-nDetectors;
for s=1:nViews % sample-loop
% for ?pp=1:h_filter
pro_left =(pe(s,1)+pe(s,2))/2.0;
pro_right=(pe(s,nDetectors)+pe(s,nDetectors-1))/2.0;
for pp=1:h_filter; ??????????????%left part
pPro(pp,1)=pro_left; ???????
end
% ???
for pp=1:nDetectors ?????????????????????%middle part
pPro(h_filter+pp,1)=pe(s,pp);
end
% ???
for pp=h_filter+nDetectors+1:length_conv
pPro(pp)=pro_right;
end
% ??result_conv ???
for n=1:nDetectors
result_conv=0;
??? for jj=1:filternum
pPmove=pPro(n+jj-1,1);
result_conv=result_conv+pPmove*filter_ramp(jj,1);
end
pfilter(s,n)=result_conv;
end
end
figure(5);
showimge(pfilter,360,512,min(min(pfilter)),max(max(pfilter)));
% %%%%%%%%%%%%%%%%%%%%%%reArrange%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%Back projection%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('BackProjection')
detecLength=proj.detecLength;
unitDis=detecLength/(nDetectors-1);
unitDis_dig=unitDis*phyRatoDig;
deltaBeta=2*pi/nViews;
M=proj.M;
N=proj.N;
fReconstruct=[];
for i=1:M
x=i-(M+1)/2;
for j=1:N
y=(N+1)/2-j;
r=sqrt(x^2+y^2);
theta=atan2(y,x);
% ????
result=0;
for s=1:nViews
beta= (s-1)*deltaBeta;
s1=D_dig*r*cos(beta-theta)/(D_dig+r*sin(beta-theta));
U=(D_dig+r*sin(beta-theta))/D_dig;
p1=sourceToDetector_dig/D_dig*s1;
if(p1>Detector(1,1)&&p1<Detector(1,512))
num=(p1-Detector(1,1)+unitDis_dig)/unitDis_dig; ???????
numlow=floor(num);
result=result+((num-numlow)*pfilter(s,numlow)+(1-num+numlow)*pfilter(s,numlow+1))/U/U*deltaBeta;
end
end
fReconstruct(i,j)=result;
if( fReconstruct(i,j)<0)
fReconstruct(i,j)=0;
end
end
end
%
figure(6)
%
final=zeros(M,N);
for i=1:M
final(i,:)=fReconstruct(:,257-i);
end
A_022
