平行束投影重建算法（MATLAB程序）——DBP、FBP、CBP

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