angs = 0:179; % angles of rotation 0, 1, 2...
359
%init_img = phantom('Modified Shepp-Logan', 180); % Initial image 2D [100 x 100]
%init_img = phantom(256); % Initial image 2D [100 x 100]
init_img=im2gray(imread('[Link]'));
%init_img = phantom(180);
%sino = radon(init_img, 0:179);
%sino = radon(init_img, angs); % Create a sinogram using radon transform
sino=im2gray(imread('[Link]'));
sino_org=uint8(sino);
%sino=sino';
% Here is your function ....
% This step is necessary so that frequency information is preserved: we pad
% the sinogram with zeros so that this is ensured.
sino = padarray(sino, floor(size(sino,1)/2), 'both');
% Rotate the sinogram 90-degree to be compatible with your codes definition of view
and radial positions
% dim 1 -> view
% dim 2 -> Radial position
sino = sino';
% diagDim should be the length of an individual row of the sinogram - don't
% hardcode this!
diagDim = size(sino, 2);
% The 2DFT (2D Fourier transform) of our image will start as a matrix of
% all zeros.
fimg = zeros(diagDim);
% Design your 1-d ramp filter.
rampFilter_1d = abs(linspace(-1, 1, diagDim))';
rowIndex = 1;
%for nn = 80:80
for nn = angs
% fprintf('rowIndex = %g => nn = %g\n', rowIndex, nn);
% Each contribution to the image's 2DFT will also begin as all zero.
imContrib = zeros(diagDim);
% Get the current row of the sinogram - use rowIndex.
curRow = sino(rowIndex,:);
% Take the 1D Fourier transform the current row - be careful, as it's
% necessary to perform ifftshift and fftshift as Matlab tends to
% place zero-frequency components of a spectrum at the edges.
fourierCurRow = fftshift(fft(ifftshift(curRow)));
% Place the Fourier-transformed sino row and place it at the center of
% the next image contribution. Add the ramp filter in Fourier domain.
imContrib(round(diagDim/2), :) = fourierCurRow;
imContrib(round(diagDim/2), :) = imContrib(round(diagDim/2), :)' .*
rampFilter_1d; % <-- NOT fft(rampFilter_1d)
% Rotate the current image contribution to be at the correct angle on
% the 2D Fourier-space image.
imContrib = imrotate(imContrib, nn, 'crop');
% Add the current image contribution to the running representation of
% the image in Fourier space!
fimg = fimg + imContrib;
rowIndex = rowIndex + 1;
end
% Finally, just take the inverse 2D Fourier transform of the image! Don't
% forget - you may need an fftshift or ifftshift here.
rcon = fftshift(ifft2(ifftshift(fimg)));
%rcon = fftshift(ifft2(fimg));
%rcon=rcon';
%rcon=rcon';
%rcon=flip(rcon,1);
rcon=imrotate(rcon,180);
icol=3;
irow=2;
I3 = iradon(double(sino_org),0:179,'linear','Ram-Lak',1,output_size);
%,'linear','none');
nimage=1;
figure,subplot(irow,icol,nimage),imshow(sino_org);title('Original
Image');nimage=nimage+1;
subplot(irow,icol,nimage),imshow(uint8(sino));title('Binary');nimage=nimage+1;
subplot(irow,icol,nimage),imshow(real(rcon));title('BP');nimage=nimage+1;
subplot(irow,icol,nimage),imshow(I3);title('iRadon');nimage=nimage+1;
subplot(irow,icol,nimage),imshow(init_img);title('init_img');nimage=nimage+1;
subplot(irow,icol,nimage),imshow(fimg);title('fimg');nimage=nimage+1;
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Sinogram by IK', 'NumberTitle', 'Off')
