Reconstructions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Author: Onkar Raut

%Title: Program to reconstruct the image using the Central Slice Theorem

%Purpose: Medical Imaging Final Exam

%Date: May 05th,2010

%Statement: The central slice theorem states that the Fourier Transform of

%the projection of an image at an angle (projection of image is the RADON 

%Transform) is nothing but the set of points lying on a line passing 

%through the origin(slice) in the 2D fourier transform of that image

%Objective:      To verify and reconstruct medical data using the central slice theorem. 

%                     The following is a very practical solution to perfectly reconstruct an image from the FT of the image

%                     We Shall do the same by doing the following:

%                    1. acquire the projections of the image using the radon

%                   transform and then use them to reconstruct the image.

%                    2. Take the fourier transform of the profiles/prjections,

%                   modulate them circular shift them and align about the origin

%                    3. interpolate the 2D fft of the image and take the 2D inverse

%                   fourier transform. 

%To Do: The code only needs get the proper

%            modulation function such that the image

%            comes into the center. I would be able to do

%            this in due time. Will update this

%            and post as and when I find the appropriate

%            modulation function.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

close all;

clc;

imsize=256;

imcenter=imsize/2;

offset=32;

angle_change=1;

del_phi=angle_change*pi/180;

% blank canvas

I=zeros(imsize,imsize);

% Creating image

I(imcenter-offset:imcenter+offset,imcenter-offset:imcenter+offset)=1;

I=imread('phantom256.png');

% Getting the profiles/projections

profile=radon(I,0:angle_change:180-angle_change);

% Prepare the profiles for reconstruction;

req_size=size(profile,1)-imsize;

new_profile=profile(ceil(req_size/2):end-ceil(req_size/2),:);

profile=new_profile;

% Taking the fft of the profiles

profile_fft=fft(profile);

% Circular shifting the profile

profile_fft=circshift(profile_fft,round(size(profile_fft,1)/2));

% define the lengths of the profiles i.e. zero paddings and store them as a

% vector

angle_vector=0:del_phi:45*pi/180;

profile_len=imsize./cos(angle_vector);

m=fliplr(profile_len);

profile_len_90=[profile_len m(2:end)];

profile_lengths=round(profile_len_90);

n=fliplr(profile_lengths);

lengths=[profile_lengths n(2:end-1)];

prof_lengths=lengths';

im=zeros(imsize,imsize);

for j=0:(180-angle_change/angle_change)

    theta=pi/180*j*angle_change;

    req_prof_length=prof_lengths(j+1,1);

    curr_prof_length=size(profile_fft(:,j+1),1);

    diff=req_prof_length-curr_prof_length;

    % Get the zero padded profile

    curr_profile=[zeros(floor(diff/2),1); profile(:,j+1); zeros(floor(diff/2),1)];

    % Get FFT of the zero padded profile

    current_profile_fft_int=fft(curr_profile);

%    current_profile_fft=fft(curr_profile);

%   Need Proper Modulation here before the circshift!!!!

    mod=ones(size(current_profile_fft_int,1),size(current_profile_fft_int,2));

    mod(1:2:end)=-1;

    %   Also need to use ram lak filter here

    current_profile_fft_int=mod.*current_profile_fft_int;

    current_profile_fft=circshift(current_profile_fft_int,round(size(current_profile_fft_int,1)/2));

    %define the rotation matrix

    rotate_matrix=[cos(theta) -sin(theta);sin(theta) cos(theta)];

    % find all the applicable rotated points

    dummy_coor=zeros(2,size(current_profile_fft,1)-1);

    dummy_coor(1,:)=(-floor(size(current_profile_fft,1)/2-1):floor(size(current_profile_fft,1)/2)-1);

    dummy_coor(2,:)=0;

    rot_coor=ceil(rotate_matrix*dummy_coor);

    true_coor=rot_coor+imcenter;

    for i=1:size(rot_coor,2);

       % if(true_coor(2,i)<imsize&&true_coor(1,i)<imsize)

            im(true_coor(1,i),true_coor(2,i))=current_profile_fft(i,1); % im(true_coor(1,i),true_coor(2,i))+

       % end

    end

end

bh=circshift(im,[128 128]);

%bh=imrotate(bh,45);

figure(1),imshow(20*log(abs(im)),[]);

final_uninterpolated_image=ifft2(bh);

figure(2),imshow(20*log(abs(circshift(fft2(I),[imcenter imcenter]))),[]);

figure(3),imshow(abs(final_uninterpolated_image),[]);

shift_im=circshift(final_uninterpolated_image,[64 64]);

figure(4), imshow(abs(shift_im),[]);