function [] = CBP(projData, filter)
% Fuction [] = CBP(projData，filter)
% projData is the projection data.
% filter is R-L Filter or S-L Filter. "filter = 1;" -- R-L & "filter = 2;"-- S--L
% Aim : 2D Convoluted Back-Projection
% Result: draw the sinogram and the restruction image
% Figure： 1111 & 2222
% Draw the sinogram
figure(111111);imshow(projData,[]);title('Sinogram');
% Restruction -- Convoluted Back-Projection
rowNum = size(projData,1);
columnNum = size(projData,2);
imageSzie = pow2(nextpow2(size(projData,1))-1);
% Selecte the style of Filter
switch filter
case 1
% R-L Filter
filterCore = zeros((rowNum*2-1),1);
for i = 0:(rowNum-1)
if i == 0
filterCore(rowNum) = 1/4;
elseif rem(i,2) == 0
filterCore(rowNum-i) = 0;
filterCore(rowNum+i) = 0;
else
filterCore(rowNum-i) = -1/(i*pi)^2;
filterCore(rowNum+i) = -1/(i*pi)^2;
end
end
case 2
% S-L Filter
filterCore = -2/pi^2./(4*(-(rowNum-1):(rowNum-1)).^2-1);
otherwise
filterCore = 0; % 有待完善
end
%Convolution
projCon = zeros(rowNum,columnNum);
for i = 1:columnNum
s = projData(:,i);
dataCon = conv(s',filterCore');
projCon(:,i) = dataCon(rowNum:2*rowNum-1);
end
% Back-Projection
imageRes = zeros(imageSzie,imageSzie);
for i = 1:imageSzie
for j = 1:imageSzie
for k = 0:179
theta = k/180*pi;
t=(j-imageSzie/2-0.5)*cos(theta)+(imageSzie/2+0.5-i)*sin(theta)+(rowNum+1)/2;
if(t < rowNum && t >= 1)
t1 = floor(t);
t2 = floor(t+1);
imageRes(i,j) = imageRes(i,j)+(t2-t)*projCon(t1,k+1)+(t-t1)*projCon(t2,k+1);
end
end
end
end
imageRes = pi/columnNum*imageRes;
figure(222222);imshow(imageRes,[]);title('Restruction Image(CBP)');
