Fourier Transform from a computer generated hologram / Phase Mask

I am trying to perform a fourier transform of a phase mask to reconstruct my original image after having performed the fourier transform and extracting its phase. I am able to get a reconstructed image but it seems that the original image is overlayed with a x and y mirror image of the original image.

Does anyone know how I can get rid of this or what I am doing wrong?

My code is as follows:

clear all; close all;

I = imread('batman.gif', 'gif');
I = double(I);
I = I./max(max(I));
avgl=mean(mean(I));
figure; imshow(mat2gray(I));
title('Original Object');
figure;
axis([0 , 101, 0, 1]);
xlabel ('Number of Iterations')
ylabel ('RMSE')
hold on
I1 = I;

for n=1:101; %Iterations to optimize the phase hologram

      H = fftshift((fft2(fftshift(I1))));
      I2 = fftshift(fft2(fftshift(exp(1j.*angle(H)))));
      avg2=mean(mean(abs(I2)));
      I2=(I2./avg2).*avgl;
      rmse=mean(mean((abs(I2)-I).^2))^0.5;
      plot(n,rmse,'o');
      pause(0.00001); %To see the error in each iteration.
      I1=I.*exp(1j*angle(I2));
  end

figure; imshow(H);
colormap gray
ti = get(gca,'TightInset')
set(gca,'Position',[ti(1) ti(2) 1-ti(3)-ti(1) 1-ti(4)-ti(2)]);
set(gca,'units','centimeters')
pos = get(gca,'Position');
ti = get(gca,'TightInset');

set(gcf, 'PaperUnits','centimeters');
set(gcf, 'PaperSize', [pos(3)+ti(1)+ti(3) pos(4)+ti(2)+ti(4)]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition',[0 0 pos(3)+ti(1)+ti(3) pos(4)+ti(2)+ti(4)]);

imwrite(H, 'SupermanPhaseMask2.gif')

hold off

I2=I2./max(max(abs(I2)));
figure; imshow(abs(I2));
title('Reconstucted Image')

Any help would be amazing. Thanks in advance.

[Patrick Bevington]

Hi, I realise you asked this question some time ago but the issue is that you are using:

imshow

If you use:

pcolor

this should fix it.

Thanks for the code you have provided, this has been very useful to me.

Answer 1

[Miloth]

Hi, you used the fast Fourier transform fft2 in the first line inside the for cicle. You should use the inverse transform ifft2 instead:

H = fftshift((ifft2(fftshift(I1))));

Source: Gerchberg–Saxton algorithm, Wikipedia