# changing circular shift in hologram reconstruction

I'm trying to write a very simple code for the reconstruction of an in-line hologram of particles in a micro channel (digital holographic microscopy), but I'm getting some results that are hard for me to understand.

clear;

lambda = 632.8e-9/1.333; %wavelength
pix_size = 7.4e-6; %camera's pixel size
dx = pix_size;
dy = pix_size;

N = size(hologram);

tukey_win = tukeywin(N(1),0.075); %tukey window
% multiply image (single precision) with the tukey window to darken edges
trimmed_holo = single(hologram) .* (tukey_win * tukey_win');

prop_distance = 0.01:0.005:0.2; % propagation distances vector
fft_hologram = fft2(trimmed_holo); %Fourier transform of hologram
fft_shifted_hologram = fftshift(fft_hologram); %Fourier shift of hologram

for k = 1:length(prop_distance) % do multiple propagation distances
% obtain transfer function
G = trans_func(lambda, N(1), dx, dy, prop_distance(k));

%multiply FT of impulse response times FT of hologram
psi = G .* fft_shifted_hologram;

% get slice of reconstructed optical field through inverse FT
recon_image = ifftshift(ifft2(psi));

%write image to file
imwrite(uint8(abs(recon_image)), [file_path_save 'out' num2str(k,'%2.2d') '.tif']);

end

function G = trans_func(lambda, N, Delta_xi, Delta_eta, z)

M = N;
[m, n] = meshgrid(-N/2:N/2-1,-M/2:M/2-1);

first_term  = (lambda^2 * (n + (N^2 * Delta_xi^2)  ./ (2 * z * lambda)).^2) ...
./ (N^2 * Delta_xi^2);
second_term = (lambda^2 * (m + (M^2 * Delta_eta^2) ./ (2 * z * lambda)).^2) ...
./ (M^2 * Delta_eta^2);
arg_root = 1 - first_term - second_term;
arg_root(arg_root<0) = 0;

G = exp(-2*pi*1i*z/lambda * sqrt(arg_root));


The actual reconstruction appears to be correct, but I am getting a circular shift that varies with the reconstruction (propagation distance). Here's the original hologram: Here's z = 0.01, which (just) looks like I missed the ifftshift(): Here's at z = 0.0105: The circular shift keeps getting smaller, and by the last propagation distance (z = 0.2) it is practically gone: I've looked at the code for quite some time and can't figure out what's causing the circular shift, so I would greatly appreciate it if anyone could explain to me what is causing/how to correct it. Thanks.

The problem is this line:

recon_image = ifftshift(ifft2(psi));


It should be

recon_image = ifft2(ifftshift(psi));


You need to shift the 0 frequency back to the corners before doing ifft2

I find fftshift / ifftshift and creating the frequency coordinates causes so many programming errors. This is my code for creating the coordinates (with 0 frequency in the corner):

function [ux,uy,r,th] = fftMesh(sizeV)
%FFTMESH - Create DFT cartesian and polar coordinates
%
% Inputs:
%
% sizeV     [width,height] of spectrum (same as image size)
%
% Outputs:
%
% ux        x in cartesian coords
% uy        y in cartesian coords
% r         radius in polar coords
% th        angle (rad CCW) in polar coords

rows = sizeV(1);
cols = sizeV(2);

[ux, uy] = meshgrid(...
([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2))...
);
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

end


If you use this, you dont have to worry about using fftshift / ifftshift. Alternatively, just remove the ifftshift lines in the function to keep the 0 freq in the middle.

Updated code:

clear;

lambda = 632.8e-9/1.333; %wavelength
pix_size = 7.4e-6; %camera's pixel size
dx = pix_size;
dy = pix_size;

N = size(hologram);

tukey_win = tukeywin(N(1),0.075); %tukey window
% multiply image (single precision) with the tukey window to darken edges
trimmed_holo = single(hologram) .* (tukey_win * tukey_win');

prop_distance = 0.01:0.005:0.2; % propagation distances vector
fft_hologram = fft2(trimmed_holo); %Fourier transform of hologram

for k = 1:length(prop_distance) % do multiple propagation distances
% obtain transfer function
G = trans_func(lambda, N(1), dx, dy, prop_distance(k));

%multiply FT of impulse response times FT of hologram
psi = G .* fft_hologram;

% get slice of reconstructed optical field through inverse FT
recon_image = ifft2(psi);

%write image to file
imagesc(real(recon_image)); pause;
imagesc(imag(recon_image)); pause;
imagesc(abs(recon_image))

end

function G = trans_func(lambda, N, Delta_xi, Delta_eta, z)

M = N;
[m, n] = meshgrid(...
([1:M]-(fix(M/2)+1))/(M-mod(M,2)), ...
([1:N]-(fix(N/2)+1))/(N-mod(N,2))...
);
m = ifftshift(m);   % Quadrant shift to put 0 frequency at the corners
n = ifftshift(n);

first_term  = (lambda^2 * (n + (N^2 * Delta_xi^2)  ./ (2 * z * lambda)).^2) ...
./ (N^2 * Delta_xi^2);
second_term = (lambda^2 * (m + (M^2 * Delta_eta^2) ./ (2 * z * lambda)).^2) ...
./ (M^2 * Delta_eta^2);
arg_root = 1 - first_term - second_term;
arg_root(arg_root<0) = 0;

• Thanks for pointing this out. It partially solved the problem, only in the first propagation distance (z = 0.01). The following reconstructed images still have this circular shift that changes depending on propagation distance z. Here are the second, third and fourth images in the propagation sequence. May 2, 2018 at 15:11
• I changed the coordinate generation to my method and removed the fftshifts and I dont get the shift. However there is another problem - G is not conjugate symmetric, so there is an imaginary part to recon_image. This doesn't sound right, why do you multiply by exp(-2*pi*1i*z/lambda)? May 2, 2018 at 15:53
• I tried your method and the circular shift does go away, but the resulting image (reconstruction) is not right. May 2, 2018 at 18:21
• @supvato Your transfer function creates an image with imaginary components. Do you have some reference material where you got the formula from? May 2, 2018 at 18:25
• There are many references for holographic reconstruction. I am using the Handbook of Holographic Interferometry by Kreis (2005). Here are some of the different alternatives for the reconstruction; I am using the third option, Transfer function G. The reference wave r = 1 for my data. Therefore, all I have to do is get the FFT of the hologram h and multiply it by G and then get the iFFT. May 5, 2018 at 2:05