Skip to content
Narrate section

Computer-Generated Holography

In this section, we introduce Computer-Generated Holography (CGH) 12 as another emerging method to simulate light. CGH offers an upgraded but more computationally expensive way to simulating light concerning the raytracing method described in the previous section. This section dives deep into CGH and will explain how CGH differs from raytracing as we go.

What is holography?

Informative

Holography is a method in Optical sciences to represent light distribution using amplitude and phase of light. In much simpler terms, holography describes light distribution emitted from an object, scene, or illumination source over a surface by treating the light as a wave. The primary difference of holography concerning raytracing is that it accounts not only amplitude or intensity of light but also the phase of light. Unlike classical raytracing, holography also includes diffraction and interference phenomena. In raytracing, the smallest building block that defines light is a ray, whereas, in holography, the building block is a light distribution over surfaces. In other terms, while raytracing traces rays, holography deals with surface-to-surface light transfer.

Did you know this source?

There is an active repository on GitHub, where latest CGH papers relevant to display technologies are listed. Visit GitHub:bchao1/awesome-holography for more.

What is a hologram?

Informative

Hologram is either a surface or a volume that modifies the light distribution of incoming light in terms of phase and amplitude. Diffraction gratings, Holographic Optical Elements, or Metasurfaces are good examples of holograms. Within this section, we also use the term hologram as a means to describe a lightfield or a slice of a lightfield.

What is Computer-Generated Holography?

Informative

It is the computerized version (discrete sampling) of holography. In other terms, whenever you can program the phase or amplitude of light, this will get us to Computer-Generated Holography.

Where can I find an extensive summary on CGH?

You may be wondering about the greater physical details of CGH. In this case, we suggest our readers watch the video below. Please watch this video for an extensive summary on CGH 3.

Defining a slice of a lightfield

Informative · Practical

CGH deals with generating optical fields that capture light from various scenes. CGH often describes these optical fields (a.k.a. lightfields, holograms) as planes. So in CGH, light travels from plane to plane, as depicted below. Roughly, CGH deals with plane to plane interaction of light, whereas raytracing is a ray or beam oriented description of light.

Image title

A rendering showing how a slice (a.k.a. lightfield, optical field, hologram) propagates from one plane to another plane.

In other words, in CGH, you define everything as a "lightfield," including light sources, materials, and objects. Thus, we must first determine how to describe the mentioned lightfield in a computer. So that we can run CGH simulations effectively.

A lightfield is a planar slice in the context of CGH, as depicted in the above figure. This planar field is a pixelated 2D surface (could be represented as a matrix). The pixels in this 2D slice hold values for the amplitude of light, \(A\), and the phase of the light, \(\phi\) at each pixel. Whereas in classical raytracing, a ray only holds the amplitude or intensity of light. With a caveat, though, raytracing could also be made to care about the phase of light. Still, it will then arrive with all the complications of raytracing, like sampling enough rays or describing scenes accurately.

Each pixel in this planar lightfield slice encapsulates the \(A\) and \(\phi\) as \(A cos(wt + \phi)\). If you recall our description of light, we explain that light is an electromagnetic phenomenon. Here, we model the oscillating electric field of light with \(A cos(wt + \phi)\) shown in our previous light description. Note that if we stick to \(A cos(wt + \phi)\), each time two fields intersect, we have to deal with trigonometric conversion complexities like sampled in this example:

\[ A_0 cos(wt + \phi_0) + A_1 cos(wt + \phi_1), \]

Where the indices zero and one indicate the first and second fields, and we have to identify the right trigonometric conversion to deal with this sum.

Instead of complicated trigonometric conversions, what people do in CGH is to rely on complex numbers as a proxy to these trigonometric conversions. In its proxy form, a pixel value in a field is converted into \(A e^{-j \phi}\), where \(j\) represents a complex number (\(\sqrt{-1}\)). Thus, with this new proxy representation, the same intersection problem we dealt with using sophisticated trigonometry before could be turned into something as simple as \(A_0 A_1 e^{-j(\phi_0 +\phi_1)}\).

In the above summation of two fields, the resulting field follows an exact sum of the two collided fields. On the other hand, in raytracing, often, when a ray intersects with another ray, it will be left unchanged and continue its path. However, in the case of lightfields, they form a new field. This feature is called interference of light, which is not introduced in raytracing, and often raytracing omits this feature. As you can tell from also the summation, two fields could enhance the resulting field (constructive interference) by converging to a brighter intensity, or these two fields could cancel out each other (destructive interference) and lead to the absence of light --total darkness--.

There are various examples of interference in nature. For example, the blue color of a butterfly wing results from interference, as biology typically does not produce blue-colored pigments in nature. More examples of light interference from daily lives are provided in the figure below.

Image title

Two photographs showin some examples of light interference: (left) thin oil film creates rainbow interference patterns (CC BY-SA 2.5 by Wikipedia user John) and a soup bubble interference with light and creates vivid reflections (CC BY-SA 3.0 by Wikipedia user Brocken Inaglory).

We have established an easy way to describe a field with a proxy complex number form. This way, we avoided complicated trigonometric conversions. Let us look into how we use that in an actual simulation. Firstly, we can define two separate matrices to represent a field using real numbers:

import torch

amplitude = torch.tensor(100, 100, dtype = torch.float64)
phase = torch.tensor(100, 100, dtype = torch.float64)

In this above example, we define two matrices with 100 x 100 dimensions. Each matrix holds floating point numbers, and they are real numbers. To convert the amplitude and phase into a field, we must define the field as suggested in our previous description. Instead of going through the same mathematical process for every piece of our future codes, we can rely on a utility function in odak to create fields consistently and coherently across all our future developments. The utility function we will review is odak.learn.wave.generate_complex_field():

Here, we provide visual results from this piece of code as below:

Definition to generate a complex field with a given amplitude and phase.

Parameters:

  • amplitude
                Amplitude of the field.
                The expected size is [m x n] or [1 x m x n].
    
  • phase
                Phase of the field.
                The expected size is [m x n] or [1 x m x n].
    

Returns:

  • field ( ndarray ) –

    Complex field. Depending on the input, the expected size is [m x n] or [1 x m x n].

Source code in odak/learn/wave/util.py
def generate_complex_field(amplitude, phase):
    """
    Definition to generate a complex field with a given amplitude and phase.

    Parameters
    ----------
    amplitude         : torch.tensor
                        Amplitude of the field.
                        The expected size is [m x n] or [1 x m x n].
    phase             : torch.tensor
                        Phase of the field.
                        The expected size is [m x n] or [1 x m x n].

    Returns
    -------
    field             : ndarray
                        Complex field.
                        Depending on the input, the expected size is [m x n] or [1 x m x n].
    """
    field = amplitude * torch.cos(phase) + 1j * amplitude * torch.sin(phase)
    return field

Let us use this utility function to expand our previous code snippet and show how we can generate a complex field using that:

import torch
import odak # (1)

amplitude = torch.tensor(100, 100, dtype = torch.float64)
phase = torch.tensor(100, 100, dtype = torch.float64)
field = odak.learn.wave.generate_complex_field(amplitude, phase) # (2)
  1. Adding odak to our imports.
  2. Generating a field using odak.learn.wave.generate_complex_field.

Propagating a field in free space

Informative · Practical

The next question we have to ask is related to the field we generated in our previous example. In raytracing, we propagate rays in space, whereas in CGH, we propagate a field described over a surface onto another target surface. So we need a transfer function that projects our field on another target surface. That is the point where free space beam propagation comes into play. As the name implies, free space beam propagation deals with propagating light in free space from one surface to another. This entire process of propagation is also referred to as light transport in the domains of Computer Graphics. In the rest of this section, we will explore means to simulate beam propagation on a computer.

A good news for Matlab fans!

We will indeed use odak to explore beam propagation. However, there is also a book in the literature, [Numerical simulation of optical wave propagation: With examples in MATLAB by Jason D. Schmidt](https://www.spiedigitallibrary.org/ebooks/PM/Numerical-Simulation-of-Optical-Wave-Propagation-with-Examples-in-MATLAB/eISBN-9780819483270/10.1117/3.866274?SSO=1)4, that provides a crash course on beam propagation using MATLAB.

As we revisit the field we generated in the previous subsection, we remember that our field is a pixelated 2D surface. Each pixel in our fields, either a hologram or image plane, typically has a small size of a few micrometers (e.g., \(8 \mu m\)). How light travels from each one of these pixels on one surface to pixels on another is conceptually depicted as a figure at the beginning of this section (green wolf image with two planes). We will name that figure's first plane on the left as the hologram plane and the second as the image plane. In a nutshell, the contribution of a pixel on a hologram plane could be calculated by drawing rays to every pixel on the image plane. We draw rays from a point to a plane because in wave theory --what CGH follows--, light can diffract (a small aperture creating spherical waves as Huygens suggested). Each ray will have a certain distance, thus causing various delays in phase \(\phi\). As long as the distance between planes is large enough, each ray will maintain an electric field that is in the same direction as the others (same polarization), thus able to interfere with other rays emerging from other pixels in a hologram plane. This simplified description oversimplifies solving the Maxwell equations in electromagnetics.

A simplified result of solving Maxwell's equation is commonly described using Rayleigh-Sommerfeld diffraction integrals. For more on Rayleigh-Sommerfeld, consult Heurtley, J. C. (1973). Scalar Rayleigh–Sommerfeld and Kirchhoff diffraction integrals: a comparison of exact evaluations for axial points. JOSA, 63(8), 1003-1008. 5. The first solution of the Rayleigh-Sommerfeld integral, also known as the Huygens-Fresnel principle, is expressed as follows:

\[ u(x,y)=\frac{1}{j\lambda} \int\!\!\!\!\int u_0(x,y)\frac{e^{jkr}}{r}cos(\theta)dxdy, \]

where the field at a target image plane, \(u(x,y)\), is calculated by integrating over every point of the hologram's area, \(u_0(x,y)\). Note that, for the above equation, \(r\) represents the optical path between a selected point over a hologram and a selected point in the image plane, theta represents the angle between these two points, k represents the wavenumber (\(\frac{2\pi}{\lambda}\)) and \(\lambda\) represents the wavelength of light. In this described light transport model, optical fields, \(u_0(x,y)\) and \(u(x,y)\), are represented with a complex value,

\[ u_0(x,y)=A(x,y)e^{j\phi(x,y)}, \]

where \(A\) represents the spatial distribution of amplitude and \(\phi\) represents the spatial distribution of phase across a hologram plane. The described holographic light transport model is often simplified into a single convolution with a fixed spatially invariant complex kernel, \(h(x,y)\) 6.

\[ u(x,y)=u_0(x,y) * h(x,y) =\mathcal{F}^{-1}(\mathcal{F}(u_0(x,y)) \mathcal{F}(h(x,y))). \]

There are multiple variants of this simplified approach:

In many cases, people choose to use the most common form of \(h(x, y)\) described as

\[ h(x,y)=\frac{e^{jkz}}{j\lambda z} e^{\frac{jk}{2z} (x^2+y^2)}, \]

where z represents the distance between a hologram plane and a target image plane. Before, we introduce you how to use existing beam propagation in our library, let us dive deep in compiling a beam propagation code following the Rayleigh-Sommerfeld integral, also known as the Huygens-Fresnel principle. In the rest of this script, I will walk you through the below code:

import sys
import odak # (1)
import torch
from tqdm import tqdm


def main(): # (2)
    length = [7e-6, 7e-6] # (3)
    for fresnel_id, fresnel_number in enumerate(range(99)): # (4)
        fresnel_number += 1
        propagate(
                  fresnel_number = fresnel_number,
                  length = [length[0] + 1. / fresnel_number * 8e-6, length[1] + 1. / fresnel_number * 8e-6]
                 )


def propagate(
              wavelength = 532e-9, # (6)
              pixel_pitch = 3.74e-6, # (7)
              length = [15e-6, 15e-6],
              image_samples = [2, 2], # Replace it with 1000 by 1000 (8)
              aperture_samples = [2, 2], # Replace it with 1000 by 1000 (9)
              device = torch.device('cpu'),
              output_directory = 'test_output', 
              fresnel_number = 4,
              save_flag = False
             ): # (5)
    distance = pixel_pitch ** 2 / wavelength / fresnel_number
    distance = torch.as_tensor(distance, device = device)
    k = odak.learn.wave.wavenumber(wavelength)
    x = torch.linspace(- length[0] / 2, length[0] / 2, image_samples[0], device = device)
    y = torch.linspace(- length[1] / 2, length[1] / 2, image_samples[1], device = device)
    X, Y = torch.meshgrid(x, y, indexing = 'ij') # (10)
    wxs = torch.linspace(- pixel_pitch / 2., pixel_pitch / 2., aperture_samples[0], device = device)
    wys = torch.linspace(- pixel_pitch / 2., pixel_pitch / 2., aperture_samples[1], device = device) # (11)
    h  = torch.zeros(image_samples[0], image_samples[1], dtype = torch.complex64, device = device)
    for wx in tqdm(wxs):
        for wy in wys:
            h += huygens_fresnel_principle(wx, wy, X, Y, distance, k, wavelength) # (12)
    h = h * pixel_pitch ** 2 / aperture_samples[0] / aperture_samples[1] # (13) 

    if save_flag:
        save_results(h, output_directory, fresnel_number, length, pixel_pitch, distance, image_samples, device) # (14)
    return True


def huygens_fresnel_principle(x, y, X, Y, z, k, wavelength): # (12)
    r = torch.sqrt((X - x) ** 2 + (Y - y) ** 2 + z ** 2)
    h = torch.exp(1j * k * r) * z / r ** 2 * (1. / (2 * odak.pi * r) + 1. / (1j * wavelength))
    return h


def save_results(h, output_directory, fresnel_number, length, pixel_pitch, distance, image_samples, device):
    from matplotlib import pyplot as plt
    odak.tools.check_directory(output_directory)
    output_intensity = odak.learn.wave.calculate_amplitude(h) ** 2
    odak.learn.tools.save_image(
                                '{}/diffraction_output_intensity_fresnel_number_{:02d}.png'.format(output_directory, int(fresnel_number)),
                                output_intensity,
                                cmin = 0.,
                                cmax = output_intensity.max()
                               )
    cross_section_1d = output_intensity[output_intensity.shape[0] // 2]
    lengths = torch.linspace(- length[0] * 10 ** 6 / 2., length[0] * 10 ** 6 / 2., image_samples[0], device = device)
    plt.figure()
    plt.plot(lengths.detach().cpu().numpy(), cross_section_1d.detach().cpu().numpy())
    plt.xlabel('length (um)')
    plt.figtext(
                0.14,
                0.9, 
                r'Fresnel Number: {:02d}, Pixel pitch: {:.2f} um, Distance: {:.2f} um'.format(fresnel_number, pixel_pitch * 10 ** 6, distance * 10 ** 6),
                fontsize = 11
               )
    plt.savefig('{}/diffraction_1d_output_intensity_fresnel_number_{:02d}.png'.format(output_directory, int(fresnel_number)))
    plt.cla()
    plt.clf()
    plt.close()


if __name__ == '__main__':
    sys.exit(main())
  1. Importing relevant libraries
  2. This is our main routine.
  3. Length of the final image plane along X and Y axes.
  4. Fresnel number is an arbitrary number that helps to get a sense if the optical configuration could be considered as a Fresnel (near field) or Fraunhofer regions.
  5. Propagating light with the given configuration.
  6. Wavelength of light.
  7. Square aperture length of a single pixel in the simulation. This is where light diffracts from.
  8. Number of pixels in the image plane along X and Y axes.
  9. Number of point light sources used to represent a single pixel's square aperture.
  10. Sample point locations along X and Y axes at the image plane.
  11. Sample point locations along X and Y axes at the aperture plane.
  12. For each, virtual point light source defined inside the aperture, we simulate the light as if divergind point light source.
  13. Normalize with the number of samples (trapezoid integration).
  14. Rest of this code is for logistics for saving images.

We start the implementation by importing necessary libraries such as odak or torch. The first function, def main, sets the length of our image plane, where we will observe the diffraction pattern. As we set the size of our image plane, we also set a arbitrary number called Fresnel Number,

\[ n_F = \frac{w^2}{\lambda z}, \]

where \(z\) is the propagation distance, \(w\) is the side length of an aperture diffracting light like a pixel's square aperture -- this is often the pixel pitch -- and \(\lambda\) is the wavelength of the light. This number helps us to get an idea if the set optical configuration falls under a certain regime like Fresnel or Fraunhofer. Fresnel number also provides a practical ease related to comparing solutions. Regardless of the optical configuration, a result with a specific Fresnel number will follow a similar pattern with different optical configuration. Thus, providing a way to verify your solutions. In the next step, we call the light propagation function, def propagate. In the beginning of this function, we set the optical configuration. For instance, we set pixel_pitch, this is the side length of a square aperture that the light will diffract from. Inside the def propagate function, we reset the distance such that it follows the input Fresnel Number and wavelength. We define the locations of the samples across X and Y axes that will represent points to calculate on the image plane, x and y. Than, we define the locations of the samples across X and Y axes that will represent the point light source locations inside the aperture, wxs and wys, which in this case a square aperture that represents a single pixel and its sidelength is provided by pixel_pitch. The nested for loop goes over the wxs and wys. Each time, we choose a point from the aperture, we propagate a spherical wave from that point using def huygens_fresnel_principle. Note that we accumulate the effect of each spherical wave on a variable called h. This is diffraction pattern in complex form from our square aperture, and we also normalize it using pixel_pitch and aperture_samples. Here, we provide visual results from this piece of code as below:

Image title

Saved 1D intensities on image plane representing diffraction patterns for various Fresnel numbers. These patterns are generated by using "test/test_diffraction_integral.py".

Image title

Saved 2D intensities on image plane representing diffraction patterns for various Fresnel numbers. These patterns are generated by using "test/test_diffraction_integral.py".

Note that beam propagation can also be learned for physical setups to avoid imperfections in a setup and to improve the image quality at an image plane:

The above descriptions establish a mathematical understanding of beam propagation. Let us examine the implementation of a beam propagation method called Bandlimited Angular Spectrum by reviewing these two utility functions from odak:

Helper function for odak.learn.wave.band_limited_angular_spectrum.

Parameters:

  • nu
                 Resolution at X axis in pixels.
    
  • nv
                 Resolution at Y axis in pixels.
    
  • dx
                 Pixel pitch in meters.
    
  • wavelength
                 Wavelength in meters.
    
  • distance
                 Distance in meters.
    
  • device
                 Device, for more see torch.device().
    

Returns:

  • H ( complex64 ) –

    Complex kernel in Fourier domain.

Source code in odak/learn/wave/classical.py
def get_band_limited_angular_spectrum_kernel(
                                             nu,
                                             nv,
                                             dx = 8e-6,
                                             wavelength = 515e-9,
                                             distance = 0.,
                                             device = torch.device('cpu')
                                            ):
    """
    Helper function for odak.learn.wave.band_limited_angular_spectrum.

    Parameters
    ----------
    nu                 : int
                         Resolution at X axis in pixels.
    nv                 : int
                         Resolution at Y axis in pixels.
    dx                 : float
                         Pixel pitch in meters.
    wavelength         : float
                         Wavelength in meters.
    distance           : float
                         Distance in meters.
    device             : torch.device
                         Device, for more see torch.device().


    Returns
    -------
    H                  : torch.complex64
                         Complex kernel in Fourier domain.
    """
    x = dx * float(nu)
    y = dx * float(nv)
    fx = torch.linspace(
                        -1 / (2 * dx) + 0.5 / (2 * x),
                         1 / (2 * dx) - 0.5 / (2 * x),
                         nu,
                         dtype = torch.float32,
                         device = device
                        )
    fy = torch.linspace(
                        -1 / (2 * dx) + 0.5 / (2 * y),
                        1 / (2 * dx) - 0.5 / (2 * y),
                        nv,
                        dtype = torch.float32,
                        device = device
                       )
    FY, FX = torch.meshgrid(fx, fy, indexing='ij')
    HH_exp = 2 * torch.pi * torch.sqrt(1 / wavelength ** 2 - (FX ** 2 + FY ** 2))
    distance = torch.tensor([distance], device = device)
    H_exp = torch.mul(HH_exp, distance)
    fx_max = 1 / torch.sqrt((2 * distance * (1 / x))**2 + 1) / wavelength
    fy_max = 1 / torch.sqrt((2 * distance * (1 / y))**2 + 1) / wavelength
    H_filter = ((torch.abs(FX) < fx_max) & (torch.abs(FY) < fy_max)).clone().detach()
    H = generate_complex_field(H_filter, H_exp)
    return H

A definition to calculate bandlimited angular spectrum based beam propagation. For more Matsushima, Kyoji, and Tomoyoshi Shimobaba. "Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields." Optics express 17.22 (2009): 19662-19673.

Parameters:

  • field
               A complex field.
               The expected size is [m x n].
    
  • k
               Wave number of a wave, see odak.wave.wavenumber for more.
    
  • distance
               Propagation distance.
    
  • dx
               Size of one single pixel in the field grid (in meters).
    
  • wavelength
               Wavelength of the electric field.
    
  • zero_padding
               Zero pad in Fourier domain.
    
  • aperture
               Fourier domain aperture (e.g., pinhole in a typical holographic display).
               The default is one, but an aperture could be as large as input field [m x n].
    

Returns:

  • result ( complex ) –

    Final complex field [m x n].

Source code in odak/learn/wave/classical.py
def band_limited_angular_spectrum(
                                  field,
                                  k,
                                  distance,
                                  dx,
                                  wavelength,
                                  zero_padding = False,
                                  aperture = 1.
                                 ):
    """
    A definition to calculate bandlimited angular spectrum based beam propagation. For more 
    `Matsushima, Kyoji, and Tomoyoshi Shimobaba. "Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields." Optics express 17.22 (2009): 19662-19673`.

    Parameters
    ----------
    field            : torch.complex
                       A complex field.
                       The expected size is [m x n].
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    zero_padding     : bool
                       Zero pad in Fourier domain.
    aperture         : torch.tensor
                       Fourier domain aperture (e.g., pinhole in a typical holographic display).
                       The default is one, but an aperture could be as large as input field [m x n].


    Returns
    -------
    result           : torch.complex
                       Final complex field [m x n].
    """
    H = get_propagation_kernel(
                               nu = field.shape[-2], 
                               nv = field.shape[-1], 
                               dx = dx, 
                               wavelength = wavelength, 
                               distance = distance, 
                               propagation_type = 'Bandlimited Angular Spectrum',
                               device = field.device
                              )
    result = custom(field, H, zero_padding = zero_padding, aperture = aperture)
    return result

Definitions for various beam propagation methods mostly in accordence with "Computational Fourier Optics" by David Vuelz.

Parameters:

  • field
               Complex field [m x n].
    
  • k
               Wave number of a wave, see odak.wave.wavenumber for more.
    
  • distance
               Propagation distance.
    
  • dx
               Size of one single pixel in the field grid (in meters).
    
  • wavelength
               Wavelength of the electric field.
    
  • propagation_type (str, default: 'Bandlimited Angular Spectrum' ) –
               Type of the propagation.
               The options are Impulse Response Fresnel, Transfer Function Fresnel, Angular Spectrum, Bandlimited Angular Spectrum, Fraunhofer.
    
  • kernel
               Custom complex kernel.
    
  • zero_padding
               Zero padding the input field if the first item in the list set True.
               Zero padding in the Fourier domain if the second item in the list set to True.
               Cropping the result with half resolution if the third item in the list is set to true.
               Note that in Fraunhofer propagation, setting the second item True or False will have no effect.
    
  • aperture
               Aperture at Fourier domain default:[2m x 2n], otherwise depends on `zero_padding`.
               If provided as a floating point 1, there will be no aperture in Fourier domain.
    
  • scale
               Resolution factor to scale generated kernel.
    
  • samples
               When using `Impulse Response Fresnel` propagation, these sample counts along X and Y will be used to represent a rectangular aperture. First two is for a hologram pixel and second two is for an image plane pixel.
    

Returns:

  • result ( complex ) –

    Final complex field [m x n].

Source code in odak/learn/wave/classical.py
def propagate_beam(
                   field,
                   k,
                   distance,
                   dx,
                   wavelength,
                   propagation_type='Bandlimited Angular Spectrum',
                   kernel = None,
                   zero_padding = [True, False, True],
                   aperture = 1.,
                   scale = 1,
                   samples = [20, 20, 5, 5]
                  ):
    """
    Definitions for various beam propagation methods mostly in accordence with "Computational Fourier Optics" by David Vuelz.

    Parameters
    ----------
    field            : torch.complex
                       Complex field [m x n].
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    propagation_type : str
                       Type of the propagation.
                       The options are Impulse Response Fresnel, Transfer Function Fresnel, Angular Spectrum, Bandlimited Angular Spectrum, Fraunhofer.
    kernel           : torch.complex
                       Custom complex kernel.
    zero_padding     : list
                       Zero padding the input field if the first item in the list set True.
                       Zero padding in the Fourier domain if the second item in the list set to True.
                       Cropping the result with half resolution if the third item in the list is set to true.
                       Note that in Fraunhofer propagation, setting the second item True or False will have no effect.
    aperture         : torch.tensor
                       Aperture at Fourier domain default:[2m x 2n], otherwise depends on `zero_padding`.
                       If provided as a floating point 1, there will be no aperture in Fourier domain.
    scale            : int
                       Resolution factor to scale generated kernel.
    samples          : list
                       When using `Impulse Response Fresnel` propagation, these sample counts along X and Y will be used to represent a rectangular aperture. First two is for a hologram pixel and second two is for an image plane pixel.

    Returns
    -------
    result           : torch.complex
                       Final complex field [m x n].
    """
    if zero_padding[0]:
        field = zero_pad(field)
    if propagation_type == 'Angular Spectrum':
        result = angular_spectrum(
                                  field = field,
                                  k = k,
                                  distance = distance,
                                  dx = dx,
                                  wavelength = wavelength,
                                  zero_padding = zero_padding[1],
                                  aperture = aperture
                                 )
    elif propagation_type == 'Bandlimited Angular Spectrum':
        result = band_limited_angular_spectrum(
                                               field = field,
                                               k = k,
                                               distance = distance,
                                               dx = dx,
                                               wavelength = wavelength,
                                               zero_padding = zero_padding[1],
                                               aperture = aperture
                                              )
    elif propagation_type == 'Impulse Response Fresnel':
        result = impulse_response_fresnel(
                                          field = field,
                                          k = k,
                                          distance = distance,
                                          dx = dx,
                                          wavelength = wavelength,
                                          zero_padding = zero_padding[1],
                                          aperture = aperture,
                                          scale = scale,
                                          samples = samples
                                         )
    elif propagation_type == 'Seperable Impulse Response Fresnel':
        result = seperable_impulse_response_fresnel(
                                                    field = field,
                                                    k = k,
                                                    distance = distance,
                                                    dx = dx,
                                                    wavelength = wavelength,
                                                    zero_padding = zero_padding[1],
                                                    aperture = aperture,
                                                    scale = scale,
                                                    samples = samples
                                                   )
    elif propagation_type == 'Transfer Function Fresnel':
        result = transfer_function_fresnel(
                                           field = field,
                                           k = k,
                                           distance = distance,
                                           dx = dx,
                                           wavelength = wavelength,
                                           zero_padding = zero_padding[1],
                                           aperture = aperture
                                          )
    elif propagation_type == 'custom':
        result = custom(
                        field = field,
                        kernel = kernel,
                        zero_padding = zero_padding[1],
                        aperture = aperture
                       )
    elif propagation_type == 'Fraunhofer':
        result = fraunhofer(
                            field = field,
                            k = k,
                            distance = distance,
                            dx = dx,
                            wavelength = wavelength
                           )
    elif propagation_type == 'Incoherent Angular Spectrum':
        result = incoherent_angular_spectrum(
                                             field = field,
                                             k = k,
                                             distance = distance,
                                             dx = dx,
                                             wavelength = wavelength,
                                             zero_padding = zero_padding[1],
                                             aperture = aperture
                                            )
    else:
        logging.warning('Propagation type not recognized')
        assert True == False
    if zero_padding[2]:
        result = crop_center(result)
    return result

Definition for calculating the wavenumber of a plane wave.

Parameters:

  • wavelength
           Wavelength of a wave in mm.
    

Returns:

  • k ( float ) –

    Wave number for a given wavelength.

Source code in odak/learn/wave/util.py
def wavenumber(wavelength):
    """
    Definition for calculating the wavenumber of a plane wave.

    Parameters
    ----------
    wavelength   : float
                   Wavelength of a wave in mm.

    Returns
    -------
    k            : float
                   Wave number for a given wavelength.
    """
    k = 2 * np.pi / wavelength
    return k

Let us see how we can use the given beam propagation function with an example:

import sys
import os
import odak
import numpy as np
import torch


def test(output_directory = 'test_output'):
    odak.tools.check_directory(output_directory)
    wavelength = 532e-9 # (1)
    pixel_pitch = 8e-6 # (2)
    distance = 0.5e-2 # (3)
    propagation_types = ['Angular Spectrum', 'Bandlimited Angular Spectrum', 'Transfer Function Fresnel'] # (4)
    k = odak.learn.wave.wavenumber(wavelength) # (5)


    amplitude = torch.zeros(500, 500)
    amplitude[200:300, 200:300 ] = 1. # (5)
    phase = torch.randn_like(amplitude) * 2 * odak.pi # (6)
    hologram = odak.learn.wave.generate_complex_field(amplitude, phase) # (7)

    for propagation_type in propagation_types:
        image_plane = odak.learn.wave.propagate_beam(
                                                     hologram,
                                                     k,
                                                     distance,
                                                     pixel_pitch,
                                                     wavelength,
                                                     propagation_type,
                                                     zero_padding = [True, False, True] # (8)
                                                    ) # (9)

        image_intensity = odak.learn.wave.calculate_amplitude(image_plane) ** 2 # (10)
        hologram_intensity = amplitude ** 2

        odak.learn.tools.save_image(
                                    '{}/image_intensity_{}.png'.format(output_directory, propagation_type.replace(' ', '_')), 
                                    image_intensity, 
                                    cmin = 0., 
                                    cmax = image_intensity.max()
                                ) # (11)
        odak.learn.tools.save_image(
                                    '{}/hologram_intensity_{}.png'.format(output_directory, propagation_type.replace(' ', '_')), 
                                    hologram_intensity, 
                                    cmin = 0., 
                                    cmax = 1.
                                ) # (12)
    assert True == True


if __name__ == '__main__':
    sys.exit(test()) 
  1. Setting the wavelength of light in meters. We use 532 nm (green light) in this example.
  2. Setting the physical size of a single pixel in our simulation. We use \(6 \mu m\) pixel size (width and height are both \(6 \mu m\).)
  3. Setting the distance between two planes, hologram and image plane. We set it as half a centimeterhere.
  4. We set the propagation type to Bandlimited Angular Spectrum.
  5. Here, we calculate a value named wavenumber, which we introduced while we were talking about the beam propagation functions.
  6. Here, we assume that there is a rectangular light at the center of our hologram.
  7. Here, we generate the field by combining amplitude and phase.
  8. Here, we zeropad and crop our field before and after the beam propagation to make sure that there is no aliasing in our results (see Nyquist criterion).
  9. We propagate the beam using the values and field provided.
  10. We calculate the final intensity on our image plane. Remember that human eyes can see intensity but not amplitude or phase of light. Intensity of light is a square of its amplitude.
  11. We save image plane intensity to an image file.
  12. For comparison, we also save the hologram intensity to an image file so that we can observe how our light transformed from one plane to another.

Let us also take a look at the saved images as a result of the above sample code:

Image title

Saved intensities before (left_ and after (right) beam propagation (hologram and image plane intensities). This result is generated using "test/test_learn_beam_propagation.py".
Challenge: Light transport on Arbitrary Surfaces

We know that we can propagate a hologram to any image plane at any distance. This propagation is a plane-to-plane interaction. However, there may be cases where a simulation involves finding light distribution over an arbitrary surface. Conventionally, this could be achieved by propagating the hologram to multiple different planes and picking the results from each plane on the surface of that arbitrary surface. We challenge our readers to code the mentioned baseline (multiple planes for arbitrary surfaces) and ask them to develop a beam propagation that is less computationally expensive and works for arbitrary surfaces (e.g., tilted planes or arbitrary shapes). This development could either rely on classical approaches or involve learning-based methods. The resultant method could be part of odak.learn.wave submodule as a new class odak.learn.wave.propagate_arbitrary. In addition, a unit test test/test_learn_propagate_arbitrary.py has to adopt this new class. To add these to odak, you can rely on the pull request feature on GitHub. You can also create a new engineering note for arbitrary surfaces in docs/notes/beam_propagation_arbitrary_surfaces.md.

Optimizing holograms

Informative · Practical

In the previous subsection, we propagate an input field (a.k.a. lightfield, hologram) to another plane called the image plane. We can store any scene or object as a field on such planes. Thus, we have learned that we can have a plane (hologram) to capture or display a slice of a lightfield for any given scene or object. After all this introduction, it is also safe to say, regardless of hardware, holograms are the natural way to represent three-dimensional scenes, objects, and data!

Holograms come in many forms. We can broadly classify holograms as analog and digital. Analog holograms are physically tailored structures. They are typically a result of manufacturing engineered surfaces (micron or nanoscale structures). Some examples of analog holograms include diffractive optical elements 13, holographic optical elements 14, and metasurfaces 15. Here, we show an example of an analog hologram that gives us a slice of a lightfield, and we can observe the scene this way from various perspectives:

Image title

A video showing analog hologram example from Zebra Imaging -- ZScape.

Digital holograms are the ones that are dynamic and generated using programmable versions of analog holograms. Typically, the tiniest fraction of digital holograms is a pixel that either manipulates the phase or amplitude of light. In our laboratory, we build holographic displays 1612, a programmable device to display holograms. The components used in such a display are illustrated in the following rendering and contain a Spatial Light Modulator (SLM) that could display programmable holograms. Note that the SLM in this specific hardware can only manipulate phase of an incoming light.

Image title

A rendering showing a standard holographic display hardware.

We can display holograms that generate images to fill a three-dimensional volume using the above hardware. We know that they are three-dimensional from the fact that we can focus on different parts of the images by changing the focus of our camera (closely observing the camera's location in the above figure). Let us look into a sample result to see what these three-dimensional images look like as we focus on different scene parts.

Image title

A series of photographs at various focuses capturing images from our computer-generated holograms.

Let us look into how we can optimize a hologram for our holographic display by visiting the below example:

import sys
import odak
import torch


def test(output_directory = 'test_output'):
    odak.tools.check_directory(output_directory)
    device = torch.device('cpu') # (1)
    target = odak.learn.tools.load_image('./test/data/usaf1951.png', normalizeby = 255., torch_style = True)[1] # (4)
    hologram, reconstruction = odak.learn.wave.stochastic_gradient_descent(
                                                                           target,
                                                                           wavelength = 532e-9,
                                                                           distance = 20e-2,
                                                                           pixel_pitch = 8e-6,
                                                                           propagation_type = 'Bandlimited Angular Spectrum',
                                                                           n_iteration = 50,
                                                                           learning_rate = 0.1
                                                                          ) # (2)
    odak.learn.tools.save_image(
                                '{}/phase.png'.format(output_directory), 
                                odak.learn.wave.calculate_phase(hologram) % (2 * odak.pi), 
                                cmin = 0., 
                                cmax = 2 * odak.pi
                               ) # (3)
    odak.learn.tools.save_image('{}/sgd_target.png'.format(output_directory), target, cmin = 0., cmax = 1.)
    odak.learn.tools.save_image(
                                '{}/sgd_reconstruction.png'.format(output_directory), 
                                odak.learn.wave.calculate_amplitude(reconstruction) ** 2, 
                                cmin = 0., 
                                cmax = 1.
                               )
    assert True == True


if __name__ == '__main__':
    sys.exit(test())
  1. Replace cpu with cuda if you have a NVIDIA GPU with enough memory or AMD GPU with enough memory and ROCm support.
  2. We will provide the details of this optimization function in the next part.
  3. Saving the phase-only hologram. Note that a phase-only hologram is between zero and two pi.
  4. Loading an image from a file with 1920 by 1080 resolution and using green channel.

The above sample optimization script uses a function called odak.learn.wave.stochastic_gradient_descent. This function sits at the center of this optimization, and we have to understand what it entails by closely observing its inputs, outputs, and source code. Let us review the function.

Definition to generate phase and reconstruction from target image via stochastic gradient descent.

Parameters:

  • target
                        Target field amplitude [m x n].
                        Keep the target values between zero and one.
    
  • wavelength
                        Set if the converted array requires gradient.
    
  • distance
                        Hologram plane distance wrt SLM plane.
    
  • pixel_pitch
                        SLM pixel pitch in meters.
    
  • propagation_type
                        Type of the propagation (see odak.learn.wave.propagate_beam()).
    
  • n_iteration
                        Number of iteration.
    
  • loss_function
                        If none it is set to be l2 loss.
    
  • learning_rate
                        Learning rate.
    

Returns:

  • hologram ( Tensor ) –

    Phase only hologram as torch array

  • reconstruction_intensity ( Tensor ) –

    Reconstruction as torch array

Source code in odak/learn/wave/classical.py
def stochastic_gradient_descent(
                                target,
                                wavelength,
                                distance,
                                pixel_pitch,
                                propagation_type = 'Bandlimited Angular Spectrum',
                                n_iteration = 100,
                                loss_function = None,
                                learning_rate = 0.1
                               ):
    """
    Definition to generate phase and reconstruction from target image via stochastic gradient descent.

    Parameters
    ----------
    target                    : torch.Tensor
                                Target field amplitude [m x n].
                                Keep the target values between zero and one.
    wavelength                : double
                                Set if the converted array requires gradient.
    distance                  : double
                                Hologram plane distance wrt SLM plane.
    pixel_pitch               : float
                                SLM pixel pitch in meters.
    propagation_type          : str
                                Type of the propagation (see odak.learn.wave.propagate_beam()).
    n_iteration:              : int
                                Number of iteration.
    loss_function:            : function
                                If none it is set to be l2 loss.
    learning_rate             : float
                                Learning rate.

    Returns
    -------
    hologram                  : torch.Tensor
                                Phase only hologram as torch array

    reconstruction_intensity  : torch.Tensor
                                Reconstruction as torch array

    """
    phase = torch.randn_like(target, requires_grad = True)
    k = wavenumber(wavelength)
    optimizer = torch.optim.Adam([phase], lr = learning_rate)
    if type(loss_function) == type(None):
        loss_function = torch.nn.MSELoss()
    t = tqdm(range(n_iteration), leave = False, dynamic_ncols = True)
    for i in t:
        optimizer.zero_grad()
        hologram = generate_complex_field(1., phase)
        reconstruction = propagate_beam(
                                        hologram, 
                                        k, 
                                        distance, 
                                        pixel_pitch, 
                                        wavelength, 
                                        propagation_type, 
                                        zero_padding = [True, False, True]
                                       )
        reconstruction_intensity = calculate_amplitude(reconstruction) ** 2
        loss = loss_function(reconstruction_intensity, target)
        description = "Loss:{:.4f}".format(loss.item())
        loss.backward(retain_graph = True)
        optimizer.step()
        t.set_description(description)
    logging.warning(description)
    torch.no_grad()
    hologram = generate_complex_field(1., phase)
    reconstruction = propagate_beam(
                                    hologram, 
                                    k, 
                                    distance, 
                                    pixel_pitch, 
                                    wavelength, 
                                    propagation_type, 
                                    zero_padding = [True, False, True]
                                   )
    return hologram, reconstruction

Let us also examine the optimized hologram and the image that the hologram reconstructed at the image plane.

Image title

Optimized phase-only hologram. Generated using "test/test_learn_wave_stochastic_gradient_descent.py".

Image title

Optimized phase-only hologram reconstructed at the image plane, generated using "test/test_learn_wave_stochastic_gradient_descent.py".
Challenge: Non-iterative Learned Hologram Calculation

We provided an overview of optimizing holograms using iterative methods. Iterative methods are computationally expensive and unsuitable for real-time hologram generation. We challenge our readers to derive a learned hologram generation method for multiplane images (not single-plane like in our example). This development could either rely on classical convolutional neural networks or blend with physical priors explained in this section. The resultant method could be part of odak.learn.wave submodule as a new class odak.learn.wave.learned_hologram. In addition, a unit test test/test_learn_hologram.py has to adopt this new class. To add these to odak, you can rely on the pull request feature on GitHub. You can also create a new engineering note for arbitrary surfaces in docs/notes/learned_hologram_generation.md.

Simulating a standard holographic display

Informative · Practical

We optimized holograms for a holographic display in the previous section. However, the beam propagation distance we used in our optimization example was large. If we were to run the same optimization for a shorter propagation distance, say not cms but mms, we would not get a decent solution. Because in an actual holographic display, there is an aperture that helps to filter out some of the light. The previous section contained an optical layout rendering of a holographic display, where this aperture is also depicted. As depicted in the rendering located in the previous section, this aperture is located between a two lens system, which is also known as 4F imaging system.

Did you know?

4F imaging system can take a Fourier transform of an input field by using physics but not computers. For more details, please review these course notes from MIT.

Let us review the class dedicated to accurately simulating a holographic display and its functions:

Internal function to reconstruct a given hologram.

Parameters:

  • hologram_phases
                         Hologram phases [ch x m x n].
    
  • amplitude
                         Amplitude profiles for each color primary [ch x m x n]
    
  • no_grad
                         If set True, uses torch.no_grad in reconstruction.
    
  • get_complex
                         If set True, reconstructor returns the complex field but not the intensities.
    

Returns:

  • reconstructions ( tensor ) –

    Reconstructed frames.

Source code in odak/learn/wave/propagators.py
def reconstruct(self, hologram_phases, amplitude = None, no_grad = True, get_complex = False):
    """
    Internal function to reconstruct a given hologram.


    Parameters
    ----------
    hologram_phases            : torch.tensor
                                 Hologram phases [ch x m x n].
    amplitude                  : torch.tensor
                                 Amplitude profiles for each color primary [ch x m x n]
    no_grad                    : bool
                                 If set True, uses torch.no_grad in reconstruction.
    get_complex                : bool
                                 If set True, reconstructor returns the complex field but not the intensities.

    Returns
    -------
    reconstructions            : torch.tensor
                                 Reconstructed frames.
    """
    if no_grad:
        torch.no_grad()
    if len(hologram_phases.shape) > 3:
        hologram_phases = hologram_phases.squeeze(0)
    if get_complex == True:
        reconstruction_type = torch.complex64
    else:
        reconstruction_type = torch.float32
    reconstructions = torch.zeros(
                                  self.number_of_frames,
                                  self.number_of_depth_layers,
                                  self.number_of_channels,
                                  self.resolution[0] * self.resolution_factor,
                                  self.resolution[1] * self.resolution_factor,
                                  dtype = reconstruction_type,
                                  device = self.device
                                 )
    if isinstance(amplitude, type(None)):
        amplitude = torch.zeros(
                                self.number_of_channels,
                                self.resolution[0] * self.resolution_factor,
                                self.resolution[1] * self.resolution_factor,
                                device = self.device
                               )
        amplitude[:, ::self.resolution_factor, ::self.resolution_factor] = 1.
    if self.resolution_factor != 1:
        hologram_phases_scaled = torch.zeros_like(amplitude)
        hologram_phases_scaled[
                               :,
                               ::self.resolution_factor,
                               ::self.resolution_factor
                              ] = hologram_phases
    else:
        hologram_phases_scaled = hologram_phases
    for frame_id in range(self.number_of_frames):
        for depth_id in range(self.number_of_depth_layers):
            for channel_id in range(self.number_of_channels):
                laser_power = self.get_laser_powers()[frame_id][channel_id]
                phase = hologram_phases_scaled[frame_id]
                hologram = generate_complex_field(
                                                  laser_power * amplitude[channel_id],
                                                  phase * self.phase_scale[channel_id]
                                                 )
                reconstruction_field = self.__call__(hologram, channel_id, depth_id)
                if get_complex == True:
                    result = reconstruction_field
                else:
                    result = calculate_amplitude(reconstruction_field) ** 2
                reconstructions[
                                frame_id,
                                depth_id,
                                channel_id
                               ] = result.detach().clone()
    return reconstructions

Function that represents the forward model in hologram optimization.

Parameters:

  • input_field
                  Input complex input field.
    
  • channel_id
                  Identifying the color primary to be used.
    
  • depth_id
                  Identifying the depth layer to be used.
    

Returns:

  • output_field ( tensor ) –

    Propagated output complex field.

Source code in odak/learn/wave/propagators.py
def __call__(self, input_field, channel_id, depth_id):
    """
    Function that represents the forward model in hologram optimization.

    Parameters
    ----------
    input_field         : torch.tensor
                          Input complex input field.
    channel_id          : int
                          Identifying the color primary to be used.
    depth_id            : int
                          Identifying the depth layer to be used.

    Returns
    -------
    output_field        : torch.tensor
                          Propagated output complex field.
    """
    distance = self.distances[depth_id]
    if not self.generated_kernels[depth_id, channel_id]:
        if self.propagator_type == 'forward':
            H = get_propagation_kernel(
                                       nu = self.resolution[0] * 2,
                                       nv = self.resolution[1] * 2,
                                       dx = self.pixel_pitch,
                                       wavelength = self.wavelengths[channel_id],
                                       distance = distance,
                                       device = self.device,
                                       propagation_type = self.propagation_type,
                                       samples = self.aperture_samples,
                                       scale = self.resolution_factor
                                      )
        elif self.propagator_type == 'back and forth':
            H_forward = get_propagation_kernel(
                                               nu = self.resolution[0] * 2,
                                               nv = self.resolution[1] * 2,
                                               dx = self.pixel_pitch,
                                               wavelength = self.wavelengths[channel_id],
                                               distance = self.zero_mode_distance,
                                               device = self.device,
                                               propagation_type = self.propagation_type,
                                               samples = self.aperture_samples,
                                               scale = self.resolution_factor
                                              )
            distance_back = -(self.zero_mode_distance + self.image_location_offset - distance)
            H_back = get_propagation_kernel(
                                            nu = self.resolution[0] * 2,
                                            nv = self.resolution[1] * 2,
                                            dx = self.pixel_pitch,
                                            wavelength = self.wavelengths[channel_id],
                                            distance = distance_back,
                                            device = self.device,
                                            propagation_type = self.propagation_type,
                                            samples = self.aperture_samples,
                                            scale = self.resolution_factor
                                           )
            H = H_forward * H_back
        self.kernels[depth_id, channel_id] = H
        self.generated_kernels[depth_id, channel_id] = True
    else:
        H = self.kernels[depth_id, channel_id].detach().clone()
    field_scale = input_field
    field_scale_padded = zero_pad(field_scale)
    output_field_padded = custom(field_scale_padded, H, aperture = self.aperture)
    output_field = crop_center(output_field_padded)
    return output_field

This sample unit test provides an example use case of the holographic display class.


Let us also examine how the reconstructed images look like at the image plane.

Image title

Reconstructed phase-only hologram at two image plane, generated using "test/test_learn_wave_holographic_display.py".

You may also be curious about how these holograms would look like in an actual holographic display, here we provide a sample gallery filled with photographs captured from our holographic display:

Photographs of holograms captured using the holographic display in Computational Light Laboratory

Conclusion

Informative

Holography offers new frontiers as an emerging method in simulating light for various applications, including displays and cameras. We provide a basic introduction to Computer-Generated Holography and a simple understanding of holographic methods. A motivated reader could scale up from this knowledge to advance concepts in displays, cameras, visual perception, optical computing, and many other light-based applications.

Reminder

We host a Slack group with more than 300 members. This Slack group focuses on the topics of rendering, perception, displays and cameras. The group is open to public and you can become a member by following this link. Readers can get in-touch with the wider community using this public group.


  1. Max Born and Emil Wolf. Principles of optics: electromagnetic theory of propagation, interference and diffraction of light. Elsevier, 2013. 

  2. Joseph W Goodman. Introduction to Fourier optics. Roberts and Company publishers, 2005. 

  3. Koray Kavakli, David Robert Walton, Nick Antipa, Rafał Mantiuk, Douglas Lanman, and Kaan Akşit. Optimizing vision and visuals: lectures on cameras, displays and perception. In ACM SIGGRAPH 2022 Courses, pages 1–66. 2022. 

  4. Jason D Schmidt. Numerical simulation of optical wave propagation with examples in matlab. (No Title), 2010. 

  5. John C Heurtley. Scalar rayleigh–sommerfeld and kirchhoff diffraction integrals: a comparison of exact evaluations for axial points. JOSA, 63(8):1003–1008, 1973. 

  6. Maciej Sypek. Light propagation in the fresnel region. new numerical approach. Optics communications, 116(1-3):43–48, 1995. 

  7. Kyoji Matsushima and Tomoyoshi Shimobaba. Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields. Optics express, 17(22):19662–19673, 2009. 

  8. Wenhui Zhang, Hao Zhang, and Guofan Jin. Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range. Optics letters, 45(6):1543–1546, 2020. 

  9. Wenhui Zhang, Hao Zhang, and Guofan Jin. Adaptive-sampling angular spectrum method with full utilization of space-bandwidth product. Optics Letters, 45(16):4416–4419, 2020. 

  10. Yifan Peng, Suyeon Choi, Nitish Padmanaban, and Gordon Wetzstein. Neural holography with camera-in-the-loop training. ACM Transactions on Graphics (TOG), 39(6):1–14, 2020. 

  11. Praneeth Chakravarthula, Ethan Tseng, Tarun Srivastava, Henry Fuchs, and Felix Heide. Learned hardware-in-the-loop phase retrieval for holographic near-eye displays. ACM Transactions on Graphics (TOG), 39(6):1–18, 2020. 

  12. Koray Kavaklı, Hakan Urey, and Kaan Akşit. Learned holographic light transport. Applied Optics, 61(5):B50–B55, 2022. 

  13. Gary J Swanson. Binary optics technology: the theory and design of multi-level diffractive optical elements. Technical Report, MASSACHUSETTS INST OF TECH LEXINGTON LINCOLN LAB, 1989. 

  14. Herwig Kogelnik. Coupled wave theory for thick hologram gratings. Bell System Technical Journal, 48(9):2909–2947, 1969. 

  15. Lingling Huang, Shuang Zhang, and Thomas Zentgraf. Metasurface holography: from fundamentals to applications. Nanophotonics, 7(6):1169–1190, 2018. 

  16. Koray Kavaklı, Yuta Itoh, Hakan Urey, and Kaan Akşit. Realistic defocus blur for multiplane computer-generated holography. In 2023 IEEE Conference Virtual Reality and 3D User Interfaces (VR), 418–426. IEEE, 2023.