[question] how to get phase correctly

[funs]

Dear All,

How is everything during COVID-19. I wish you guys be good in health.

I have a question about how to get phase distribution. I follow the example and use the following code:

# calculate amp and phase for all range
amp = np.abs(recons[:,:,:])
phase = np.angle(recons[:,:,:])
phase_unwrap = unwrap_phase(phase,wrap_around=True,seed=1)

hp.show(amp) # display result of intensity
hp.show(phase) # display result of phase
hp.show(phase_unwrap) # display result of unwrap phase

# leveling
leveling = cv2.blur(phase_unwrap,ksize=(100,100))    # leveling
phase_leveling = phase_unwrap - leveling
hp.show(phase_leveling) # display result of leveling unwrap phase

However, the phase diagram is crash. Please help to figure out how to reconstruct phase correctly.

Thank you very much.

The figure from left is: amplitude, phase, unwrap phase.

[briandleahy]

Hi @funs , can you provide a bit more information? Where are you getting recons from? Is that something you generated with holopy, or is it data you collected?

If it is something you generated with holopy, can you post a minimum working example so I can reproduce the problem on my machine?

If it is data that you collected, can you tell me a little bit about the recorded hologram?

[funs]

Dear @briandleahy ,

Thank you for your reply. The following is the code and the images. ps. image from https://unal-optodigital.github.io/DLHM/

Please help me to reconstruct the phase diagram.

# -*- coding: utf-8 -*-

# import packages
import os
import holopy as hp
import numpy as np
from holopy.core.io import get_example_data_path
from holopy.core.io import save_images
from holopy.propagation import ps_propagate
from scipy.ndimage.measurements import center_of_mass
from skimage.restoration import unwrap_phase
import cv2

# set image path

rawholoimgspath = os.getcwd()+"\\rawholoimgs"
imagepath = "\\20210423 test imgs\\diatom_hologram.jpg"
bgpath = "\\20210423 test imgs\\diatom_reference.jpg"

imagepath = rawholoimgspath + imagepath
bgpath = rawholoimgspath + bgpath

rimagepath = os.getcwd()+'\\temp'+'\\resizeholo.bmp'
rbgpath = os.getcwd()+'\\temp'+'\\resizebg.bmp'

# read image and resize
zoom = 1
image = cv2.imread(imagepath)
image = cv2.resize(image, (round(image.shape[1]*zoom),round(image.shape[0]*zoom)), interpolation=cv2.INTER_CUBIC)
cv2.imwrite(rimagepath,image)
bg = cv2.imread(bgpath)
bg = cv2.resize(bg, (round(bg.shape[1]*zoom),round(bg.shape[0]*zoom)), interpolation=cv2.INTER_CUBIC)
cv2.imwrite(rbgpath,bg)

# holopy parameters
L = 10e-3
cam_spacing = 6.1444e-6
rcam_spacing = cam_spacing/zoom
mag = 30
npix_out =  1024*zoom
zstack = np.arange(300e-6, 305e-6, 1e-6)
print("zstack number is:"+str(zstack.shape))
illum_wavelen = 405e-9
medium_index = 1.0
med_wavelen = illum_wavelen / medium_index
xmin = 400#600
xmax = 600#1000
ymin = 400#900
ymax = 600#1300

# load image
rawholo = hp.load_image(rimagepath, spacing=rcam_spacing, illum_wavelen=illum_wavelen, medium_index=medium_index,channel=0) # load hologram
bg = hp.load_image(rbgpath, spacing=rcam_spacing,channel=0) # load background image

# back ground correction
#holo = rawholo
#holo = hp.core.process.bg_correct(rawholo, bg+1, bg) # subtract background (not divide)
holo = hp.core.process.bg_correct(rawholo, bg) # divide background

beam_c = center_of_mass(bg.values.squeeze()) # get beam center
#beam_c = (1000,1300)
#beam_c = center_of_mass(holo.values.squeeze()) #get beam center by holo

#print('load image done!')
out_schema = hp.core.detector_grid(shape=npix_out, spacing=rcam_spacing/mag) # set output shape
#print('set output shape done!')

# propagete process    
recons = ps_propagate(holo, zstack, L, beam_c, out_schema) # do propagation

# calculate amp and phase for ROI between x y max and min
#hp.show(recons[:,xmin:xmax,ymin:ymax])
amp = np.abs(recons[:,xmin:xmax,ymin:ymax])
phase = np.angle(recons[:,xmin:xmax,ymin:ymax])
#phase = np.arctan2(recons[:,xmin:xmax,ymin:ymax].imag, recons[:,xmin:xmax,ymin:ymax].real)
phase_unwrap = unwrap_phase(phase,wrap_around=False,seed=None)

hp.show(amp) # display result of intensity
hp.show(phase) # display result of phase
hp.show(phase_unwrap) # display result of unwrap phase

# leveling
leveling = cv2.blur(phase_unwrap,ksize=(100,100))    # leveling
phase_leveling = phase_unwrap - leveling
hp.show(phase_leveling) # display result of leveling unwrap phase

[briandleahy]

Hi @funs

There is a lot going on in the code snippet that you included above; it is not easy for me to follow and much of it is outside holopy -- phase unwrapping from skimage, blurring from cv2, etc. Can you please include a shorter description that points to just the problems you are having with holopy?

To me, the problem seems to be a physics or imaging problem and not a holopy problem per se, but I think I can help you with it once you can give me a better idea of what is going on.

[funs]

Dear @briandleahy ,

The steps of this idea is:

1. load raw image and background
2. resize image by interpolation to increase resolution
3. calculate complex 2D-array by holopy
4. calculate amplitude by np.abs()
5. calculate phase by np.angle() 5-1. unwrap phase by skimage 5-2. subtract background which was process by blurring from opencv (cv2)

My problem is : the result of step 5-1, the unwraped phase image is crashed.

There is a reference paper denotes the phase retrieval algorithm. Should I implement this algorithm?

Anand, V., Katkus, T., Linklater, D. P., Ivanova, E. P., & Juodkazis, S. (2020). Lensless three-dimensional quantitative phase imaging using phase retrieval algorithm. Journal of Imaging, 6(9), 99.