%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% EE 684 - Spring 2016
%%%% Simple Code for Image Manipulation 
%%%%
%%%% Your project should try to 
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%(A)****************************************************************

% Reading and reconstructing MRI image from k-space
fileToRead1= 'MRI_kspace_data.txt';
rawdata=  importdata(fileToRead1);
ksp= reshape(rawdata(:,1)+j*rawdata(:,2), 256, 256);    % compose complex k-space

% display k-space
figure(1); 
image(abs(ksp)*1e6);    % different ways of displaying an image
colormap(gray);         % needed to display in grayscale not color

% reconstruction from both magnitude and phase of k-space
im= fftshift(ifft2(fftshift(ksp)));     % inverse 2D FFT

% display image 
figure(2); 
imagesc(abs(im));       % remember that image is complex
%imagesc(real(im));
%imagesc(imag(im));
colormap(gray);         % needed to display in grayscale not color


%(B)****************************************************************

% reconstruction from only magnitude of k-space
im_mag= fftshift(ifft2(fftshift(abs(ksp))));     % inverse 2D FFT

% display image 
figure(3); 
imagesc(abs(im_mag));       % remember that image is complex
%image (abs(im_mag)*10e7)
colormap(gray);         % needed to display in grayscale not color

% reconstruction from only phase of k-space
im_ph= fftshift(ifft2(fftshift(exp(j*angle(ksp)))));     % inverse 2D FFT

% display image 
figure(4); 
imagesc(abs(im_ph));       % remember that image is complex
colormap(gray);         % needed to display in grayscale not color


%(C)****************************************************************

% Note that the image reconstructed above has real numbers
% Now we will quantize that to different levels (number of bits)
% We should compare visual appearance at different level
nBits= 8;

im_min= min(min(abs(im)));
im_max= max(max(abs(im)));
im_quantized= round((abs(im)-im_min)/(im_max-im_min) * (2^nBits - 1));
figure (5);
imagesc(im_quantized);       % remember that image is complex
colormap(gray);         % needed to display in grayscale not color


%(D)****************************************************************


% Effect of k-space sampling 
undersampling_factor= 2;
size1= size(ksp);               % detect size of k-space data
%ksp_undersampled= ksp(1:undersampling_factor:size1(1), 1:undersampling_factor:size1(2));  % undersample in both kx and ky
%ksp_undersampled= ksp(1:undersampling_factor:size1(1), :);          % undersample only in kx
ksp_undersampled= ksp(:, 1:undersampling_factor:size1(2));          % undersample only in ky

im_undersampled= fftshift(ifft2(fftshift(ksp_undersampled),size1(1), size1(2)));     % inverse 2D FFT

% display image 
figure(6); 
imagesc(abs(im_undersampled));       % remember that image is complex
colormap(gray);         % needed to display in grayscale not color

% -----

% Effect of k-space bandwidth
BW_factor= 2;
size1= size(ksp);               % detect size of k-space data
BW_kx= round(size1(1) / BW_factor);
BW_ky= round(size1(2) / BW_factor);
%ksp_redBW= ksp( size1(1)/2-BW_kx/2+1 : size1(1)/2+BW_kx/2, size1(2)/2-BW_ky/2+1 : size1(2)/2+BW_ky/2-1);  % reduce BW in both kx and ky
%ksp_redBW= ksp( size1(1)/2-BW_kx/2 : size1(1)/2+BW_kx/2-1, :);           % reduce BW only in kx
ksp_redBW= ksp( : , size1(2)/2-BW_ky/2 : size1(2)/2+BW_ky/2-1);          % reduce BW only in ky

im_redBW= fftshift(ifft2(fftshift(ksp_redBW),size1(1), size1(2)));     % inverse 2D FFT

% display image 
figure(7); 
imagesc(abs(im_redBW));       % remember that image is complex
colormap(gray);         % needed to display in grayscale not color


%(E)****************************************************************

%% Read in CT image (DICOM format)
[CT, MAP, ALPHA, OVERLAYS]= dicomread('CT_000043.dcm');

% display image 
figure(7); 
imagesc(CT);
%imtool(CT)
info = dicominfo('CT_000043.dcm')           % display DICOM image info


%(F)****************************************************************

% Use imtool to check how enlargement makes image blocky
% Here, we try Fourier interpolation for better results
size_factor= 4;
size1= size(CT);
%fCT= fft2(double(CT));
fCT= zeros(size1(1)*size_factor , size1(2)*size_factor);
fCT(1:size1(1), 1:size1(2))= fftshift(fft2(fftshift(double(CT+1024))));;

CT_interp= abs(fftshift(ifft2(fftshift(fCT))))-1024; 
figure(8);
imagesc(CT_interp);
colormap(gray);


%(G)****************************************************************



%(H)****************************************************************


%(I)****************************************************************



