clear all
close all

printResults = 1;

% load a slice of the data

dataA = double( nrrdLoad('caseD00920-dwi-slice40.nhdr') );

% just grab one of the gradient images

data = squeeze( dataA(:,:,1,20 ) );

% plot the original image

figure, imagesc(data), colormap(gray), title('Original image');

% Fourier transform it

f = fft2( data );
fs = fftshift( f );

% this will be complex, so let's look at the magnitude images

as = abs(f);
afs = abs(fs);

% plot the results

figure, imagesc( as ), title('|F(data)|'), colorbar, grid on;
figure, imagesc( afs ), title('|F(data)| centered'), colorbar, grid on;

figure, imagesc( log(as) ), title('log |F(data)|'), colorbar, grid on;
figure, imagesc( log(afs) ), title('log |F(data)| centered'), colorbar, grid on;

% cut out the part that should be zero in a perfect world

fc = fs;
fc(1:56,:) = 0;
fc(end-55:end,:) = 0;
fc(:,1:56) = 0;
fc(:,end-55:end) = 0;

figure, imagesc( log( abs(fc) ) ), title('log |F(data)| centered and cut'), colorbar, grid on

% reverse everything

fci = ifftshift( fc );

tsti = ifft2( fci );
atsti = abs(tsti);

% round it, so we can deal with integer values again
% since data-size seems to be our main enemy

ratsti = round( atsti );

di = atsti-data;

dri = ratsti-data;

% plot the reconstructed image

figure, imagesc( ratsti ), colormap( gray ), title( 'Reconstructed image' )

% and the difference image

figure, imagesc( dri ), colormap( gray ), title( 'Reconstruction difference' );

if ( printResults )
  
  figure( 1 ); print -dpng original.png
  figure( 2 ); print -dpng fouriermagnitude.png
  figure( 3 ); print -dpng fouriermagnitude-centered.png
  figure( 4 ); print -dpng logfouriermagnitude.png
  figure( 5 ); print -dpng logfouriermagnitude-centered.png
  figure( 6 ); print -dpng logfouriermagnitude-centered-and-cut.png
  figure( 7 ); print -dpng reconstructed.png
  figure( 8 ); print -dpng reconstruction-diffference.png
  
end