changing circular shift in hologram reconstruction

Question

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;

hologram = imread([file_path file_name]);
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:

At z = 0.011:

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.

Answer 1

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;

hologram = double(rgb2gray(imread('test.jpg')));
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;