scico.linop.optics#

Optical propagator classes.

This module provides classes that model the propagation of a monochromatic waveform between two parallel planes in a homogeneous medium. The corresponding linear operators are referred to here as “propagators”, which represents a departure from standard terminology, in which “propagator” refers specifically to the Fourier domain component of the linear operator, i.e. if the full linear operator can be written as \(F^{-1} D F\) where \(F\) is the Fourier transform, then \(D\) is usually referred to as the propagator.

The following notation is used throughout the module:

\[\begin{split}\begin{align} \Delta x, \Delta y & \quad \text{Sampling intervals in } x \text{ and } y \text{ axes}\\ z & \quad \text{Propagation distance} \;\; (z \geq 0) \\ N_x, N_y & \quad \text{Number of samples in } x \text{ and } y \text{ axes}\\ k_0 & \quad \text{Illumination wavenumber corresponding to } 2\pi / \text{wavelength} \;. \end{align}\end{split}\]

Variables \(\Delta x, \Delta y, z,\) and \(k_0\) represent physical quantities. Any units may be chosen, but they must be consistent across all of these variables, e.g. m (metres) for \(\Delta x, \Delta y, z,\) and \(\mathrm{m}^{-1}\) for \(k_0\), as well as with the units for the physical dimensions of the source wavefield.

Subscripts \(S\) and \(D\) are used to refer to the source and destination planes respectively when it is necessary to distinguish between them. In the absence of subscripts, the variables refer to the source plane (e.g. both \(\Delta x\) and \(\Delta x_S\) refer to the \(x\)-axis sampling interval in the source plane, while \(\Delta x_D\) refers to it in the destination plane).

Note that \(x\) corresponds to axis 0 (rows, increasing downwards) and \(y\) to axis 1 (columns, increasing to the right).

Functions

radial_transverse_frequency(input_shape, dx)

Construct radial Fourier coordinate system.

Classes

`AngularSpectrumPropagator`(input_shape, dx, k0, z)	Angular spectrum propagator.
`FraunhoferPropagator`(input_shape, dx, k0, z)	Fraunhofer (far-field) propagator.
`FresnelPropagator`(input_shape, dx, k0, z[, ...])	Fresnel (small-angle/paraxial) propagator.
`Propagator`(input_shape, dx, k0, z[, pad_factor])	Base class for angular spectrum and Fresnel propagators.

scico.linop.optics.radial_transverse_frequency(input_shape, dx)[source]#

Construct radial Fourier coordinate system.

Parameters:

input_shape (Tuple[int, ...]) – Tuple of length 1 or 2 containing the number of samples per dimension, i.e. \((N_x,)\) or \((N_x, N_y)\)
dx (Union[float, Tuple[float, ...]]) – Sampling interval at source plane. If a float and len(input_shape)==2 the same sampling interval is applied to both dimensions. If dx is a tuple, it must have same length as input_shape, and corresponds to either \((\Delta x,)\) or \((\Delta x, \Delta y)\).

Return type:

ndarray

Returns:

If len(input_shape)==1, returns an ndarray containing corresponding Fourier coordinates. If len(input_shape) == 2, returns an ndarray containing the radial Fourier coordinates \(\sqrt{k_x^2 + k_y^2}\,\).

class scico.linop.optics.Propagator(input_shape, dx, k0, z, pad_factor=1, **kwargs)[source]#

Bases: LinearOperator

Base class for angular spectrum and Fresnel propagators.

Parameters:

input_shape (Tuple[int, ...]) – Shape of input array as a tuple of length 1 or 2, corresponding to \((N_x,)\) or \((N_x, N_y)\).
dx (Union[float, Tuple[float, ...]]) – Sampling interval at source plane. If a float and len(input_shape)==2 the same sampling interval is applied to both dimensions. If dx is a tuple, it must have same length as input_shape, and corresponds to either \((\Delta x,)\) or \((\Delta x, \Delta y)\).
k0 (float) – Illumination wavenumber, \(k_0\), corresponding to \(2 \pi\) / wavelength.
z (float) – Propagation distance, \(z\).
pad_factor (int) – The padded input shape is the input shape multiplied by this integer factor.

k0: float#: Illumination wavenumber; 2𝜋/wavelength

padded_shape: Tuple[int, ...]#: Shape of input after padding

L: Tuple[float, ...]#: Padded source plane side length (dx[i] * padded_shape[i])

kp#: Transverse Fourier coordinates (radial)

dx: Union[float, Tuple[float, ...]]#: Source plane sampling interval

z: float#: Propagation distance

class scico.linop.optics.AngularSpectrumPropagator(input_shape, dx, k0, z, pad_factor=1, jit=True, **kwargs)[source]#

Bases: Propagator

Angular spectrum propagator.

Propagates a planar source field with coordinates \((x, y, z_0)\) to a destination plane at a distance \(z\) with coordinates \((x, y, z_0 + z)\). The action of this linear operator is given by (Eq. 3.74, [27])

\[(A \mb{u})(x, y, z_0 + z) = \frac{1}{2 \pi} \iint_{-\infty}^{\infty} \mb{\hat{u}}(k_x, k_y) e^{j \sqrt{k_0^2 - k_x^2 - k_y^2} \, z} e^{j (x k_x + y k_y) } d k_x \ d k_y \;,\]

where the \(\mb{\hat{u}}\) is the Fourier transform of the field \(\mb{u}(x, y)\) in the plane \(z=z_0\), given by

\[\mb{\hat{u}}(k_x, k_y) = \iint_{-\infty}^{\infty} \mb{u}(x, y) e^{- j (x k_x + y k_y)} d k_x \ d k_y \;,\]

where \((k_x, k_y)\) are the \(x\) and \(y\) components respectively of the wave-vector of the plane wave, and \(j\) is the imaginary unit.

The angular spectrum propagator can be written

\[A\mb{u} = F^{-1} D F \mb{u} \;,\]

where \(F\) is the Fourier transform with respect to \((x, y)\), \(F^{-1}\) is the inverse transform with respect to \((k_x, k_y)\), and the propagator term is given by

\[D = \exp \left( j \sqrt{k_0^2 - k_x^2 - k_y^2} \, z \right) \;.\]

Aliasing of the wavefield at the destination plane is avoided when the propagator term is adequately sampled according to [52]

\[(\Delta x)^2 \geq \frac{\pi}{k_0 N_x} \sqrt{ (\Delta x)^2 N_x^2 + 4 z^2} \quad \text{and} \quad (\Delta y)^2 \geq \frac{\pi}{k_0 N_y} \sqrt{ (\Delta y)^2 N_y^2 + 4 z^2} \;.\]

Parameters:

input_shape (Tuple[int, ...]) – Shape of input array. Can be a tuple of length 2 or 3.
dx (Union[float, Tuple[float, ...]]) – Sampling interval, \(\Delta x\), at source plane. If a float and len(input_shape)==2 the same sampling interval is applied to both dimensions. If dx is a tuple, must have same length as input_shape.
k0 (float) – Illumination wavenumber, \(k_0\), corresponding to \(2 \pi\) / wavelength.
z (float) – Propagation distance, \(z\).
pad_factor (int) – The padded input shape is the input shape multiplied by this integer factor.
jit (bool) – If True, call jit on this LinearOperator to jit the forward, adjoint, and gram functions. Same as calling jit after the LinearOperator is created.

adequate_sampling()[source]#

Verify the angular spectrum kernel is not aliased.

Checks the condition for adequate sampling [52],

\[(\Delta x)^2 \geq \frac{\pi}{k_0 N_x} \sqrt{ (\Delta x)^2 N_x^2 + 4 z^2} \quad \text{and} \quad (\Delta y)^2 \geq \frac{\pi}{k_0 N_y} \sqrt{ (\Delta y)^2 N_y^2 + 4 z^2} \;.\]

Returns:: True if the angular spectrum kernel is adequately sampled, False otherwise.

pinv(y)[source]#: Apply pseudoinverse of Angular Spectrum propagator.

k0: float#: Illumination wavenumber; 2𝜋/wavelength

padded_shape: Tuple[int, ...]#: Shape of input after padding

L: Tuple[float, ...]#: Padded source plane side length (dx[i] * padded_shape[i])

dx: Union[float, Tuple[float, ...]]#: Source plane sampling interval

z: float#: Propagation distance

class scico.linop.optics.FresnelPropagator(input_shape, dx, k0, z, pad_factor=1, jit=True, **kwargs)[source]#

Bases: Propagator

Fresnel (small-angle/paraxial) propagator.

\[(A \mb{u})(x, y, z + z_0) = e^{j k_0 z} \frac{1}{2 \pi} \iint_{-\infty}^{\infty} \mb{\hat{u}}(k_x, k_y) e^{-j \frac{z}{2 k_0}\left(k_x^2 + k_y^2\right) } e^{j (x k_x + y k_y) } d k_x \ d k_y \;,\]

where the \(\mb{\hat{u}}\) is the Fourier transform of the field in the source plane, given by

\[\mb{\hat{u}}(k_x, k_y) = \iint_{-\infty}^{\infty} \mb{u}(x, y) e^{- j (x k_x + y k_y)} d k_x \ d k_y \;.\]

This linear operator is valid when \(k_0^2 << k_x^2 + k_y^2\). The Fresnel propagator can be written

\[A\mb{u} = F^{-1} D F \mb{u} \;,\]

where \(F\) is the Fourier transform with respect to \((x, y)\), \(F^{-1}\) is the inverse transform with respect to \((k_x, k_y)\), and the propagator term is given by

\[D = \exp \left( -j \frac{z}{2 k_0}\left(k_x^2 + k_y^2 \right) \right) \;,\]

where \((k_x, k_y)\) are the \(x\) and \(y\) components respectively of the wave-vector of the plane wave, and \(j\) is the imaginary unit.

The propagator term is adequately sampled when [51]

\[(\Delta x)^2 \geq \frac{2 \pi z }{k_0 N_x} \quad \text{and} \quad (\Delta y)^2 \geq \frac{2 \pi z }{k_0 N_y} \;.\]

Parameters:

input_shape (Tuple[int, ...]) – Shape of input array as a tuple of length 1 or 2, corresponding to \((N_x,)\) or \((N_x, N_y)\).
dx (float) – Sampling interval at source plane. If a float and len(input_shape)==2 the same sampling interval is applied to both dimensions. If dx is a tuple, it must have same length as input_shape, and corresponds to either \((\Delta x,)\) or \((\Delta x, \Delta y)\).
k0 (float) – Illumination wavenumber, \(k_0\), corresponding to \(2 \pi\) / wavelength.
z (float) – Propagation distance, \(z\).
pad_factor (int) – The padded input shape is the input shape multiplied by this integer factor.

adequate_sampling()[source]#

Verify the Fresnel propagation kernel is not aliased.

Checks the condition for adequate sampling [51],

\[(\Delta x)^2 \geq \frac{2 \pi z }{k_0 N_x} \quad \text{and} \quad (\Delta y)^2 \geq \frac{2 \pi z }{k_0 N_y} \;.\]

Returns:: True if the Fresnel propagation kernel is adequately sampled, False otherwise.

k0: float#: Illumination wavenumber; 2𝜋/wavelength

padded_shape: Tuple[int, ...]#: Shape of input after padding

L: Tuple[float, ...]#: Padded source plane side length (dx[i] * padded_shape[i])

dx: Union[float, Tuple[float, ...]]#: Source plane sampling interval

z: float#: Propagation distance

class scico.linop.optics.FraunhoferPropagator(input_shape, dx, k0, z, jit=True, **kwargs)[source]#

Bases: LinearOperator

Fraunhofer (far-field) propagator.

Propagates a source field with coordinates \((x_S, y_S)\) to a destination plane at a distance \(z\) with coordinates \((x_D, y_D)\).

The action of this linear operator is given by (Eq. 4.25, [27])

\[(A \mb{u})(x_D, y_D) = \underbrace{\frac{k_0}{2 \pi} \frac{e^{j k_0 z}}{j z} \mathrm{exp} \left( j \frac{k_0}{2 z} (x_D^2 + y_D^2) \right)}_{\triangleq P(x_D, y_D)} \int \mb{u}(x_S, y_S) e^{-j \frac{k_0}{z} (x_D x_S + y_D y_S) } dx_S \ dy_S \;.\]

This is valid when \(N_F << 1\), where \(N_F\) is the Fresnel number (Sec. 1.5, Sec. 4.7.2.1) [40]. Writing the Fourier transform of the field \(\mb{u}\) as

\[\hat{\mb{u}}(k_x, k_y) = \int e^{-j (k_x x + k_y y)} \mb{u}(x, y) dx \ dy \;,\]

the action of this linear operator can be written

\[(A \mb{u})(x_D, y_D) = P(x_D, y_D) \ \hat{\mb{u}} \left({\frac{k_0}{z} x_D, \frac{k_0}{z} y_D}\right) \;.\]

Ignoring multiplicative prefactors, the Fraunhofer propagated field is the Fourier transform of the source field, evaluated at coordinates \((k_x, k_y) = (\frac{k_0}{z} x_D, \frac{k_0}{z} y_D)\).

In general, the sampling intervals (and thus plane lengths) differ between source and destination planes. In particular, (Eq. 5.18, [51])

\[\Delta x_D = \frac{2 \pi z}{k_0 L_{Sx} } \quad \text{and} \quad L_{Dx} = \frac{2 \pi z}{k_0 \Delta x_S } \;,\]

and similarly for the \(y\) axis.

The Fraunhofer propagator term \(P(x_D, y_D)\) is adequately sampled when

\[\Delta x_S \geq \sqrt{\frac{2 \pi z}{N_x k_0}} \quad \text{and} \quad \Delta y_S \geq \sqrt{\frac{2 \pi z}{N_y k_0}} \;.\]

Parameters:

input_shape (Tuple[int, ...]) – Shape of input array as a tuple of length 1 or 2, corresponding to \((N_x,)\) or \((N_x, N_y)\).
dx (Union[float, Tuple[float, ...]]) – Sampling interval at source plane. If a float and len(input_shape)==2 the same sampling interval is applied to both dimensions. If dx is a tuple, it must have same length as input_shape, and corresponds to either \((\Delta x,)\) or \((\Delta x, \Delta y)\).
k0 (float) – Illumination wavenumber, \(k_0\), corresponding to \(2 \pi\) / wavelength.
z (float) – Propagation distance, \(z\).
jit (bool) – If True, jit the evaluation, adjoint, and gram functions of this LinearOperator. Default: True.

k0: float#: Illumination wavenumber

z: float#: Propagation distance

L: Tuple[float, ...]#: Source plane side length (dx[i] * input_shape[i])

dx: Tuple[float, ...]#: Source plane sampling interval

dx_D: Tuple[float, ...]#: Destination plane sampling interval

L_D: Tuple[float, ...]#: Destination plane side length

adequate_sampling()[source]#

Verify the Fraunhofer propagation kernel is not aliased.

Checks the condition for adequate sampling [51],

\[\Delta x_S \geq \sqrt{\frac{2 \pi z}{N_x k_0}} \quad \text{and} \quad \Delta y_S \geq \sqrt{\frac{2 \pi z}{N_y k_0}} \;.\]

Returns:: True if the Fraunhofer propagation kernel is adequately sampled, False otherwise.