Does it work with phase-only hologram

[yellow-force]

Hi Fred (and others),

How can I use my own phase/amplitude pattern for calculating the Fresnel and Fraunhofer diffraction?

In my project I have a phase-only hologram (a 2D hologram with amplitude = 1 and phase distribution calculated from an iterative Fourier transform algorithm, called Gerchberg Saxton algorithm). The hologram can be expressed as H = 1 x exp{j*h(x,y)}, where h(x,y) represents the phase angle and is already known. The intensity of the Fourier transform of H is my target image.

Could you please suggest how I can use the known h(x,y) to reproduce the target image? How can I use h(x,y) as the "Field"?

Many thanks!

[yellow-force]

To make the testing even simpler. Forget about the phase-only hologram. Simply take an image that contains a letter "A". Assume that this image is called I. Applying the inverse FFT: H(x,y) = ifft(I).

If I then use H(x,y) as my input aperture, illuminated by a plane wave. In the far-field diffraction I should get the original image "A". Alternatively, if I add a lens and should be able to reproduce the image "A" at the focal plane.

I tried different ways but just can not reproduce the image. Please help! Thanks!!!

[ldoyle]

Hi yellow-force, this should be possible with LightPipes. For the forward direction, start with an empty field (LP.Begin) and overlay your phase with MultPhase or SubPhase. Propagate this to the far field by a large distance or into the focal plane of a lens. See the manual and examples to get startet. [1]

An example of a Gerchberg-Saxton iteration with known intensity and unkown phase can be found at [2] or in the Examples directory of the git repository. Maybe the Fourier-Optics example [3] can be of inspiration for you, too.

[1] https://opticspy.github.io/lightpipes/description_of_lightpipes_for_python.html#spherical-coordinates [2] https://opticspy.github.io/lightpipes/examples_of_lightpipes_for_python.html#phase-recovery [3] https://opticspy.github.io/lightpipes/examples_of_lightpipes_for_python.html#fourier-optics

[yellow-force]

Hi ldoyle,

Thank you for your reply. Unfortunately, I still can not reproduce the image. Please see below my short script.

If the image is 512 x 512, should I use "Begin" to create a field of exactly 512 by 512, or should I make it bigger (zero-padding)?

import numpy as np import matplotlib.pyplot as plt import matplotlib.image as mpimg from scipy.fftpack import fftn, ifftn, fftshift from LightPipes import *

image = mpimg.imread("an example image...") ### image size 512 x 512 imgray = np.sum(image,axis=2) ### get grayscale

h = ifftn(imgray)

Field = Begin(16mm, 0.52um, 512) Field = SubIntensity(np.abs(h),Field) Field = SubPhase(np.angle(h),Field) Field = Lens(0.75m,0,0,Field) Field = Fresnel(0.75m, Field) Im_Field = Intensity(0,Field) plt.figure() plt.imshow(Im_Field, cmap='gray') plt.colorbar()

[ldoyle]

I gave it a quick try, you were very close. The output image looked very dark and wrong, but viewing in log scale revealed the problem: You need to use fftshift to "center" the image as seen by the FFT. Use the modified script to see the difference for yourself!

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from scipy.fftpack import fftn, ifftn, fftshift, ifftshift
from LightPipes import *

image = mpimg.imread('letterA_512.png') ### image size 512 x 512
imgray = np.sum(image,axis=2) ### get grayscale

h = ifftn(imgray)  #wrong, image barely visible in corners
# h = fftshift(ifftn(ifftshift(imgray))) #correct, letter visible

size = 16*mm
lam = 520*nm
N=512

f=0.75*m

Field = Begin(size, lam, N)
Field = SubIntensity(np.abs(h),Field)
Field = SubPhase(np.angle(h),Field)
Field = Lens(f,0,0,Field)
Field = Fresnel(f, Field)
Im_Field = Intensity(0,Field)
# Ph = Phase(Field, unwrap=True)
plt.figure()
fig, axs = plt.subplots(1, 2)
ax1, ax2 = axs
im1 = ax1.imshow(Im_Field, cmap='gray')
ax1.set_title('lin scale')
ax2.imshow(np.log(Im_Field), cmap='gray')
ax2.set_title('log scale')
# fig.colorbar(im1, ax=ax1)

Attachments: source image, output plot for wrong version, right version.

Very useful comment on fftshift and ifftshift I recently found in [1]:

# using the convention
x = ifft(fft(x)) # x is in real space
X = fft(ifft(X)) # X is in Fourier space
# and both 0-centered in middle of array:
X = fftshift(fft(ifftshift(x)))  # correct magnitude and phase
x = fftshift(ifft(ifftshift(X)))  # correct magnitude and phase
X = fftshift(fft(x))  # correct magnitude but wrong phase !
x = fftshift(ifft(X))  # correct magnitude but wrong phase !

[1] https://github.com/numpy/numpy/issues/13442

[yellow-force]

Many thanks for your reply. Much appreciated!

Why does the image look so dark, or why is the contrast so low? Compared with the original image, the reconstructed image has much poor quality and signal-noise-ratio.

I suppose if we insert a lens, then at the back focal plane we should be able to get the Fourier transform of the input which should be exactly, or at least very close to, the original image. But why is the quality degraded so badly?

[ldoyle]

This bugged me also and I played around once more. The answer is a missing square between intensity and field magnitude!

• The image you load is usually an intensity distribution. Take the sqrt() for a correct scaling. In an image with 0 and 1 only, of course this does not have an effect, but once you have grayscale you basically view the square or root of the reconstructed image which will look similar but slightly wrong.
• However, much more imporantly, the iFFT only works on the E-field, not the intensity. I'm sure one can write down the math why that makes sense (the lens and Fraunhofer far field are similar to FFT which is why the trick works, but only in E-field, not E^2). So interpreting h as the E-field and then squaring this in SubIntensity yields almost perfect results.

Check this code below. As a bonus I included a CircAperture to create a lowpass filter, making the A smoother when really constraining the input plane.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from scipy.fftpack import fftn, ifftn, fftshift, ifftshift
from LightPipes import *

image = mpimg.imread('letterA_512.png') ### image size 512 x 512
imgI = np.sum(image,axis=2) ### get grayscale
imgE = np.sqrt(imgI) #need the E-field before FFT
# of course with 1 and 0 this does not make a difference but for 
# completeness make this correct

h = fftshift(ifftn(ifftshift(imgE)))
hInt = np.abs(h)**2 #the sqrt above will just give wrong scaling, but
# the square here is essential for FFT to be defined in correct quantity
hPhase = np.angle(h) # probably because otherwise there is a discrepancy
# between the magnitude abs() and the angle() here...

size = 16*mm
lam = 520*nm
N=512

f=0.75*m

Field = Begin(size, lam, N)
Field = SubIntensity(hInt,Field)
Field = SubPhase(hPhase,Field)
#Field = CircAperture(size/20, 0, 0, Field) #uncomment to enable lowpass filter
I0 = Intensity(0,Field)
Field = Lens(f,0,0,Field)
Field = Forvard(f, Field)
Im_Field = Intensity(0,Field)
Ph = Phase(Field, unwrap=True)
plt.figure()
fig, axs = plt.subplots(1, 2)
ax1, ax2 = axs
im1 = ax1.imshow(np.log(I0), cmap='jet') #switch to log to see filter better
# im1 = ax1.imshow((I0), cmap='gray') #switch to lin to see most information
                                      # is actually in center of input plane
ax1.set_title('I_0 at lens plane')
ax2.imshow((Im_Field), cmap='gray')
ax2.set_title('Ifield')

I was wondering how the code should "know" the difference whether h was given by intensity or E-field magnitude, but I think the key is that SubIntensity() will take the square root to convert intensity to E-field, while SubPhase will use the angle without taking the route, since it is not conventional (or even sensibly defineable) to have the phase angle of the intensity. So if you do not use the same convention, the angles and magnitudes will not match.

For completeness here are the outputs with and without spatial filter.