1. Introduction

Since the past decade, the production methods of the semiconductor industry allow the structuring of surfaces in the nanometer scale of the visible wavelengths and below. This capability introduced a new type of optical components under the generic term of diffractive optics. These include foliations, gratings, diffractive lenses [1] and more. One of the main benefits of diffractive components is the introduction of new concepts of optical functionality, which are not attainable with conventional optical systems, while minimizing the system’s size. However, these structures may have an impact on the overall part’s geometry and thus aesthetic design (e.g. in headlamps) resulting in unintended color fringes or reflections.

The concept of rendering represents the key approach for the determination of an optical component’s light interaction with its surrounding scenery. The basic principle of rendering is to define all elements of the scenery as geometric objects and calculate the intersections of a geometric line (ray) with the scenery’s objects. Depending on the luminance representation, a scalar factor or a wavelength spectrum as well as the ray’s outgoing direction are calculated at the intersection point in correlation to the geometry’s material. With further determined intersections a ray path through the scenery is defined. Therefore, this process is commonly known as “path tracing”. One characteristic aspect of path tracing is that, for computation efficiency, the rays are propagated inversely from the camera to the source, which is referred to as ”backward raytracing” [2] .

The main challenge is now to render the diffractive optical elements within their later scenery. To achieve this goal a method for defining the wave-optical effects induced by the diffractive component needs to be introduced into the rendering algorithm. However, with the already mentioned scale of the diffractive element’s surface structure, complex field approaches, which fully solve Maxwell’s equations, are necessary to define the light’s interaction with these surfaces, instead of applying the much simpler geometric path tracing approach.

There is a wide variety of simulation approaches already available for a wave-optical solution. These include the Finite Difference Time Domain (FDTD - [3]), Rigorous Coupled Wave Analysis (RCWA - [4]), Beam Propagation Method (BPM - [5]) or the Method of Lines (MOL - [6]). However, in [7], [8] and [9] it is pointed out that a full-field calculation at runtime for every intersection point would cause a non manageable amount of computing time for rendering macro-scale sceneries.

Therefore, several authors [789789789] have developed approximate approaches for simulating wave-optical effects and combining them with path tracing, whose reduced computational effort allows these methods to be performed at the corresponding surface intersection point within the rendering algorithm. In [7], a scattering function is developed by the use of the Wigner distribution function to account for the phase propagation and therefore the wave interference. This function is only valid for a paraxial approach and is included in a state of the art rendering algorithm. It approximates rudimentary wave effects to visualize some of the color artifacts of an optical disk. In [8], wave-optical effects are induced by solving the Helmholtz equation in planar slices and using plane wave expansions to visualize randomly distributed scratches on metallic surfaces. The approach is strongly limited in geometric freedom of the diffractive structure. In [9], effects are respected by weighted samples within a spectral rendering process along with random point sampling. The key criterion for evaluating these rendering methods is to compare color shapes and gradients induced by wave-optical effects of e.g a Compact Disk (CD) between the rendering result and the corresponding photograph of a test scenery. The published approximation methods differ significantly from the corresponding photographs in a way which makes them not usable for the simulation of more complex diffractive elements like foliations and gratings, since these elements need to be simulated by full-field approaches.

There have also been first attempts to unify the backward ray-tracing concept in rendering with rigorous diffractive simulation in which the wave-optics simulation has been performed through FDTD [10,11] or scalar diffraction theory [12]. These methods enable the rendering of diffractive elements such as metasurfaces or objects with structural color (e.g. the Morpho butterfly) and similarly treat the diffractive simulation as a surface material model within raytracing. However, there are subtle but important limitations: In the case of [10], the authors only model light incidence within one plane of incidence in their diffraction model, leading to irregularities for three-dimensional rendering. In [11], the authors simplify their material model by rescaling and bicubic interpolation of samples, to reduce the amount of data the model needs. Furthermore, because of the comparably low computing efficiency of FDTD, only a low resolution of the angular sampling is used. In [12], the diffraction simulation is only scalar and not entirely rigorous, due to a ray-based integration of light scattering within the microstructure. This results in limitations as missing knife-edge diffraction within the microstructure or the inability to simultaneously account for reflection and transmission.

On the other hand, several authors combine the RCWA with the concept of forward raytracing and focus mainly on optical systems or component design without path rendering. In these works, the RCWA is used to derive a function or a look-up table describing the spectrally resolved properties of the diffractive surface which is then used in a standard forward raytracing. Investigated properties are for example the diffraction efficiency, the transmittance / reflectance or the scattering behaviour. This idea has applications in the design of high-efficiencey Light Emitting Diodes (LEDs) [13,14], in solar cells [15], in holograms [16] or in projectors [17].

Thus, no backward raytracing for path rendering has up to date been performed using RCWA [18]. As RCWA provides some significant advantages over other methods as FDTD - e.g. a direct assessment of the far field through the diffraction orders, a thorough near field investigation and a higher computation efficiency -, the full potential of rigorous wave simulations has up to date not been exploited. Furthermore, all existing methods combining wave-optical approaches with rendering provide and present only limited possibilities for simulating the optical interaction between different components in complex sceneries.

Therefore, the goal of the presented paper is to develop an inclusion process of a full-field simulation approach based on RCWA into the backward path rendering process. The focus is set on determining and testing the diffractive optical elements’ performance, setup and aesthetic design within its planned application scenery. The ray’s outgoing direction as well as the resulting wavelength spectrum of an intersection should be entirely defined by the field simulation, but should additionally be performed in a Monte Carlo setup to fully use the already available methods for rendering complex scenes. To this end, the field simulation is first performed separately with discrete angles and wavelengths and then the result is coupled to the ray approach within the rendering algorithm. Herein, the specific characteristics of the RCWA simulation of evaluating diffraction orders is exploited for increased efficiency. This tool can be of use in a wide range of capabilities of optical functionality within diffractive optics as well as to omit the immense cost of prototyping in planar technology. The entire design, verification and testing steps along with a virtual image of a wide range of diffractive optical elements can be carried out with the help of the presented methods before manufacturing the first prototype.

This paper is structured as follows: Section 2 gives an overview of the field simulation component of this work followed by the derivation of the developed rendering method. Section 3 presents results and comparisons of a rendered scenery with the correlating photograph. Finally, Section 4 sums up the key elements of the proposed method as well as a few outlooks for adjoining work.

2. Methods

Concerning the rendering component of this work, color and luminance need to be represented in a wavelength dependent form, since the inclusion of wavelength dependent effects is necessary. This type of rendering is commonly known as spectral rendering [2]. The book Physically Based Rendering Technique (pbrt-v3 [2]) by Matt Pharr, Wenzel Jakob and Greg Humphreys and the corresponding freely accessible implementation pbrt-v3 present state of the art methods of the physically formulated interaction of light with arbitrary surfaces and set the base point of the presented method for rendering diffractive structures. The key benefit of pbrt concerning the goal of this work is the possibility to define a wavelength dependent spectrum for each ray.

The challenge is to identify a physically accurate correspondence between the ray model and the field representation. A self-evident possibility is the coupling of the ray as a geometric vector with an amplitude or spectrum to a homogeneous wave within the field simulation. Therefore, a simulation method based on homogeneous wave representations is necessary, where the field interaction can be determined on a periodic structure with its geometry in the scale of the incident wavelength and below.

2.1 Rigorous Coupled Wave Analysis - RCWA

With focus on computing time, Rigorous Coupled Wave Analysis (RCWA - [19]) has been chosen for the field simulation component of the presented method. Especially its direct solution character rather than an iterative one makes it suitable for performing efficient simulations of large input parameter ranges. Figure 1 shows the basic principle of RCWA for a simplified geometry.

 

Fig. 1. The left shows the basic principle of RCWA, where a homogeneous wave of amplitude $E_0$ and wave vector $k_{inc}$ is incident on an infinite 1D-periodic grating with identical cells of period $\Lambda _x$. The graph on the right represents the periodic multilayer cell-representation approximating a CD’s surface as used for verification in Section 3. Here, $n_A$, $n_P$ and $n_M$ indicate the index of refraction for air, polycarbonate and aluminum respectively and $d_1 = d_3 = d_4 = 100 \mathrm {nm}$ and $d_2 = 1 \mathrm {\mu m}$.

Download Full Size | PDF

The excitatory field of this setup consists of a homogeneous wave defined by incidence and polarization angles as well as a wavelength and an input amplitude. At first, the field in all three regions $G_{1,2,3}$ is expressed in terms of Fourier field expansions. The regions $G_1$ and $G_3$ represent infinite half spaces filled with a homogeneous material of refractive index $n_1$ and $n_3$ respectively. The region $G_2$ can be described as a cell of a varying refractive index along the x-axis, that repeats itself infinitely with period $\Lambda _x$ along the x-axis. The geometry in the left of Fig. 1 is simplified for illustration and, as the used RCWA-package follows no restrictions concerning the geometric complexity or the used materials within the grating’s periodic cell in $G_2$, more complex structures, including multilayers, can be considered. The right part of Fig. 1 shows exemplarily a cell decomposition for the layer setup of the CD used in section 3. Inserting the field representations for the regions $G_{1,2,3}$ into Maxwell’s equations results in an eigenvalue problem, whose solution defines the propagation terms of the corresponding field orders. By using the continuity condition of electromagnetic fields on lateral faces, the final result of the RCWA method consists of several diffraction orders, where each order defines an outgoing homogeneous wave with correlating output angles and amplitudes $R_i$ ($G_1$), $T_i$ ($G_3$) and $c_i$ ($G_2$). This basic principle of RCWA has been developed in [19] and further extensions have been published in the following years. These include a multilayer approach [20] to account for the grating’s variation along the z-axis. Furthermore, there are published concepts to handle conical wave incidence [4], 2D gratings [21] as well as several convergence optimizations. In the latter, the concept of normal vector transformation [22] in correlation to the inverse formulation [23] is a key method for the simulation of large input parameter ranges for 2D grating structures. All noted extensions as well as other self developed ones have been implemented into a 3D RCWA Python package. This package is used to perform simulations of the diffractive element’s periodic cell over input parameters ranges of the incidence angles as well as the visible wavelength spectrum. The results are then saved to a local file, which is later read in by the rendering algorithm at runtime. Therefore, in the rendering step, the RCWA’s results need only to be read from a local array instead of performing a full RCWA simulation at runtime. Concerning the incident azimuthal angle, the symmetry of the grating is utilized in the corresponding rendering function.

With the results of the RCWA simulation, there are now direction combinations between the incident wave and the outgoing diffraction orders available to define the ray’s path at the intersection on the diffractive element’s surface. Furthermore, the corresponding diffraction efficiencies represent the transition coefficient for every direction combination between incoming and outgoing ray.

2.2 BWDF derivation

The key challenge of rendering is to determine the subsequent ray direction as well as the corresponding wavelength spectrum after an intersection point between the current ray and the scenery is identified. The concept of the bidirectional distribution function (BxDF - [2]) is a widely used model for meeting these two challenges. In the following, a short introduction to the BxDF is given before a new BxDF concept is described that includes wave-optical effects based on RCWA simulation results.

In the general case of rendering, apart from special formulations like photon mapping [2], the ray path originates at the camera point and propagates towards the light source. This backward raytracing shows key benefits in the efficiency of the rendering algorithm. At the intersection point, the camera’s ray direction $\omega _o$ is flipped to point back towards the camera point. The current camera ray would now propagate further into the scene with the yet unknown direction $\omega _i$. Both directions are formulated within the local coordinate system of the intersecting geometry. It is now physically relevant, which luminance $L_o(p,\omega _o)$ is emitted into the camera’s direction $\omega _o$ at the current intersection point $p$. Let $L_i(p,\omega _i)$ be the yet unknown luminance at the same point $p$ incident from direction $\omega _i$. $L_o(p,\omega _o)$ can be formulated according to [2] by

Calculating this integral constitutes one of the main challenges of physically based rendering. Considering a simple material like a perfectly reflective mirror, only one direction $\omega _i$ in correlation to the given camera direction $\omega _o$ is evident. This follows the basic principle of the perfect reflection of mirrors. However, in the general case, all spatial directions in relation to the unit sphere contribute to the integral in Eq. (1). The Monte Carlo concept [2] is a favored method for the numerical estimation of such higher dimensional integrals. The fundamental principle for estimating the integral in Eq. (1) is to propagate only one ray from the current intersection point. The direction $\omega _i$ is chosen by a probability density function and the resulting luminance factor from direction $\omega _i$ is multiplied by a normalization factor as the integral of the chosen density function. With the use of more rays and intersections, this formulation approximates the exact solution of the integral as an estimation of Eq. (1). This principle is the key point of the luminance’s estimation in this work.

One key problem of formulating the BxDF for wave-optical effects is the already mentioned “backward tracing”. This tracing direction causes a problem for the field inclusion since $L_o(p,\omega _o)$ as well as the direction $\omega _i$ shall be determined by the RCWA, but the input parameters of the necessary RCWA simulation have to be derived from the yet unknown direction $\omega _i$.

In the following, a BxDF for the integration of wave-optical effects is presented. In accordance to pbrt’s implementation and nomenclature, the function is declared as BWDF (bidirectional wave density function). For simplicity, the situation of a merely reflective 1D-grating is considered, where according to Fig. 2 the intersection point of the current ray originating from the camera point is known.

 

Fig. 2. The intersection of a spectral ray with a purely reflective diffraction grating.

Download Full Size | PDF

Its direction $\omega _o$ in the local RCWA coordinate system of the intersected geometry’s surface is represented by the black vector. From the grating’s point of view, this direction must be an arbitrary outgoing diffraction order of the grating structure, for which the corresponding azimuthal $\phi _o$ and polar $\theta _o$ angles are calculated. At this point, it is assumed, that the direction $\omega _o$ corresponds to one out of a total amount of $N_{x,tot}$ considered diffraction orders. By neglecting higher orders with low diffraction efficiency, $N_{x,tot}$ may be smaller than the amount of diffraction orders used in the RCWA simulation itself. Nevertheless, all $N_{x,tot}$ have to be arranged symmetrically around the zeroth order of total reflection. One diffraction order $n_{x,rand}$ of $N_{x,tot}$ is chosen by using a uniformly distributed density function and the direction $\omega _o$ is assigned to this order. Furthermore, one sampled wavelength $\lambda _{rand}$ is arbitrarily selected from the entire discretized spectrum of $L_o(p,\omega _o)$, since the wavelength dependence of the current path can only be expressed for one wavelength. The continuity of electromagnetic fields on boundaries defines consequentially the incident ray (green - $\omega _i$) corresponding to the outgoing ray (black - $\omega _o$). This correspondence can be expressed by the use of the RCWA wave vector components in Eq. (2).

With the knowledge of the incoming direction the discretized spectrum of $f(p,\omega _o,\omega _i)$ needs to be iterated since the direction $\omega _i$ is wavelength dependent. For every iterated wavelength $\lambda _s$ of this spectrum the resulting diffraction order is derived to

The presented BWDF concept can be extended for transmitted field components as well. In this case, the negative z-component of $\omega _o$ in Fig. 2 needs to be considered in the calculation of the corresponding angles. Furthermore, a uniformly distributed choice test at every intersection is necessary to determine if the ray path is reflected or transmitted. Therefore, every coefficient $DE_s$ is additionally multiplied by a factor of 2. Moreover, the structure needs to be simulated over the same input parameter ranges in the RCWA with a z-axis inverted setup. This extended concept is denoted as the “see through" principle since all four possible ray paths of the structure are respected.

In addition, the current BWDF implementation allows the inclusion of 2D periodic structures. In this case, another diffraction order $n_{y,rand}$ of $N_{y,tot}$ is chosen randomly within a uniform distribution at the beginning. The y-component’s equation of Eq. (2) is therefore derived to $k_{y,i} = k_{y,o} + n_{y,rand} (2 \pi / \Lambda _y)$ with the grating period $\Lambda _y$ in the y-axis. The diffraction order with respect to the y-axis is chosen in the same way as in Eq. (4) with the corresponding $\Delta _y$ test. The resulting integer order of $\overline {N_{x,y,s,int}}$ with two indices is used in the bilinear interpolation. The factor $N_{y,tot}$ is multiplied additionally to $N_{x,tot}$ onto the returned spectrum in order to account for the uniform distribution along which $n_{y,rand}$ has been chosen.

3. Results and discussion

For testing and verification purposes of the presented unification concept, the widely used Cornell Box is appropriate. Its idea has been first presented by the Computer Science department of the name giving Cornell University [26]. The verification box in this work has been built from simeropor plates, where several color spraying layers maintain a homogeneous color and smooth surfaces of the adjoining walls, as seen in Fig. 3.

 

Fig. 3. The photograph of the used Cornell Box.

Download Full Size | PDF

An array of LEDs with a wide emission angle and a nearly constant emission rate over the visible wavelength spectrum is used as the upper light source. The choice for an appropriate diffraction element is rather difficult. In a first comparison, a periodic diffractive element should be placed within the Cornell Box. The periodic cells need to be defined in the RCWA simulation. However, for commercial diffractive optical elements there is not enough information present about the periodic unit cell to be able to define the cell accurately in the RCWA simulation. For this reason, a CD has been used for testing and verification purposes, since there is a lot more information available about the explicit structuring of the CD’s surface. At this point, the model assumes that the CD’s data cells may be represented by a single periodic cell. Furthermore, conventional refractive glass elements are positioned in the box to induce more reflection and transmission artifacts for a precise verification.

Figure 3 shows the photograph of the Cornell Box and Fig. 4 the corresponding rendering results. The pictures have been calculated with 2045 samples per pixel at a resolution of 978 to 734 pixels. The BWDF-method uses $N_{x,tot} = 11$ diffractive orders while 31 input wavelength samples evenly distributed between 400nm and 700nm have been simulated in the RCWA. In addition, the parameters $\Delta _x$ and $\Delta _y$ used for the case comparisons in the following examples are set to 0.1. The graph on the right of Fig. 1 shows the periodic cell approximation of the CD’s surface used. The refractive indices of $n_A = 1.0$ for air, $n_P = 1.5$ for polycarbonate and $n_M$ for aluminum are used herein where $n_M$ is wavelength dependent and complex. Furthermore, $n_A$ and $n_M$ surround the given cell from above and below respectively. The angle discretization step in the RCWA simulation is $2^{\circ }$ for the polar and azimuthal angle. Concerning the azimuthal angle, the structure’s symmetry is utilized and consequently, only the input azimuthal angle range from $0^{\circ }$ to $90^{\circ }$ is simulated in the RCWA package. The RCWA-simulation itself is set up with 61 diffractive orders and its entire simulation time is up to four hours on a desktop personal computer with Intel processor I7-3770k (16 GB RAM, 4 CPU kernels). The rendering process itself, using the precomputed RCWA-values in a 1.2GB array stored on RAM during runtime, requires 22 minutes on this PC architecture.

 

Fig. 4. The rendering result with the BWDF function for a 1D-grating shows a suitable approximation of the light’s interaction onto the diffractive optical element, herein the CD. The rendered picture on the left uses pbrt’s “direct lighting”-integrator, the picture on the right pbrt’s “path”-integrator.

Download Full Size | PDF

For the purpose of comparison, pbrt’s “direct lighting”-integrator, where the direction of $\omega _i$ of a continuous BxDF material is dependent on the returned luminance weight, and pbrt’s “path”-integrator without such spatial dependence are used. Furthermore, the presented BDWF-method, as an extension of pbrt, uses pbrt’s texture architecture, which, in the case of the CD, is used for the transition into the polar coordiante system of the CD. The texture concept of pbrt allows the user to define an arbitrary rotation mapping within the input script.

It is evident that the color gradients within the use of both integrators are of identical shape and position compared to the photograph. Even the color gradient in the top left corner is visible in the rendered image created with the “path”-integrator. Therefore, both results substantiate the consistency between the rendered images and the real-world photograph.

However, the effects of the different integrators are visible as well. The “direct lighting”-integrator causes the dull ceiling to appear rather dark since it is only illuminated from ambient light and the corresponding ray paths are significantly less frequently chosen due to the already mentioned weight dependence. This effect is not existent with the use of the “path”-integrator, however this integrator induces more noise on all surfaces. It is generally presumed that the integration of the introduced BWDF as a further random process would induce a higher degree of noise in the rendered image. Nevertheless, comparing the rendering results, the degree of noise in the shadow next to the glass cube is not lower than the noise on the CD’s surface. This proves that the presented BWDF method doesn’t induce more noise than any other material sampling function of pbrt. Furthermore, the inclusion of the CD’s RCWA sampling function doesn’t require significantly more computation time than for e.g. using the sampling for a mirror or glass material in place of the CD’s surface.

Figure 5 demonstrates the comparison between photograph and rendering in a second setup und view. The rendered image uses the same amount of samples and pixels as in the previous case of Fig. 4. The presented result strengthens the converging approximation of the visible wave-optical effects since the color artifacts and positions of the rendered image converge towards the photograph. However, this setup shows a few problems concerning the exact definition of the scenery in the rendering algorithm. At first, the rendered image depicts a generally higher luminance of the scene as can be seen by a comparison of the sidewall’s brightness. Secondly, the photograph shows only one irradiated artifact on the CD in the lower center of the image, whereas the rendered result depicts an additional color gradient in the same region. This discrepancy is presumably caused by the camera, of which not all necessary properties are known up to this point. It is not clear, in which way the used camera handles very bright light incidence considering any scaling or clamping. Furthermore, the photograph depicts a light green artifact on the lower right corner of the CD, which is not existent in the corresponding rendered image. A possible explanation for this discrepancy is that there is a slight deviation in camera orientation or field of view between photograph and rendered image which results in the fact that the corresponding diffraction order is not visible in the rendering.

 

Fig. 5. The Cornell Box photograph is shown on the left, the rendering result with the BWDF function of the CD is shown on the right. The rendered result uses pbrts “path”-integrator.

Download Full Size | PDF

Up to this point, only a purely reflective 1D-periodic surface has been considered. As described in section 2.2 the BWDF concept can be extended to 2D periodic structures as well as to reflecting and transmitting components from both sides. An example of this “see through” principle is displayed in Fig. 6. The shown structure is a purely imaginary model. Each periodic cell with period $\Lambda _x = 1\mathrm{\mu}$m and $\Lambda _y = 1\mathrm{\mu}$m respectively and thickness 1$\mathrm{\mu}$m consists of a centrally positioned gold plate with feed size 0.5$\mathrm{\mu}$m $\times$ 0.5$\mathrm{\mu}$m $\times$ 1$\mathrm{\mu}$m surrounded by air. The rendered picture uses the same sampling, discretization and pixel modalities as in Fig. 4. Concerning the azimuthal angle, the structure’s symmetry is utilized and consequently, only the input azimuthal angle range from $0^{\circ }$ to $45^{\circ }$ is simulated in the RCWA package. The gold cells are herein arranged laterrally along the shown grating field. It shall be singled out that the gold gradient of the grating is purely the result of the RCWA field simulation instead of defining an object’s color explicitly before runtime. This shows the capabilities of including a full-field approach solving Maxwell’s equations.

 

Fig. 6. The rendering result with the inclusion of the “see through” 2D-BWDF function. The “path”-integrator of pbrt is used.

Download Full Size | PDF

4. Summary and outlook

A method for unifying spectral rendering and field simulations for complex sceneries is presented in this work. Compared to previous works, this method uses RCWA solutions to Maxwell’s equations in a backward rendering and path tracing approach, which proves, as demonstrated on the example of a CD, a high degree of conformity in the presented comparisons of photograph and rendered result. This makes the presented BWDF method suitable for optics design including the physically correct rendering of complex periodic structures with arbitrary cell design and materials. The BWDF concept induces a more accurate representation of the diffractive element and the light interaction with its surrounding scenery giving a real-life and see through image of the element. Its reduced computational effort due to the pre-computational approach in combination with the state of the art renderer pbrt suffices as a useful tool for determining the diffractive optical elements performance in its later application as well as experiencing the real life aesthetics of the diffractive optical element as a major design criterion.

Further work will contain the fine tuning of the camera parameters and the investigated diffractive element (see Section 3). Optimizations with respect to the BWDF’s speedup as well as replacing the $\Delta _x$ and $\Delta _y$ test with an appropriate interpolation function are in progress. Furthermore, the inclusion of inhomogeneous luminance, which cannot be approximated by a single homogeneous wave, of only partially or non periodic structures is contemplated. With this modified field approach metalenses [27] and other diffractive structures or even entire optical systems will be traceable for design, verification and testing purposes.