search

kaanaksit / odak

Scientific computing library for optics, computer graphics and visual perception.
https://kaanaksit.com/odak
Mozilla Public License 2.0
177 stars 52 forks source link

Band-limited angular spectrum method #10

Closed askaradeniz closed 3 years ago

askaradeniz commented 4 years ago

It would be good if we had an alternative beam propagation method in wave and learn submodules.

Band-Limited Angular Spectrum Method is described in this paper and also discussed in here. It seems to work well with practical applications.

kaanaksit commented 4 years ago

After reading this paper that you have shared earlier. I believe IR Fresnel, which stands for impulse response Fresnel propagation is the same as Angular spectrum method. Commit a34bd62d9100bb90f248e098e56d1bd45265d2e7 enables the key word Angular Spectrum for propagation type string.

Bandlimited angular spectrum method takes the fourier transform of an impulse response in the Angular Spectrum method, and multiplies it with a mask that low passes frequencies. The cut-off for that low pass operation is a function of wavelength and hologram size and calculated using Equation 13 in their paper. I have also added an implementation of Bandlimited angular spectrum method for wave. With da947041fa7509924b858bb5146867b401d4987e, you should now be able to use the key word Bandlimited Angular Spectrumwhile using odak.wave.propagate_beam() definition.

@askaradeniz may I kindly ask you to verify if this is providing a better result for your example? If so, we should together plan on migrating that to learn submodule. But first we have to verify. As a reminder, in the test folder, there is a test script called test_beam_propagation.py, I left a commented line there for you to switch to this new beam propagation method in case you want to use that test script.

Also thanks to @praneethc as he verified my understanding of the paper.

askaradeniz commented 4 years ago

@kunguz Thank you very much for your time. It works nice when I try the example test case in test_beam_propagation.py.

However, I am having some trouble with obtaining results for images. These are the results: IR Fresnel fresnel Bandlimited IR Fresnel band_limited

This is the code I used for the image example:

from odak import np
from odak.wave import wavenumber, calculate_amplitude
from odak.tools import load_image, resize_image
import cv2
from skimage.measure import compare_psnr

from odak.wave import propagate_beam

## Load image.
im = load_image("data/0103.png")
if np.__name__ == 'numpy':
    im = np.float64(im/255.)
else:
    im = im/255.

im = im[:,:, 0]
im = im[:512,:512]

## Construct field.
field = np.zeros((im.shape[0],im.shape[1]),dtype=np.complex64)
field[:,:] = im[:,:]
#random_phase = np.pi*np.random.random(field.shape)
random_phase = 0
field = field*np.cos(random_phase) + 1j*field*np.sin(random_phase)

## Prop. variables.
#propagation_type = "Bandlimited IR Fresnel"
propagation_type = "IR Fresnel"
wavelength                 = 0.5*pow(10,-6)
pixeltom                   = 6*pow(10,-6)
distance                   = 0.2
k                          = wavenumber(wavelength)

## propagate to slm plane.
hologram = propagate_beam(field,k,distance,pixeltom,wavelength,propagation_type)

## propagate to image plane.
reconstruction = propagate_beam(hologram,k,-distance,pixeltom,wavelength,propagation_type)
reconstruction_amplitude = calculate_amplitude(reconstruction)

ratio = (np.sum(reconstruction_amplitude * im, (-2, -1))
            / np.sum(reconstruction_amplitude * reconstruction_amplitude, (-2, -1)))
reconstruction_amplitude *= ratio
reconstruction_amplitude = np.clip(reconstruction_amplitude, 0, 1)
print(compare_psnr(im, reconstruction_amplitude))
reconstruction_amplitude = np.uint8(reconstruction_amplitude*255)

cv2.namedWindow("recon_amp", cv2.WINDOW_NORMAL)
cv2.imshow("recon_amp", reconstruction_amplitude)
cv2.namedWindow("target_amp", cv2.WINDOW_NORMAL)
cv2.imshow("target_amp", im)
cv2.waitKey(0)
kaanaksit commented 4 years ago

Can you please rerun your code using 1021a28cc869822e5b6d4e3702fb3d609f586117?

Hopefully you will get this with bandlimiting: resim and this is without bandlimiting: resim

askaradeniz commented 4 years ago

It works great!

Using the above code with distance=0.1, I got the following results. IR Fresnel no_bandlimit no_bandlimit2

Bandlimited IR Fresnel bandlimit bandlimit_0 2

kaanaksit commented 4 years ago

This is encouraging. I am trying to invite one person to translate this piece of code from odak.wave to odak.learn. Before we proceed any further, let's wait for that person to reply. She may be able to able to assist us (@praneethc knows about this person). But if you think we can port it super quick in 5 mins, let's do that --which I believe we can using the same amount of time that I dedicated to write this reply :) --. This heavily depends on your need at this point, so it is your call. We can either wait or you (@askaradeniz) can port.

askaradeniz commented 4 years ago

Okay, I think we can wait for the reply while I work on the other issues.

kaanaksit commented 4 years ago

On a second note, Angular spectrum and Impulse Response methods for beam propagation are different. Previously, I stated that they are the same, this is wrong. I checked the kernel in Computational Fourier Optics book by David Voelz, Impulse response slightly differs. Later, I saw that this paper also identifies them as two separate papers.

This commit bcfce8c2667a7553ae404615661df354acba67e1 introduces two separate keys for propagation_type for the odak.wave.propagate_beam definition. These are Angular Spectrum and IR Fresnel.

rongduo commented 4 years ago

I wrote the pytorch version of the band-limited angular spectrum function and created a pull request to #14 .....

kaanaksit commented 4 years ago

Thanks @rongduo , I will now go ahead and merge your version

kaanaksit commented 4 years ago

With b674184480a2a95cbb6925a3228f774daef182f0, I confirm that the test/test_learn_propagate_beam.py returns a meaningful outcome, thanks to @rongduo and @askaradeniz .

Input field: resim Propagated field: resim Backpropagated field: resim

askaradeniz commented 3 years ago

It looks like we get different results from torch and numpy versions of Bandlimited Angular Spectrum method (https://github.com/kunguz/odak/commit/6a7470798b401d4934522ec7ebddc78d2c3b2e38).

This is an example output of compare function written in test_learn_beam_propagation.py:

Mismatched elements: 250000 / 250000 (100%)
Max absolute difference: 394.44279377
Max relative difference: 474.19269395
 x: array([[-61.937 +89.58j ,  32.14 -102.063j,  36.588 +44.577j, ...,
         78.629 +52.555j, -75.779 -20.299j,  63.486 -35.612j],
       [ 54.121 +19.182j, -27.26  +16.655j, -59.53   -9.736j, ...,...
 y: array([[  -6.57 +296.346j, -105.022-301.575j,  261.531+215.636j, ...,
         -11.587  -8.539j,  -43.761 +84.96j ,   55.72 -199.442j],
       [ -24.799-185.679j,  118.832+198.915j, -265.152-154.982j, ...,...
kaanaksit commented 3 years ago

Looks like we need help here. @askaradeniz are you available to revise the code and identify the origin of this difference?

askaradeniz commented 3 years ago

If I'm not missing something, the only difference I noticed are these two lines: https://github.com/kunguz/odak/blob/6a7470798b401d4934522ec7ebddc78d2c3b2e38/odak/wave/classical.py#L48 https://github.com/kunguz/odak/blob/6a7470798b401d4934522ec7ebddc78d2c3b2e38/odak/learn/wave/classical.py#L44 I tried to change this line but the results still do not match.

kaanaksit commented 3 years ago

That is a good catch. Let me rephrase this once again as I created new definitions to avoid a stand-alone beam_propagation code with tons of code.

So here is the latest and greatest for h in the band-limited angular spectrum method:

https://github.com/kunguz/odak/blob/4304e956809188694ce74d00d3f54b42e3e8cf24/odak/wave/classical.py#L109

We need two things here, one is to migrate to the same h in the odak.learn.wave submodule. Plus, we need to move band-limited-angular-spectrum method to a definition just as in odak.wave : https://github.com/kunguz/odak/blob/4304e956809188694ce74d00d3f54b42e3e8cf24/odak/wave/classical.py#L82-L119 and we need to reflect that change to the if statement in odak.learn.wave by following odak.wave again: https://github.com/kunguz/odak/blob/4304e956809188694ce74d00d3f54b42e3e8cf24/odak/wave/classical.py#L35-L36

@rongduo do you have availability for us to help fix this so that we get the same result from both in odak.wave and odak.learn.wave? Thanks!

kaanaksit commented 3 years ago

I have recently merged @rongduo 's update by rebasing. Torch implementation still does not match with numpy/cupy implementation. But the absolute difference decreased from 300s to 15. The bottleneck here is slight differences in odak.learn.wave.set_amplitude and torch's ffts. @askaradeniz I will close this issue, but before I do that, please verify if this produces decent holograms at your end.

Current:

Mismatched elements: 250000 / 250000 (100%)
Max absolute difference: 14.61796575
Max relative difference: 104.98545789
 x: array([[-233.767+825.579j,  403.702-631.955j,  -13.526-349.845j, ...,
        -316.94 +131.304j, -120.413 +29.046j,  431.921-208.002j],
       [ 595.648+113.056j,   88.465-388.643j,  146.171-622.236j, ...,...
 y: array([[-232.49 +8.239e+02j,  400.054-6.306e+02j,  -10.546-3.497e+02j,
        ..., -318.999+1.294e+02j, -114.324+2.835e+01j,
         429.296-2.082e+02j],...

Previous

Mismatched elements: 250000 / 250000 (100%)
Max absolute difference: 394.44279377
Max relative difference: 474.19269395
 x: array([[-61.937 +89.58j ,  32.14 -102.063j,  36.588 +44.577j, ...,
         78.629 +52.555j, -75.779 -20.299j,  63.486 -35.612j],
       [ 54.121 +19.182j, -27.26  +16.655j, -59.53   -9.736j, ...,...
 y: array([[  -6.57 +296.346j, -105.022-301.575j,  261.531+215.636j, ...,
         -11.587  -8.539j,  -43.761 +84.96j ,   55.72 -199.442j],
       [ -24.799-185.679j,  118.832+198.915j, -265.152-154.982j, ...,...
askaradeniz commented 3 years ago

It works well. Thanks @kunguz and @rongduo .