SimOptDevice: a library for virtual optical experiments

Virtual experiments have become an indispensable tool for the design and the accuracy assessment of novel measurement procedures and instruments. Virtual experiments are particularly relevant in modern optics due to its challenging demands for highly accurate measurements. This paper introduces SimOptDevice, a flexible library for opto-mechanical virtual experiments. After describing the scope and general structure of the library, its underlying mathematical tools used for solving the related numerical tasks are described. Finally, the application of SimOptDevice to a recent interferometric measurement procedure is presented.

1 Introduction

Following the advance of technology and the demand for highly accurate measurements, optical instruments and experiments have become very complex in recent years. In addition, sophisticated data analysis has become an important part of modern optical measurement devices. To ensure that a measurement principle is fit for its purpose, it is beneficial to first test it in a virtual environment prior to building the physical setup. That way, experimenters can save time and costs in the development of novel procedures. Furthermore, virtual experiments are often essential for the assessment of accuracies that can be reached.

For these reasons, virtual experiments have become an important tool in optics. Examples of applications are non-null interferometer calibration (Hao et al., 2016), validation of new data analysis techniques (Shen et al., 2015), absolute flatness measurements of optical surfaces (Bouillet and Morin, 2014), and accuracy evaluations in interferometric measurements (Wiegmann et al., 2011) or deflectometric flatness measurements (Schulz et al., 2010). Virtual experiments have become essential also in many other scientific fields, e.g. simulations of X-ray optics experiments (Knudsen et al., 2013), neutron scattering (Lieutenant et al., 2004), uncertainty assessment in computer tomography (Hiller and Reindl, 2012) or coordinate measurement machines (Heißelmann et al., 2017; Trenk et al., 2004), cross-borehole imaging (Donato and Crocco, 2015), or error quantification of CNC milling machines (Soori et al., 2013).

The Physikalisch-Technische Bundesanstalt (PTB) has developed the SimOptDevice software library for optical virtual experiments. SimOptDevice is a flexible library implemented in MATLAB^® (MATLAB, 2018) that covers a large range of applications. It facilitates the design of experiments and allows one to develop and test measurement procedures. Furthermore, SimOptDevice is used to optimize existing measurement procedures with respect to measurement time and uncertainty. The latter is particularly relevant for PTB, which aims at measurements at the highest level of accuracy. Successful examples using SimOptDevice include accuracy evaluation of interferometric measurements of a synchrotron mirror (Wiegmann et al., 2011), development of a deflectometric flatness reference at PTB (Schulz et al., 2010; Ehret et al., 2012), and accuracy tests for multi-spectral imaging systems (Dierl et al., 2018).

In this paper we will describe the structure of SimOptDevice and its underlying mathematical methods. We will address the key issues of optical virtual experiments and refer to the corresponding solutions implemented in SimOptDevice. The purpose of this paper is to describe the mathematical methods necessary for implementing a complex optical simulation environment, enabling readers to set up a software solution of their own.

The general structure of SimOptDevice and its basic principles are described in Sect. 2. In Sect. 3 we introduce the mathematical methods used for the solution of related numerical problems. A prominent example is then presented in Sect. 4. Finally, some conclusions are given.

2 Scope and general description

The SimOptDevice library can be used to perform virtual optical measurements with complex beam paths running through a series of optical elements. It considers nested scanning stages accounting for translations and rotations, and supports the use of various sensors such as cameras. SimOptDevice is based on the application of ray optics.

The basic principle of SimOptDevice is a system of hierarchical coordinate systems, combined with ray tracing routines. Within each coordinate system, optical elements can be placed. A related local coordinate system is assigned to each considered optical element, along with a superordinate coordinate system relative to which the local coordinate system is defined. In this way, a tree structure of coordinate systems is built. Each local coordinate system can undergo individual rotations and translations with respect to its superordinate system. The coordinates of each element can be transformed into any of the other employed coordinate systems. Those transformations are made simple by using homogeneous coordinates which are introduced in Sect. 3.1.

The power of SimOptDevice lies in tracing rays and performing ray aiming accurately and efficiently according to the laws of refraction and reflection while being in control of all optimization parameters and the applied algorithms. Using the library for our experiments, we view and verify all intermediate results, which is very helpful for tuning the algorithms for each specific problem. Whereas ray tracing follows a ray through the optical system when start point and direction are given, ray aiming seeks the path through the system for given start and end points. The latter is a highly nonlinear optimization problem. Details are described in Sect. 3.

The accurate modelling of all elements of a measurement setup is another advantage of SimOptDevice. This includes not only optical elements like lenses, mirrors and sensors but also linear stages and rotary tables. Ensembles of elements can be saved and reused in other virtual experiments.

During development, the software has been successfully checked against ray tracing results obtained by ZEMAX. This included comparisons of optical path lengths and of points reached by ray tracing. The differences were in the sub-nanometre range.

3 Mathematical methods

3.1 Coordinate transformations

One key feature of SimOptDevice is its easy way to transform coordinates from any local system to any other coordinate system within the defined structure (see Fig. 1). Coordinates are defined as homogeneous coordinates (see e.g. Cox et al., 2007, p. 357 ff.). Rotations and translations can be represented as matrices, and it is simple to connect transformations and compute a transformation's inverse. Homogeneous coordinates for describing points in ℝ³ have four components. A vector v=(x, y, z)^T is represented in homogeneous coordinates as ${\mathit{v}}^{\prime }=(x^{\prime }$, $y^{\prime }{z}^{\prime }$, w)^T, where w≠0 is arbitrary and $x^{\prime }=x/w$, $y^{\prime }=y/w$, and $z^{\prime }=z/w$. Usually, w is set to one.

Figure 1 Example of the hierarchical structure of coordinate systems in SimOptDevice. (a) The plot of a coordinate system tree example is shown. The left column lists the element names; the right column shows a sketch of the hierarchical structure of the coordinate systems. (b) The corresponding instrument is illustrated. The position and orientation of each of the instrument's subsystems are defined with respect to the superordinate system through a transformation of homogeneous coordinates. A transformation matrix H can be computed for transformations from one system to another (see Eq. 2). It is a function of the source system and the destination system (in this example Sensor and Topo, respectively). The common superordinate system Ts_Frame is needed for the computation of the transformation matrix (cf. Eq. 2). Subsequently, coordinates of a point in one system can be transformed to any other system by multiplication by the corresponding transformation matrix H.

With this coordinate definition, a translation or rotation is represented by a 4×4 matrix:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{Translation}:\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{\Delta }x\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{\Delta }y\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{\Delta }z\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right),\\ {\displaystyle }& {\displaystyle }\mathrm{Rotation}\mathrm{about}x\mathrm{axis}:\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \cos \mathit{\theta }& -\sin \mathit{\theta }& \mathrm{0}\\ \mathrm{0}& \sin \mathit{\theta }& \cos \mathit{\theta }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right).\end{array}$$

Rotations about the other axes are defined analogously. The composition of several homogeneous transformations H₁, H₂, …, H_N is equal to a single homogeneous transformation

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathbf{H}={\mathbf{H}}_N\cdot {\mathbf{H}}_{N-\mathrm{1}}\cdot \mathrm{\dots }\cdot {\mathbf{H}}_{\mathrm{2}}\cdot {\mathbf{H}}_{\mathrm{1}}.\end{array}$$

Using Eq. (1), the translation and rotation of a coordinate system with respect to its superordinate system are represented by a single matrix. The inverse of a homogeneous transformation is represented by the inverse of the corresponding transformation matrix. The aforementioned transformation from one local source system S_a to another destination system S_b is performed with the help of a common superordinate system S_c:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{\mathbf{H}}_{a\to b}={\mathbf{H}}_{b\to c}^{-\mathrm{1}}\cdot {\mathbf{H}}_{a\to c}.\end{array}$$

For a more comprehensive introduction to homogeneous coordinates, see e.g. Cox et al. (2007) or Hartley and Zisserman (2003).

3.2 Ray tracing

Ray propagation through the defined object is performed by sequential ray tracing. For this tracing method the order of elements passed by a ray is known in advance and optical paths are computed element by element. SimOptDevice computes the ray paths locally in each optical element system. At the boundary T_i between two elements (e.g. a surface of a lens), the current intersection point p_i of the ray with that boundary along with the new ray direction e_i are calculated in the new local system (cf. Sect. 3.1). The ray is then traced until it intersects with the next boundary, and so on. The intersection point with T_i and the new direction are a function of the previous intersection point and ray direction:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{\left({\mathit{p}}_i^T,{\mathit{e}}_i^T\right)}^T=f_i\left({\mathit{p}}_{i-\mathrm{1}},{\mathit{e}}_{i-\mathrm{1}}\right).\end{array}$$

Function f_i entails performing the following steps (cf. Fig. 2 for an illustration):

1. transformation of p_i−1 and e_i−1 from system T_i−1 to T_i,

2. calculation of the ray's geometrical path length l_i between T_i−1 and T_i,

3. determination of next intersection point ${\mathit{p}}_i={\mathit{p}}_{i-\mathrm{1}}+l_i\cdot {\mathit{e}}_{i-\mathrm{1}}$,

4. calculation of normal vector n_i on T_i in p_i, and

5. calculation of the new ray direction e_i according to the laws of refraction and reflection.

Figure 2 Schematics of a ray tracing step. Each topography T has its own coordinate system. The intersection point p_i with topography T_i is computed from the previous intersection point p_i−1 and the previous ray direction e_i−1. Subsequently, the new ray direction e_i at T_i is computed according to the laws of reflection and refraction.

This is classical ray tracing with the particular feature that at each boundary the intersection point and the new direction are calculated in the new local coordinate system. Step 3 is calculated analytically for planes and spherical surfaces and has to be calculated numerically for more complex surfaces like Zernike surfaces, aspheres, or surfaces described by a Gauss function. If the ray passes N surfaces, there are N different functions f_i. The last intersection point and ray direction can be expressed as a function of the start point and direction by concatenation of functions f₁ to f_N:

$$\begin{array}{ll}{\displaystyle }{\left({\mathit{p}}_N^T,{\mathit{e}}_N^T\right)}^T& {\displaystyle }=f_N\left(f_{N-\mathrm{1}}\left(\mathrm{\dots }{f}_{\mathrm{2}}\left(f_{\mathrm{1}}\left({\mathit{p}}_{\mathrm{0}},{\mathit{e}}_{\mathrm{0}}\right)\right)\mathrm{\dots }\right)\right)\\ \text{(4)}& {\displaystyle }& {\displaystyle }=f\left({\mathit{p}}_{\mathrm{0}},{\mathit{e}}_{\mathrm{0}}\right).\end{array}$$

Furthermore, for each crossing of a boundary T_i a Jacobian matrix containing the partial derivatives of the intersection point p_i and the ray direction e_i with respect to the previous intersection point p_i−1 and the previous ray direction e_i−1 is calculated. It is associated with Eq. (3) and has the following form:

$${\displaystyle }{\displaystyle }{\mathbf{J}}_i=\left(\begin{array}{cccc}{\displaystyle \frac{\partial x_i}{\partial x_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial x_i}{\partial y_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial x_i}{\partial {\widehat{x}}_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial x_i}{\partial {\widehat{y}}_{i-\mathrm{1}}}}\\ {\displaystyle \frac{\partial y_i}{\partial x_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial y_i}{\partial y_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial y_i}{\partial {\widehat{x}}_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial y_i}{\partial {\widehat{y}}_{i-\mathrm{1}}}}\\ {\displaystyle \frac{\partial {\widehat{x}}_i}{\partial x_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial {\widehat{x}}_i}{\partial y_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial {\widehat{x}}_i}{\partial {\widehat{x}}_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial {\widehat{x}}_i}{\partial {\widehat{y}}_{i-\mathrm{1}}}}\\ {\displaystyle \frac{\partial {\widehat{y}}_i}{\partial x_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial {\widehat{y}}_i}{\partial y_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial {\widehat{y}}_i}{\partial {\widehat{x}}_{i-\mathrm{1}}}}& {\displaystyle \frac{\partial {\widehat{y}}_i}{\partial {\widehat{y}}_{i-\mathrm{1}}}}\end{array}\right),$$

where x_i and y_i are the x and y components of the point p_i, and ${\widehat{x}}_i$ and ${\widehat{y}}_i$ are the x and y components of the direction e_i. We only need the derivatives with respect to the x and y components, since the intersected topographies are known and the z component is fixed at the intersection point. Since the direction vector has unit length, its z component can be calculated by $\widehat{z}=\sqrt{\mathrm{1}-{\widehat{x}}^{\mathrm{2}}-{\widehat{y}}^{\mathrm{2}}}$. The Jacobian matrix J corresponding to the concatenation of N−1 interfaces i=1 … N−1 associated with Eq. (4) is obtained by multiplication of the corresponding matrices:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathbf{J}={\mathbf{J}}_{N-\mathrm{1}}\cdot {\mathbf{J}}_{N-\mathrm{2}}\mathrm{\cdots }{\mathbf{J}}_{\mathrm{2}}\cdot {\mathbf{J}}_{\mathrm{1}}.\end{array}$$

The inverse Jacobian matrix J⁻¹ describes partial derivatives of the ray's behaviour along the same path in the opposite direction. These Jacobian matrices are used extensively when performing the ray aiming. Their usage is explained in more detail in Sect. 3.3. A detailed description of the ray tracing and ray aiming procedures can also be found in Fortmeier (2016, in German).

3.3 Ray aiming

Given a start point p₀ and an end point ${\mathit{p}}_N^{\mathrm{dest}}$ the task of ray aiming is to find the start direction e₀ from p₀ in order to minimize the distance between p_N and the desired end point ${\mathit{p}}_N^{\mathrm{dest}}$ (see Fig. 3):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{\widehat{\mathit{e}}}_{\mathrm{0}}=\arg \underset{{\mathit{e}}_{\mathrm{0}}}{\min }\left|\left|{\mathit{p}}_N\left({\mathit{e}}_{\mathrm{0}}\right)-{\mathit{p}}_N^{\mathrm{dest}}\right|\right|.\end{array}$$

where $\left|\right|\cdot \left|\right|$ denotes the Euclidean norm and p_N is reached through application of ray tracing for a chosen start direction e₀. For successful ray aiming the resulting norm in Eq. (6) is close to zero. This is usually a highly nonlinear problem. It can be solved using one of MATLAB^®'s parallel nonlinear solver routines, e.g. lsqnonlin. The latter is a solver for nonlinear least-squares problems in which the trust-region-reflective or Levenberg–Marquard algorithm can be applied. For better convergence of the solver, the ray tracing Jacobian matrix from Eq. (5) is utilized in this step. The Jacobian matrices calculated during ray tracing are used to update the start direction e₀ when minimizing the distance between point p_N and point ${\mathit{p}}_N^{\mathrm{dest}}$. Furthermore, SimOptDevice delivers the Jacobian matrices for the change in total optical path lengths with respect to changes in a point p_i or a ray direction e_i for each surface T_i (Fortmeier et al., 2014). It is calculated analytically, thereby omitting the computationally expensive extra ray tracing steps for the numerical differentiation. The Jacobian matrices for ray aiming are important in many interferometric applications, e.g. the tilted-wave interferometer (cf. Sect. 4), where the change in the optical path length of a ray is a significant quantity of the experiment.

Figure 3 Ray aiming principle: the start and end points p₀ and ${\mathit{p}}_N^{\mathrm{dest}}$ of a ray are given. The task is to find the correct start angle e₀ such that the desired end point is hit. This is achieved by minimizing the distance between points p_N and ${\mathit{p}}_N^{\mathrm{dest}}$, when p_N is reached through application of ray tracing for a chosen start direction, e₀.

A requirement for successful ray aiming is that the destination point ${\mathit{p}}_N^{\mathrm{dest}}$ can be reached. Therefore, the valid area on the CCD (or any other final surface) is determined in a preceding step. A ray aiming between the light source and some characteristic points at the smallest aperture in the optical system is performed. The final ray directions at the aperture are used to trace the rays to the CCD, thereby defining the area that can be reached by rays from the sources. A ray aiming is regarded as successful if the norm of the difference in Eq. (6) is smaller than 1 nm. Local minima that prevent the algorithm from converging are detected and such rays are masked. Using the Jacobian matrix, the ray aiming usually converges in about four steps. Without the Jacobian matrix this takes approximately 12 to 15 steps, depending on the particular situation.

Figure 4 Examples from the tilted-wave interferometer (TWI) simulation. (a) Ray paths through the instrument. Source array, reference and test arms, and detector (CCD). (b) The sources of the laser array are switched on in four groups (S1 to S4); (c) simulated interferograms on the CCD for each of the source groups for an asphere example.

4 Application example: tilted-wave interferometer

The tilted-wave interferometer (TWI) is an interferometric measurement device for form measurements of asphere (Garbusi et al., 2008) and freeform surfaces. Figure 4 illustrates the simulation of interferograms for a measurement of an asphere. For SimOptDevice, the TWI is a particular application example since virtual experiments are not only used for designing the instrument and measurement setup, but are also an integrated part of data analysis. More precisely, we use SimOptDevice to obtain a numerical model of the experiment in dependence on the unknown form of the asphere or freeform surface under test. The latter is retrieved as the solution of a nonlinear inverse problem where in our implementation SimOptDevice constitutes the model. The inverse problem is solved iteratively, which is facilitated through the parallel computing capabilities of SimOptDevice. The nonlinear problem is linearized with the help of Jacobian matrices (Fortmeier et al., 2014). The solution of the inverse problem is found with the MATLAB^® routine fsolve. It is possible to use different solver algorithms with fsolve. We typically choose trust-region-dogleg or trust-region-reflective. Trust region algorithms locally approximate the cost function by a quadratic function and progress towards the minimum of this approximation.

In order to accomplish an accurate recovery of the specimen surface, it is necessary to obtain an accurate model of the TWI (Garbusi et al., 2008). Differences between the numerical model of the instrument and the actual instrument have to be characterized very accurately prior to a measurement. For this calibration process, well-known test specimens (typically spheres) are measured with the TWI (Baer et al., 2014). The numerical model is then updated by adjusting the Zernike polynomial parametrization of two reference surfaces in order to account for the remaining discrepancies.

SimOptDevice has also been used for uncertainty evaluation (Fortmeier et al., 2017) of the TWI and for the exploration of new measurement concepts with the TWI (Fortmeier et al., 2016).

5 Conclusions

SimOptDevice is a versatile library for conducting virtual opto-mechanical experiments that has been applied successfully in several projects and studies. SimOptDevice can model a large number of optical elements and sensors which can be combined flexibly to cover a wide range of experimental setups.

We have explained the mathematical concepts within the library. A detailed description of our ray tracing and ray aiming procedures and of the determination of Jacobian matrices, needed for efficiently solving the nonlinear inverse problems, was given. This will be useful for readers interested in implementing virtual optical experiments. We also conclude that SimOptDevice can be used to simulate very complex opto-mechanical systems.