manoharan-lab / holopy

Hologram processing and light scattering in python
GNU General Public License v3.0
131 stars 50 forks source link

[question] Question about mag #397

Closed funs closed 7 months ago

funs commented 3 years ago

Hi all,

In the example, the mag means magnification. out_schema = hp.core.detector_grid(shape=npix_out, spacing=cam_spacing/mag)

I think it should be (the distance between point source to detector / the distance between point source to sample). In my case, the mag = 1.2. However, the MemoryError appears.

In the detector_grid() method, the definition of spacing is _spacing : int or list-like (2) If int, distance between square detector pixels. If arraylike, \ spacing\ [0] between adjacent rows and \ spacing\ [1] between adjacent columns.

Please help me to clarify this.

Thank you very much.

briandleahy commented 3 years ago

For forward modeling, the spacing parameter to detector_grid is the space between pixels in the object space of the microscope. For instance, if I have pixels that are physically separated on the CCD by 10 micrometers, and the microscope's magnification is 100x, then I would want spacing to be 0.1 micrometers (like you have in your example). I could measure this by taking an image of an object of known width W, measuring its width in pixels P, and defining spacing as W / P.

Does that answer your question?

If instead you were asking about the MemoryError, can you please post a minimum working example of code that I can run on my machine that would give me the same MemoryError?

funs commented 3 years ago

Dear @briandleahy ,

Thank you for your reply. Here is my code and experiment images(20210426 BDL live cell.zip in the following). In my experiment setup, L = 61.83 mm, Z = 22.65 mm, Mag ~ 61.83/22.65=2.73 If mag = 2.73, the MemoryError message is: ps. if mag=0.5, The results is show in the following. However, the resolution is not real. ps. memory: 8G, CPU: i5-6400

Traceback (most recent call last):

  File "D:\workingspace\LensFreeDHM\holopycms.py", line 115, in <module>
    recons = ps_propagate(holo, zstack, L, beam_c, out_schema) # do propagation

  File "C:\ProgramData\Anaconda3\lib\site-packages\holopy\propagation\point_source_propagate.py", line 63, in ps_propagate
    old_Ip, npix_plane = ps_propagate_plane(

  File "C:\ProgramData\Anaconda3\lib\site-packages\holopy\propagation\point_source_propagate.py", line 207, in ps_propagate_plane
    result = np.fromfunction(lambda i,j: Ip_calc(i,j), (npix, npix), dtype=int) #result is I'

  File "C:\ProgramData\Anaconda3\lib\site-packages\numpy\core\numeric.py", line 1809, in fromfunction
    return function(*args, **kwargs)

  File "C:\ProgramData\Anaconda3\lib\site-packages\holopy\propagation\point_source_propagate.py", line 207, in <lambda>
    result = np.fromfunction(lambda i,j: Ip_calc(i,j), (npix, npix), dtype=int) #result is I'

  File "C:\ProgramData\Anaconda3\lib\site-packages\holopy\propagation\point_source_propagate.py", line 187, in Ip_calc
    i_X = np.array( (X-X0)/Dx )

MemoryError: Unable to allocate 2.54 GiB for an array with shape (18471, 18471) and data type float64
# -*- 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 = "\\20210426 BDL live cell\\7-2.bmp"
bgpath = "\\20210426 BDL live cell\\bg1.bmp"

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)

# CMS parameters
L = 61.83e-3
cam_spacing = 2.2e-6
rcam_spacing = cam_spacing/zoom
mag = 0.5
npix_out =  1600*zoom
zstack = np.arange(34e-3, 37e-3, 0.5e-3)
print("zstack number is:"+str(zstack.shape))
illum_wavelen = 532e-9
medium_index = 1.0
med_wavelen = illum_wavelen / medium_index
xmin = 600#600
xmax = 1000#1000
ymin = 900#900
ymax = 1300#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

20210426 BDL live cell.zip

image

briandleahy commented 3 years ago

Hi @funs,

This is not a minimum working example. Can you please post a minimum working example, that contains only the lines of code necessary to reproduce the problem? Ideally this should be 2 or 3 lines, without loading or saving any files to the disk. I ask because I do not get an error when I run this code on my machine, and I am not even sure which line of code is giving you the error.

That being said, mechanically I think the MemoryError arises because as an intermediate step point_source_propagate allocates larger arrays. Your code tries to allocate a 18471 x 18471 array, which requires 2.5 GB. This isn't big, but if your computer does not have a lot of RAM, or if it has a 32 bit architecture, you may have problems.

funs commented 3 years ago

Dear @briandleahy,

Thank you for reply. I solved this problem by enlarge the RAM size from 8 Gb to 32 Gb.

I still have two questions:

  1. Is there any method to reduce the need of RAM size when I apply larger "mag" parameter ?
  2. How to speed up the process time ? Because the setup of my system is designed to mag=6, it cost lots of time to running code.

Thank you very much.

vnmanoharan commented 7 months ago

Late response but if others have this question then the following might be helpful. The number of pixels to reconstruct depends on the output spacing, which in turn depends on the magnification, so I don't think there's any way around the RAM requirements or processing time. I think the simplest way to avoid running out of RAM is to reduce the number of z-distances at which the hologram is reconstructed, or to do the reconstruction in a loop. The lines

zstack = np.arange(34e-3, 37e-3, 0.5e-3)

and

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

will return a 3D array with the reconstructed image at each z-distance, which can take up a lot of memory. One way to avoid running out of memory is to do a loop to calculate the reconstruction and save it at each z-distance. You can then process those reconstructions into a volume using a separate script or program that uses, for example, the lazy-loading functionality of xarray.