FBP.rar_FBP matlab _ct_ct复原_fbp算法_图像复原

%form the image E = [ 2.00 .9200 .6900 0 0 90 -.98 .8740 .6624 0 -.0184 90 -.62 .3100 .1100 .22 0 72 -.62 .4100 .1600 -.22 0 108 .21 .2500 .2100 0 .35 90 .21 .0460 .0460 0 .1 0 .21 .0460 .0460 0 -.1 0 .21 .0460 .0230 -.08 -.605 0 .21 .0230 .0230 0 -.605 0 .21 .0460 .0230 .06 -.605 90 3.01 .0750 .0500 0 -.480 0 ]; I=phantom(E,256); figure,imshow(I); IMG=double(I); THETA=0:179; % This MATLAB code takes an image matrix and vector of angles and then % finds the 1D projection (Radon transform) at each of the angles. It returns % a matrix whose columns are the projections at each angle. % pad the image with zeros so we don't lose anything when we rotate. [iLength, iWidth] = size(IMG); iDiag = sqrt(iLength^2 + iWidth^2); LengthPad = ceil(iDiag - iLength) + 2; WidthPad = ceil(iDiag - iWidth) + 2; padIMG = zeros(iLength+LengthPad, iWidth+WidthPad); padIMG(ceil(LengthPad/2):(ceil(LengthPad/2)+iLength-1), ... ceil(WidthPad/2):(ceil(WidthPad/2)+iWidth-1)) = IMG; % loop over the number of angles, rotate 90-theta (because we can easily sum % if we look at stuff from the top), and then add up. Don't perform any % interpolation on the rotating. n = length(THETA); PR = zeros(size(padIMG,2), n); for i = 1:n tic tmpimg = imrotate(padIMG, 90-THETA(i), 'bilinear', 'crop'); PR(:,i) = (sum(tmpimg))'; THETA(i) toc end n = size(PR,1); sideSize = n; % filter the projections a = 1; [Length, Count] = size(PR); w = [-pi:(2*pi)/Length:pi-(2*pi)/Length]; rn1 = abs(2/a*sin(a.*w./2)); rn2 = sin(a.*w./2); rd = (a*w)./2; r = rn1*(rn2/rd)^2; f = fftshift(r); for i = 1:Count IMG = fft(PR(:,i)); fimg = IMG.*f'; g(:,i) = ifft(fimg); end g = real(g); filtPR=g; % convert THETA to radians th = (pi/180)*THETA; % set up the image m = length(THETA); BPI = zeros(sideSize,sideSize); % find the middle index of the projections midindex = (n+1)/2; % set up x and y matrices x = 1:sideSize; y = 1:sideSize; [X,Y] = meshgrid(x,y); xpr = X - (sideSize+1)/2; ypr = Y - (sideSize+1)/2; % loop over each projection for i = 1:m tic disp(['On angle ', num2str(THETA(i))]); filtIndex = round(midindex + xpr*sin(th(i)) - ypr*cos(th(i))); BPIa = zeros(sideSize,sideSize); spota = find((filtIndex > 0) & (filtIndex <= n)); newfiltIndex = filtIndex(spota); BPIa(spota) =filtPR(newfiltIndex(:),i); BPI = BPI + BPIa; toc end BPI = BPI./m; figure,imshow(BPI,[]);