%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% P_FILE.m
%
% Matlab program to take praw_re.dat and praw_im.dat files obtained
% from Bud Wendt and create phase and magnitude MRI images
%
% Run P_FILE at the prompt. The four images are saved as 
% Aimagdata.bmp, Aimagimage.bmp,  Amagnitude.bmp, and Aphase.bmp 
%
% Author: Orlando Cardoso
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;  % clear all variables

% Setup input display image size, colormap, and Unix Files

npts1 = 256;
npts2 = 256;
map = 'gray';
! touch Aimagdata.bmp
! touch Arealdata.bmp
! touch Amagnitude.bmp
! touch Aphase.bmp

%initialize and read filenames for the the data files

REAL_FILE =input('Real Channel Datafile:  ','s');
disp(' ');
IMAG_FILE = input('Imaginary Channel Datafile: ','s');
disp(' ');

%check to see if user wants to phase-modulate raw data for shifts

center_x_flag = 0; %input('Shift along x? (0 = no, 1 = yes)  ');
disp(' ');
center_y_flag = 0; %input('Shift along y? (0 = no, 1 = yes)  ');
disp(' ');

%read in  the real data

fid_r    = fopen(REAL_FILE, 'r', 'n');
header_r = fread(fid_r, [1,128], 'float32');
real_data = fread(fid_r, [256,256], 'float32');
real_data = rot90(real_data,2);

%read in  the imaginary data
fid_i    = fopen(IMAG_FILE,'r','n');
header_i = fread(fid_i, [1,128], 'float32');
imag_data = fread(fid_i, [256,256], 'float32');
imag_data = rot90(imag_data,2);


%%%%%%%%%%%%%%  Display Real & Imaginary Data
figure(1)
colormap(map);
image(real_data)
title('Real Data')
imwrite(real_data,gray,'Arealdata.bmp','bmp');
figure(2)
colormap(map);
image(imag_data)
title('Imag Data')
map = colormap;
imwrite(imag_data,gray,'Aimagdata.bmp','bmp');
%%%%%%%%%%%%%%%

%if necessary, modulate k-space data to shift center-of-FOV to center-of-image

if center_x_flag==1
  for i=1:2:255
    real_data(:,i) = -real_data(:,i);
    imag_data(:,i) = -imag_data(:,i);
  end
end

if center_y_flag==1
  for i=1:2:255
    real_data(i,:) = -real_data(i,:);
    imag_data(i,:) = -imag_data(i,:);
  end
end


%process imaginary and real data

imag_data = imag_data.*j;
combined_data = imag_data + real_data;
mri_image = fft2(combined_data);

% Adjust Image position
%
% A  B   C        D   E   F
%            ->   
% D  E   F        A   B   C
%

for i = 1:1:128;
  for j = 1:1:256;
  temp(i,j) = mri_image(i+128,j);
  end
end
for i = 128:1:255;
  for j = 1:1:256;
  temp(i,j) = mri_image(i-127,j);
  end
end
mri_image = temp;


%display magnitude image

scale = 100/(max(max(abs(mri_image)))); %Set color scale
mri_image = mri_image.*scale;
figure(3)
colormap(map);
image(abs(mri_image))
title('Magnitude Image')
imwrite(abs(mri_image),gray,'Amagnitude.bmp','bmp');

%display phase image

figure(4)
colormap(map);
scale = 100/(max(max(angle(mri_image))));
angle_mri = (angle(mri_image)).*scale;
image(angle_mri)
title('Phase Image')
imwrite(angle_mri,gray,'Aphase.bmp','bmp');
