5.3. Heisenberg Uncertainty Ray Bending

5.3.1. Overview

Heisenberg Uncertainty Ray Bending (HURB), also known as Heisenberg Uncertainty Monte Carlo (HUMC), is a statistical method designed to approximate edge diffraction within Monte Carlo raytracing. This approach perturbs the direction of rays near aperture edges. The magnitude of this directional deviation is determined by the ray’s proximity to the edge and is governed by the position-momentum uncertainty principle.

The concept of applying this principle to raytracing was initially proposed by Heinisch and Chou [1], with subsequent modifications introduced by Likeness [2] and Coffey et al. [3]. Numerous studies have demonstrated HURB’s effectiveness in approximating both near-field and far-field diffraction phenomena [1][4][5].

HURB has been implemented in commercial raytracers such as TracePro [6], as well as open-source projects like ISET3D [7]. While extensions exist that incorporate the light’s phase to model interference effects [4][5][8], these are outside the scope of this work.

5.3.2. Directional Distribution

The principles of position and momentum can be analyzed by separating them into their respective \(x\) and \(y\) components. For the \(x\)-direction, the uncertainty principle is expressed as [6]:

(5.68)\[\Delta x \Delta p_x \geq \frac{\hbar}{2}\]

In this equation, \(\Delta x\) represents the distance from the ray to the nearest aperture edge, \(\Delta p_x\) is the momentum uncertainty in the \(x\)-direction, and \(\hbar\) denotes the reduced Planck constant.

The relative directional uncertainty can be characterized by a tangent relationship [6]:

(5.69)\[\tan \sigma_x = \frac{\Delta p_x}{p}\]

Here, \(p\) is the total momentum, defined as:

(5.70)\[p = \hbar k = \hbar \frac{2 \pi}{\lambda} = \frac{n(\lambda_0) h}{\lambda_0}\]

Substituting Equation (5.68) and Equation (5.70) into Equation (5.69) yields the following expression for the tangent of the uncertainty in the \(x\)-direction:

(5.71)\[\tan \sigma_x = \frac{1}{2\Delta x k}\]

Similarly, for the \(y\)-direction, the tangent of the uncertainty is given by:

(5.72)\[\tan \sigma_y = \frac{1}{2\Delta y k}\]

A key assumption in this model is that the tangent of the perturbed angles, \(\tan \theta_x\) and \(\tan \theta_y\), follow a bivariate normal distribution with standard deviations \(\tan \sigma_x\) and \(\tan \sigma_y\) respectively [3][5]. The probability density function is therefore:

(5.73)\[p\left(\tan \theta_x, \tan \theta_y\right)=\frac{1}{2 \pi \tan \sigma_x \tan \sigma_y} \exp\left(-\frac{\tan^2 \theta_x}{2 \tan^2 \sigma_x} -\frac{\tan^2 \theta_y}{2 \tan^2 \sigma_y}\right)\]

It is worth noting that some authors [1][6] propose that the angles themselves, rather than their tangents, are normally distributed. For small angles, the small angle approximation renders both approaches equivalent. However, for larger angles, the tangent-based approach is considered more accurate for two primary reasons:

1. Using the tangent intrinsically limits the new ray direction to the forward half-space.

2. From a physical perspective, the uncertainty lies in \(\Delta p_x\) and \(\Delta p_y\). Since momentum \(p\) is constant, the uncertainty in direction is proportional to the tangent of the angle (as shown in Equation (5.69)), rather than the angle directly.

5.3.3. Aperture Edge Distances

Overview

The fundamental procedure, as outlined by Heinisch and Chou [1], involves identifying the direction associated with the minimum distance to an aperture edge. This shortest distance defines the minor axis of a conceptual distance ellipse. Conversely, the major axis of this ellipse is determined by fitting the largest possible area ellipse within the aperture boundaries.

These ellipse axes will be denoted as \(a\) and \(b\), and their corresponding uncertainties as \(\Delta a\) and \(\Delta b\).

Rectangle

Consider a rectangular aperture with a width of \(2B\) and a height of \(2A\), centered at coordinates \((x_0, y_0, z_0)\) and orthogonal to the z-axis. For a ray traversing the aperture at position \((x, y, z_0)\), the distances to the respective edges are given by [6]:

(5.74)\[\begin{split}\Delta a &= \lvert y - y_0 \rvert - A\\ \Delta b &= \lvert x - x_0 \rvert - B\end{split}\]

Fig. 5.11 Example for the HURB ellipse inside a rectangle aperture. The ray position is marked in cyan color, the ellipse in orange.

The direction vectors aligned with the edges are equivalent to the unit vectors in the x and y directions:

(5.75)\[\begin{split}b = \left(\begin{array}{c} 1 \\ 0 \\ 0 \end{array}\right)\end{split}\]

(5.76)\[\begin{split}a = \left(\begin{array}{c} 0 \\ 1 \\ 0 \end{array}\right)\end{split}\]

For a rotated rectangular aperture, the ray’s position must first be transformed into the rotated rectangle’s local coordinate system via a rotation operation. Conversely, the direction vectors \(a\) and \(b\) must be converted from the rotated system back into the absolute coordinate system.

Circle

For a circular aperture with radius \(R\) centered at \((x_0, y_0, z_0)\), the relative polar coordinates of a ray position \((x, y, z_0)\) within the aperture’s coordinate system are defined as:

(5.77)\[\begin{split}r &= \sqrt{(x-x_0)^2 + (y-y_0)^2}\\ \phi &= \text{atan2}(y-y_0, x-x_0)\\\end{split}\]

Fig. 5.12 Example for the HURB ellipse inside a circular aperture. The ray position is marked in cyan color, the ellipse in orange.

The minor axis length of the ellipse, \(\Delta b\), is determined by the distance from the ray to the circular edge:

(5.78)\[\Delta b = R - r\]

The major axis length, \(\Delta a\), is calculated by matching the curvature of the ellipse to that of the aperture. This relationship leads to the formula \(R = \Delta a^2 / \Delta b\) (as derived here), from which \(\Delta a\) is obtained:

(5.79)\[\Delta a = \sqrt{\Delta b R}\]

The unit vectors defining the directions of these axes are:

(5.80)\[\begin{split}b = \left(\begin{array}{c} \cos \phi \\ \sin \phi \\ 0 \end{array}\right)\end{split}\]

(5.81)\[\begin{split}a = \left(\begin{array}{c} \sin \phi \\ -\cos \phi \\ 0 \end{array}\right)\end{split}\]

It is noteworthy that the implementation within ISET3D, as observed here, does not utilize the largest-ellipse approach for defining the major axis. Instead, it directly employs the distance to the circular edge as the size of the major ellipse axis. Compared to the presented method, this alternative approach results in smaller uncertainties of \(\tan \sigma_a\).

5.3.4. New Direction Vectors

Composing a new direction vector for light propagating perpendicular to the aperture plane is a straightforward process. However, this becomes significantly more challenging when considering arbitrary ray angles.

A tilted incident ray effectively “sees” a projected aperture, which introduces several complexities:

• The direction vectors \(a\) and \(b\) are no longer orthogonal to the ray’s propagation vector \(s\).

• The calculated \(\Delta a\) and \(\Delta b\) no longer represent the true shortest distances to the aperture edges.

• The actual shortest distance to the aperture edge does not necessarily lie within the aperture plane.

• The vectors representing the actual shortest distances are generally not orthogonal to each other.

• These actual shortest distances can occur at different starting points along the ray’s path.

Our current implementation aims to address the first three of these issues. We define new orthogonal basis vectors, \(a'\) and \(b'\), relative to the ray’s direction vector \(s\):

(5.82)\[\begin{split}a' &= \frac{b \times s}{\lvert b \times s\rvert}\\ b' &= \frac{s \times a'}{\lvert s \times a'\rvert}\end{split}\]

In this new basis, the vectors \(s\), \(b\), and \(b'\) are coplanar, and \(b'\) aligns with the direction of the shortest distance to the aperture edge. While \(s\) and \(a'\) also lie within a plane, this plane is orthogonal to the one containing \(s, b, b'\). It is important to note that the original vector \(a\) typically does not lie within this second plane. This non-alignment is a consequence of the fact that the true shortest distances to the aperture edges are typically not orthogonal to each other in 3D space when the ray is tilted.

Fig. 5.13 Example of a tilted ray passing through a slit. The ray is tilted around \(a\), so \(a' = a\).

The lengths of the actual shortest distances are adjusted by the cosine of the angles between the ray direction and the original axis vectors:

(5.83)\[\begin{split}\Delta a' &= \Delta a \cos \psi_a\\ \Delta b' &= \Delta b \cos \psi_b\\\end{split}\]

Here, \(\psi_a\) and \(\psi_b\) represent the angles lying in the planes defined by \(s, a\) and \(s, b\) respectively. Assuming \(s\), \(a\), and \(b\) are unit vectors, the cosine terms can be expressed through a dot product, utilizing the identity \(\cos \psi = \sqrt{1 - \sin^2 \psi}\):

(5.84)\[\begin{split}\cos \psi_a &= \sqrt{1 - (s \cdot a)^2}\\ \cos \psi_b &= \sqrt{1 - (s \cdot b)^2}\\\end{split}\]

This leads to modified tangent expressions for the directional uncertainties:

(5.85)\[\begin{split}\tan \sigma_{a'} &= \frac{1}{2 k \Delta a \cos \psi_a}\\ \tan \sigma_{b'} &= \frac{1}{2 k \Delta b \cos \psi_b}\\\end{split}\]

These modified tangents, \(\tan \sigma_{a'}\) and \(\tan \sigma_{b'}\), are then used to determine the values \(\tan \theta_{b'}\) and \(\tan \theta_{a'}\). With this approach, the distances \(\Delta a'\) and \(\Delta b'\), as well as the vector \(b'\), are calculated correctly. However, the vector \(a'\) is generally not entirely accurate due to the non-orthogonality of the true shortest distance vectors for tilted rays. An exception occurs when the ray is tilted exclusively along the \(a\) or \(b\) direction, in which case both vectors \(a'\) and \(b'\) are correctly determined.

The new perturbed direction vector is constructed by summing the initial ray direction with components scaled by \(\tan \theta_{b'}\) and \(\tan \theta_{a'}\) along the \(b'\) and \(a'\) unit vectors, respectively:

(5.86)\[s'_\text{un} = s + b' \tan \theta_{b'} + a' \tan \theta_{a'}\]

This resultant vector, \(s'_\text{un}\), must then be normalized to obtain the final perturbed unity direction vector \(s'\):

(5.87)\[s' = \frac{s'_\text{un}}{\lvert s'_\text{un} \rvert}\]

The implementation in ISET3D (as seen here) constructs the new vector directly from \(s\), \(a\), and \(b\). This method can introduce errors for significant ray tilts, as the components \(a\) and \(b\) are no longer orthogonal to each other in the projected plane. While our implementation addresses some of the aforementioned issues, it does not fully resolve all complexities arising from tilted rays. Nevertheless, the residual issues typically manifest noticeably only for rays with large tilt angles, as the projection effects exhibit a cosine dependency that can be neglected for small angles.

5.3.5. Selecting the Uncertainty Factor

To provide further parametrization, a scaling factor \(\gamma\) has been introduced into the uncertainty equations (5.71) and (5.72). The modified equations become:

(5.88)\[\begin{split}\tan \sigma_{a'} &= \frac{\gamma}{2 k \Delta a \cos \psi_a}\\ \tan \sigma_{b'} &= \frac{\gamma}{2 k \Delta b \cos \psi_b}\\\end{split}\]

Most authors typically employ a value of \(\gamma = 1\). However, Lin (2015) proposed a scaling factor of \(\gamma = \sqrt{2}\) ([9], page 17) to achieve a better fit for diffraction patterns observed at a circular aperture. This value is also implemented in the PBRT fork used by ISET3D, as can be seen here. It’s worth noting that while these sources might present the fraction as \(\tan \sigma = \frac{1}{\sqrt{2} \ldots}\), this is mathematically equivalent to our formulation of \(\tan \sigma = \frac{\gamma}{2\ldots}\) with \(\gamma = \sqrt{2}\) in equation (5.88).

In optrace, \(\gamma\) can be customized via a parameter, as detailed in Section Section 4.5.3, allowing for fine-tuning to specific simulation scenarios. By default, optrace uses \(\gamma = \sqrt{2}\). The results of this parametrization can be found in the subsequent section.

A value of \(\gamma = 2\) has been observed to provide a superior fit for the main lobe in comparison to graphs generated with \(\gamma = 1\), which are commonly found in existing literature. With \(\gamma = 2\), the simulation asymptotically matches the side lobes for slit, pinhole, and pupil geometries. Conversely, for edge diffraction, the fit is less optimal, producing more pronounced glare to the left of the edge. In such cases, \(\gamma = 1\) is likely to yield better results.

5.3.6. Comparison to Theoretical Curves

The following curves were generated using the test scripts located under tests/hurb/. All simulations were conducted with an uncertainty factor of \(\gamma = \sqrt{2}\).

Fig. 5.14 Comparison of HURB and the theoretical curve for the far field diffraction of a slit.

Fig. 5.15 Comparison of HURB and the theoretical curve for the far-field diffraction of a pinhole.

Fig. 5.16 Comparison of HURB and the theoretical curve for the focus profile of an ideal lens behind a pupil.

Fig. 5.17 Comparison of HURB and the theoretical curve for the profile of edge diffraction.

---

References

[1] (1,2,3,4)

R. P. Heinisch and T. S. Chou. Numerical experiments in modeling diffraction phenomena. Applied Optics, 10(10):2248, October 1971. URL: http://dx.doi.org/10.1364/AO.10.002248, doi:10.1364/ao.10.002248.

[2]

Barry K. Likeness. Stray light simulation with advanced monte carlo techniques. In John D. Lytle and Howard E. Morrow, editors, Stray Light Problems in Optical Systems, volume 0107, 80–88. SPIE, September 1977. URL: http://dx.doi.org/10.1117/12.964599, doi:10.1117/12.964599.

[3] (1,2)

Katherine L. Coffey, Kory J. Priestley, J. Robert Mahan, and Maria C. Sanchez. Diffraction models of radiation entering an aperture for use in a monte carlo ray-trace environment. In Robert E. Fischer, Lawrence M. Germann, Alson E. Hatheway, Malachy McConnell, and Warren J. Smith, editors, Current Developments in Optical Design and Engineering VII, volume 3429, 213. SPIE, October 1998. URL: http://dx.doi.org/10.1117/12.328547, doi:10.1117/12.328547.

[4] (1,2)

Edward R. Freniere. Numerical experiments in modeling diffraction phenomena with monte carlo ray tracing. In José M. Sasian and Mary G. Turner, editors, Novel Optical Systems Design and Optimization IX, volume 6289, 62890N. SPIE, August 2006. URL: http://dx.doi.org/10.1117/12.681797, doi:10.1117/12.681797.

[5] (1,2,3)

Vladimir V. Serikov and Shinji Kawamoto. Numerical experiments in monte carlo modeling of polarization, diffraction, and interference phenomena. In Frank Wyrowski, editor, Wave-Optical Systems Engineering, volume 4436, 80–88. SPIE, December 2001. URL: http://dx.doi.org/10.1117/12.451286, doi:10.1117/12.451286.

[6] (1,2,3,4,5)

Edward R. Freniere, G. Groot Gregory, and Richard A. Hassler. Edge diffraction in monte carlo ray tracing. In Richard C. Juergens, editor, Optical Design and Analysis Software, volume 3780, 151–157. SPIE, September 1999. URL: http://dx.doi.org/10.1117/12.363773, doi:10.1117/12.363773.

[7]

Trisha Lian, Kevin J. MacKenzie, David H. Brainard, Nicolas P. Cottaris, and Brian A. Wandell. Ray tracing 3d spectral scenes through human optics models. Journal of Vision, 19(12):23, October 2019. URL: http://dx.doi.org/10.1167/19.12.23, doi:10.1167/19.12.23.

[8]

Yijia Yuan, Kuanfang Ren, and Claude Rozé. Fraunhofer diffraction of irregular apertures by heisenberg uncertainty monte carlo model. Particuology, 24:151–158, February 2016. URL: http://dx.doi.org/10.1016/j.partic.2015.03.006, doi:10.1016/j.partic.2015.03.006.

[9]

Andy Lai Lin. The Computational Image Systems Evaluation Toolbox: a software environment for rapid prototyping and evaluating computational imaging hardware and algorithms. PhD thesis, Stanford University, Stanford, CA, 2015. Available online via Stanford Digital Repository. URL: https://purl.stanford.edu/xz083nw1432.