Freeform optics, which refers to refractive or reflective optical elements with freeform surfaces, is considered a revolution in modern optics^[1]. Freeform surfaces are those surfaces that do not have rotational or translational symmetries^[2]. Compared with traditional optical surfaces, including spherical, aspherical, and cylindrical surfaces, freeform surfaces offer much greater design freedom, allowing optical systems to achieve previously unimaginable imaging^[3–5], illumination^[6,7], and laser beam shaping performances^[8]. The performance evaluation of these freeform optical systems is highly dependent on the light transport simulation from source to target. An efficient light transport computation is crucial in finding the optimal freeform optical systems through optimization techniques.

Ray tracing is currently the most popular and powerful technique for light transport simulation. The main operation of ray tracing is to acquire the intersection points between light rays and surfaces. The computation efficiency of this process is strongly influenced by the surface complexities. Freeform surfaces for imaging can be conveniently described by polynomials, e.g.,  $XY$ polynomials^[3]. However, plenty of illumination problems, e.g.,  picture-generating freeform lens design^[6], are solved numerically. The resulting freeform surfaces, which contain many local curvature features, are hard to describe by using polynomial expansions. Time and memory consumption of ray tracing could increase rapidly as the surface representation becomes more and more complex. Generally, iterative methods have to be employed to obtain the intersection points, which heavily rely on the initial solutions. The meshing method is commonly used to find appropriate initial solutions^[9–11]. In this method, the grids are checked one by one to judge if the intersection occurs on the current grid. However, this method becomes less efficient when the required number of mesh grids is large.

Optical surfaces that are highly freeform and contain both global and local structures can be represented by B-splines, especially by its subset, non-uniform rational B-splines (NURBS)^[12,13]. However, it is not easy to determine the intersection points, even their initial approximations, between the rays and the NURBS surface with high efficiency. The clipping method can be employed to accelerate the ray-tracing process of NURBS surfaces^[14,15]. This method adaptively divides the surface parameters into uniform parts on the 2D parameter space to find the region where the light ray intersects with the surface. Then, the region is gradually subdivided to determine an appropriate initial solution. However, this method suffers from reporting wrong intersections because the sampling points on the surface are not uniformly distributed when parameters are uniformly sampled^[15,16].

Aside from the intersection problem, another important issue is the way of implementing ray tracing. Although some deterministic methods have been discussed, it is restricted in many applications^[17]. In the evaluation and optimization process, the Monte Carlo (MC) ray-tracing method gains wider acceptance as the gold standard for calculating light transport through optical systems^[18–20]. MC ray tracing randomly samples light rays on a light source with a specific probability and traces the light path through the optical system by Snell’s laws. Then, it accumulates the light energy transmitted to the receiver to calculate the irradiance distribution.

Some commercial softwares can perform MC ray tracing of freeform optical elements and systems, but the related algorithms are rarely discussed publicly. Many self-developed ray-tracing programs are successfully applied to artificial image simulation^[21–27], aberration analysis^[28–32], and optical inspection^[33–36]. Some entire ray-tracing pipelines for nonimaging systems have been described in detail^[11,20]. Most of the current MC ray-tracing algorithms trace rays one by one, based on the same law and formula.

The ray-tracing process could be accelerated if plenty of rays are traced simultaneously. However, it is a high-dimensional problem that the parallel intersection process requires simultaneous acquisition of the initial solutions of all current rays. The traditional way of obtaining an initial solution to the intersection problem of a single ray and the freeform surface requires simultaneous computation on mesh grids, which is not suitable for parallel ray tracing. In addition, when the initial solution is not accurate enough, the parameters of the intersection could easily escape from the defined range as the iteration progresses, resulting in many calculation errors.

We implement MC parallel ray tracing (PRT) for freeform illumination lenses with NURBS surfaces to realize a fast irradiance evaluation. We achieve uniform sampling by inverse transform sampling and fast surface interpolation by adopting the analytical form of the B-spline basis function. We propose a parameter transform method to reduce the accuracy requirement of the initial solution during the intersection process using Newton’s method, which could improve the percentage of successfully traced light rays by parallel computation. This fast MC simulation program can facilitate the evaluation and optimization of freeform lenses. Section 2 describes in detail the parallel ray-tracing method for freeform NURBS surfaces. Section 3 gives two examples to demonstrate the effectiveness of the proposed method.

The purpose of the ray-tracing process is to trace the changes in the spatial parameters $(x,y)$ and orientation parameters $(r_x,r_y)$, which uniquely characterize a light ray in a three-dimensional space. $(x,y,r_x,r_y)$ form a four-dimensional parameter space, which degenerates to a two-dimensional one for a special case when all the rays are emitted from a point light source. Figure 1 illustrates the procedures of our parallel ray-tracing algorithm. After defining the optical system parameters, we sample the light source randomly in both spatial locations and orientations, generating a random set of rays. We then implement MC ray propagation from the source, through the freeform lenses, to the receiver. The difficulty arises in simultaneously determining all the parameters $(x,y,r_x,r_y)$ for the current set of rays transporting through a complex freeform surface. In this decisive and important step, we propose a parameter transformation method to reduce the ray-surface intersection errors in the iterative process of Newton’s method. Finally, we count the number of rays that reach different positions of the receiver to determine the irradiance distribution. The three steps are described in detail in the following section.

The light source is required to be sampled randomly in both spatial locations and orientations. Different sampling probability density functions $P(x,y,r_x,r_y)$ should ensure that $w(x,y,r_x,r_y)·P(x,y,r_x,r_y)$ is proportional to $L(x,y,r_x,r_y)·\cos \theta$, where $L$ is the radiance distribution of the source, $w$ denotes the weights of the sampled rays, and $\theta$ denotes the angles from the source surface normal to the sampled rays. The initial values of $w$ can be set as 1, and then $P$ is proportional to the intensity distribution of the light source. However, this can add computation complexity when $L$ is not a constant.

Here, we first sample the position coordinates and solid angle of the light source uniformly and then assign the initial weights proportional to $L·\cos \theta$. In this way, there is an equivalent relation between the radiance distribution $L$ and the weights of the sampled light rays, $$L(x,y,r_x,r_y)·\cos \theta \delta S\delta \mathrm{\Omega }\propto E_{\mu }\left(\sum _{i=0}^k{w}_i\right).$$(1)

Herein, $\delta S$ denotes a small area on the surface around $(x,y)$, and $\delta \mathrm{\Omega }$ denotes a small solid angle around $(x,y,r_x,r_y)$. $E_{\mu }({\sum }_{i=0}^k{w}_i)$ is the mathematical expectation of the summed weight of the light rays inside the small space $\delta S\delta \mathrm{\Omega }$, and $k$ denotes the number of rays inside the small space.

In a spherical coordinate system, $(r_x,r_y)$ can be replaced with the zenith angle $\theta$ and the azimuth angle $\phi$, as shown in Fig. 2(a). Note that the light rays sampled uniformly with $\theta$ and $\phi$ are not uniformly distributed in space^[19]. Since the magnitude of the small solid angle $\delta \mathrm{\Omega }$ around $(\theta ,\phi )$ is proportional to $\sin \theta \delta \theta \delta \phi$, uniform sampling of the solid angle should guarantee that the probability is proportional to $\phi$ and $\sin \theta$.

We employ the inverse transform sampling method^[37] to obtain random variables with arbitrary probability distribution $f(x)$, where $x\in [0,1]$. This method first calculates the cumulative probability distribution function $F(x)={\int }_{-\infty }^{\infty }f(x)\mathrm{d}x$ according to the required probability distribution function $f(x)$. A uniformly distributed random variable $Y=F(x)$ can be generated. Finally, the random variable $X$ with probability distribution $f(X)$ can be obtained as $X=F^{-1}(x)$. For the sampling of the zenith angle here, we first obtain a random uniformly distributed variable $Y$ between 0 and 1; then the probability distribution function $X=\arccos (1-Y)$ is a sine function. Thus, uniform sampling of random rays in space can be generated based on setting the probability distribution functions as $\theta \sim \sin \theta$, $\phi \sim U(0,2\pi )$.

Figures 2(b) and 2(c) illustrate the ray distributions of a Lambertian light source sampled in two strategies: $\theta \sim \sin \theta$ and $\theta \sim U(\mathrm{0,0.5}\pi )$, where $\phi$ is uniformly sampled for each strategy. The pseudo-code of the random rays sampling process is shown in Algorithm 1.

After generating a sequence of random rays from the light source, we trace the rays propagating through the freeform optical surfaces represented with NURBS. The key operation is the determination of the intersection points between rays and surfaces. Thus, we accelerate the intersection process by modifying the surface interpolation strategy and introducing a parameter transformation method. After obtaining the intersection points, we use Snell’s law in vector form to calculate the outgoing ray vectors.

An NURBS surface is obtained as the tensor product of two NURBS curves parameterized by two parameters $u$ and $v$, $$\mathbf{Q}(u,v)=\sum _{i=0}^{m_s-1}\sum _{j=0}^{n_s-1}{N}_{i,p}(u)N_{j,q}(v)R_{i,j}\mathbf{P}(i,j),$$(2)where $u$ and $v$ lie in $[0,1]$, $\mathbf{P}(i,j)$ denotes the control point, $R_{i,j}$ is the weight of $\mathbf{P}(i,j)$, and $N_{i,p}(u)$ and $N_{j,q}(v)$ are the basis functions of the B-splines of degree $p$ and $q$ in the $u$ and $v$ directions, respectively. The basis function $N_{i,p}(u)$ is defined by the De Boor Cox recursion formula^[12], $$\{\begin{array}{l}{N}_{i,0}=\{\begin{array}{ll}1\hfill & \text{if}u\in [u_i,u_{i+1})\hfill \\ 0\hfill & \text{if}\mathrm{else}\hfill \end{array},\hfill \\ N_{i,p}=\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(u),\hfill \\ \text{define:}\frac{0}{0}=0.\hfill \end{array}$$(3)

Herein, $u_i$ is the selected node which divides $[0,1]$ into $m_s+p$ segments. We employ the quasi-uniform nodes here, which can simplify the formula for the B-spline basis functions.

The computation volume required for interpolation increases exponentially with the number of control points and the interpolation degree. The major reason is that the recursive method repeatedly computes the formula with a large storage of accumulative data and no zero $N(u_i)$ only when $u_i\in [u_i,u_{i+p+1}]$. These unnecessary calculations slow down the intersection point-finding process. Abert et al. transferred the computationally expensive recursion into single instruction multiple-data (SIMD) suitable form to reduce the time cost of executed commands^[13]. Jester et al. proposed the B-spline quasi-interpolation scheme, which improves the computational speed of B-spline surfaces within a certain error range^[38]. To improve the speed of the interpolation process, we employ the analytical form of the B-spline basis function of degree 2^[13]. In addition, only the control points whose corresponding basis functions are not 0 are considered in the interpolation process to avoid meaningless operations. Such a strategy can promote calculation efficiency greatly when the control point number is much larger than the interpolation order.

After specifying the interpolation strategy, let us now turn to the calculation of the intersection between the rays and the NURBS surface. Suppose we have a ray parameterized as $(x(t),y(t),z(t))=\widehat{\mathbf{s}}+\widehat{\mathbf{r}}·t$, where $\widehat{\mathbf{s}}=(x_0,y_0,z_0)$ is the starting point, and $\widehat{\mathbf{r}}=(r_x,r_y,r_z)$ denotes the ray direction. The intersection can be formulated by the following equation: $$\mathbf{Q}(u,v)=\widehat{\mathbf{s}}+\widehat{\mathbf{r}}·t,$$(4)where $\mathbf{Q}=(Q_1,Q_2,Q_3)$ denotes a point on the NURBS surface. After expressing the ray direction with $\theta$ and $\phi$, Eq. (4) can be rewritten as $$\{\begin{array}{l}{Q}_1(u,v)-x_0=\sin \theta \cos \phi ·t\hfill \\ Q_2(u,v)-y_0=\sin \theta \sin \phi ·t\hfill \\ Q_3(u,v)-z_0=\cos \theta ·t\hfill \end{array}.$$(5)

We rotate $\phi$ and $\theta$ around the $z$-axis and $y$-axis, respectively, so that the new $z$-axis coincides with the direction of the light ray. This way can simplify the intersection-solving process^[39]. The transformation matrices describing the rotations are $${\mathbf{R}}_z=\left[\begin{array}{ccc}\cos \phi & \sin \phi & 0\\ -\sin \phi & \cos \phi & 0\\ 0& 0& 1\end{array}\right],{\mathbf{R}}_y=\left[\begin{array}{ccc}\cos \theta & 0& -\sin \theta \\ 0& 1& 0\\ \sin \theta & 0& \cos \theta \end{array}\right],$$(6)where ${\mathbf{R}}_z$ describes the rotation around the $z$-axis, and ${\mathbf{R}}_y$ describes the rotation around the $y$-axis.

After applying the rotation matrices to Eq. (5), we obtain $${\mathbf{R}}_y·{\mathbf{R}}_z·\left[\begin{array}{c}{Q}_x(u,v)-x_0\\ Q_y(u,v)-y_0\\ Q_z(u,v)-z_0\end{array}\right]={\mathbf{R}}_y·{\mathbf{R}}_z·\left[\begin{array}{c}\sin \theta ·\cos \phi ·t\\ \sin \theta ·\sin \phi ·t\\ \cos \theta ·t\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ t\end{array}\right].$$(7)

We take the first two equations, which could be represented as $f_1(u,v)=0\text{and}{f}_2(u,v)=0$. Then, the parameters $(u,v)$ of the intersection point on the surface can be solved by using Newton’s method. Newton’s iteration step is given by $$\left[\begin{array}{c}{u}_{k+1}\\ v_{k+1}\end{array}\right]=\left[\begin{array}{c}{u}_k\\ v_k\end{array}\right]-{\mathbf{J}}^{-1}\left[\begin{array}{c}{f}_1(u_k,v_k)\\ f_2(u_k,v_k)\end{array}\right],$$(8)where $(u_k,v_k)$ is the approximate solution of the $k$th iteration, and J is the Jacobian matrix, $$\mathbf{J}=\left[\begin{array}{cc}\frac{\partial f_1}{\partial u}& \frac{\partial f_1}{\partial v}\\ \frac{\partial f_2}{\partial u}& \frac{\partial f_2}{\partial v}\end{array}\right].$$(9)

Newton’s method needs a good initial estimate. Most of the methods divide the surface into many triangular meshes to find a suitable initial solution. If the meshing is sufficiently detailed, then the initial solution could lead to a good iteration result. Such a strategy is not suitable for parallel computation here.

A slight deviation from the exact initial solution could cause the value of $u$ or $v$ to exceed the range of [0,1], resulting in surface interpolation errors. We adopt a parameter transformation method to solve this problem. In this method, a new pair of parameters, $(\alpha ,\beta )$, are employed to replace $(u,v)$, and their relationships can be described as $$u=\frac{1}{e^{-4\alpha }+1},v=\frac{1}{e^{-4\beta }+1}.$$(10)

The parameter transformation has the following characteristics: $$\{\begin{array}{l}\alpha ,\beta \in (-\infty ,+\infty ),u,v\in (0,1)\hfill \\ \frac{\mathrm{d}u}{\mathrm{d}\alpha }>0,\frac{\mathrm{d}v}{\mathrm{d}\beta }>0,{\frac{\mathrm{d}u}{\mathrm{d}\alpha }|}_{\alpha =0}=1,{\frac{\mathrm{d}v}{\mathrm{d}\beta }|}_{\beta =0}=1\hfill \\ {u|}_{\alpha =0}=\mathrm{0.5,}{v|}_{\beta =0}=0.5\hfill \end{array}.$$(11)

Such a transformation can avoid the values of $(u,v)$ exceeding [0, 1]. The intersection equations now become $f_1(\alpha ,\beta )=0\text{and}{f}_2(\alpha ,\beta )=0$. The new Jacobian matrix J in Newton’s method can be acquired based on the chain rule $$\{\begin{array}{ccc}\frac{\partial f}{\partial \alpha }& =& \frac{\partial f}{\partial u}·\frac{\mathrm{d}u}{\mathrm{d}\alpha }\\ \frac{\partial f}{\partial \beta }& =& \frac{\partial f}{\partial v}·\frac{\mathrm{d}v}{\mathrm{d}\beta }\end{array}.$$(12)

Herein, $f$ can be $f_1$ or $f_2$. We obtain the initial solution using a coarse mesh on the surface. For the $k$th iteration, $({\alpha }_k,{\beta }_k)$ can be calculated from $(u_k,v_k)$ by the inverse operation of Eq. (10). The Newton iteration formula Eq. (8) is changed into $$\left[\begin{array}{c}{\alpha }_{k+1}\\ {\beta }_{k+1}\end{array}\right]=\left[\begin{array}{c}{\alpha }_k\\ {\beta }_k\end{array}\right]-{\mathbf{J}}^{-1}\left[\begin{array}{c}{f}_1({\alpha }_k,{\beta }_k)\\ f_2({\alpha }_k,{\beta }_k)\end{array}\right].$$(13)

Such a parameter transformation could reduce the accuracy requirement of the initial solution to the iterative computation. Thus, we can adopt surface points with a relatively coarse density as the initial solution, which can save a lot of time.

Once we obtain the intersection points between the rays and the surface, we determine the refraction directions of the outgoing rays based on Snell’s law. As illustrated in Fig. 3, $\widehat{\mathbf{i}}$ denotes the unit incident ray vectors, $\widehat{\mathbf{t}}$ denotes the unit outgoing ray vectors, $\widehat{\mathbf{n}}$ denotes the unit normal vectors to the intersection points on the surface, and $n_1$ and $n_2$ denote the refractive indices of the incident and exit media. Then, $\widehat{\mathbf{t}}$ can be obtained based on Snell’s law in vector form, $$\widehat{\mathbf{t}}=\frac{1}{n_2}\{n_1\widehat{\mathbf{i}}+[\sqrt{n_2^2-n_1^2+n_1^2{(\widehat{\mathbf{n}}·\widehat{\mathbf{i}})}^2}-n_1\widehat{\mathbf{n}}·\widehat{\mathbf{i}}]\widehat{\mathbf{n}}\}.$$(14)

Fresnel loss occurs when a light ray passes through a refractive optical surface. The reflectivities $r_s$ and $r_p$ for s-polarized light and p-polarized light, respectively, are written as $$r_s=\frac{n_1\cos {\theta }_i-n_2\cos {\theta }_t}{n_1\cos {\theta }_i+n_2\cos {\theta }_t},r_p=\frac{n_1\cos {\theta }_t-n_2\cos {\theta }_i}{n_1\cos {\theta }_t+n_2\cos {\theta }_i},$$(15)where ${\theta }_i$ and ${\theta }_t$ denote the incident and refracted angles, respectively. The reflectances $R_s$ and $R_p$ for s-polarized light and p-polarized light, respectively, can be obtained as $R_s=r_s^2\text{\hspace{0.17em}and\hspace{0.17em}}{R}_p=r_p^2$. The transmittances $T_s$ and $T_p$ for s-polarized light and p-polarized light, respectively, can be simply obtained as $T_s=1-R_s$ and $T_p=1-R_p$. For natural light or typical LED light, the total transmittance $T$ can be taken as the average of $T_s$ and $T_p$. We adjust the weights $w(x,y,r_x,r_y)$ of the light rays according to the transmittance $T$ during ray tracing.

The pseudo-code for simulating the light transport through a freeform surface is provided in Algorithm 2. Herein, we consider ray tracing as successful when the cosine distance $d=1-\mathrm{dot}(\mathbf{Q}-\widehat{\mathbf{s}},\widehat{\mathbf{r}})/{\Vert \mathbf{Q}-\widehat{\mathbf{s}}\Vert }_2$ is less than the allowed error $\epsilon$, where Q is the found intersection point.

After finishing the ray tracing through the entrance surface, the ray tracing through the second surface can be implemented based on the same procedure, where the intersection points on the entrance surface are set as the starting points of the next ray tracing. Such a process is repeated until the light rays hit the receiver.

After finishing the ray tracing through all the freeform surfaces and determining the intersection points between the rays and the target plane, we count the number of rays for each small grid of the discretized receiver.

Suppose that the total radiant power emitted from the light source is $W_0$. The weight of the $i$th ray is denoted as $w_{0,i}$, and the number of rays is $n$. The radiant power of a ray with the unit weight can be acquired by $$W_{\mathrm{unit}}=\frac{W_0}{\sum _{i=1}^n{w}_{0,i}}.$$(16)

Suppose that the range of the receiver is $\{(x,y)|x_{\min }<x<x_{\max },y_{\min }<y<y_{\max }\}$, which is discretized into $m_c\times n_c$ cells. The area of each cell can be obtained as $$S_{\mathrm{cell}}=\frac{(x_{\max }-x_{\min })·(y_{\max }-y_{\min })}{m_c·n_c}.$$(17)

The energy weight of the $i$th light ray that hits the receiver is changed into $w_{i^{\prime }}$. The total number of rays of the $j$th cell is $n_j$. The irradiance value of the $j$th cell, $j=1,2,\dots$, $m_c\times n_c$, can be acquired by $$E_j=\frac{W_{\mathrm{unit}}\times \sum _{i=1}^{n_j}{w}_{i^{\prime }}}{S_{\mathrm{cell}}}.$$(18)

The pseudo-code of the energy statistic process and the arithmetic flow of our PRT are shown in Algorithm 3. Note that the more the total traced rays, the higher the signal-to-noise ratio (SNR) of the simulated irradiance distribution. A discussion of the appropriate number of rays to be traced for a certain setting can be found in Ref. [40].

To demonstrate the effectiveness of our method, we implement irradiance evaluation of the picture-generating freeform lenses designed with the iterative wavefront tailoring (IWT) method^[41]. As shown in Fig. 4, the light rays emitted from a light source propagate through the entrance and exit optical surfaces and then generate a complex irradiance distribution of the picture type at the receiver. The entrance surface can be a spherical, aspherical, or freeform surface, and the exit surface is a freeform one. Note that the freeform surface profiles of the picture-generating freeform lenses can be very complex, which is difficult to be described by conventional polynomials, e.g.,  $XY$ or Zernike polynomials.

The ray tracing performances are evaluated by the time consumption, the intersection precision, and the success rate. The intersection precision is measured by the maximum or average value of the Euclidean distances from the calculated intersection points on the exit surface to the light rays. The success rate $\tau$ is defined as the ratio of the number of successfully traced rays to the total number of traced rays.

Our computer is equipped with an Intel i9-12900 k CPU with x86-64 architecture and an NVIDIA RTX3090 GPU with 24 GB video memory and NVIDIA Ampere architecture. The programming language is Python 3.11. We use the GPU to trace batches of light rays simultaneously and the CPU to control the calculation flow. We use the CuPy package of Python for surface modeling, intersection solving, and other parallel scientific calculations. The computing time is mainly affected by the GPU performance.

The first example is the ray tracing of a spherical-freeform lens as shown in Fig. 5(a). The refractive index of the lens is 1.4932. This lens is aimed at converting the light distribution of an LED source into a $2000\mathrm{mm}\times 2000\mathrm{mm}$ irradiance distribution at $z=1000\mathrm{mm}$, forming the letters “BIT”. The maximum dimensions of the lens along the $x$, $y$, and $z$ directions are 16.97 mm, 16.28 mm, and 12.01 mm, respectively. It can be seen from Fig. 5(b) that the freeform exit surface has complex local features which are directly related to the target irradiance pattern. The entrance surface is $50\times 50$ 3D grids sampled from a semi-sphere of 5 mm radius, which is also expressed with NURBS. The complex exit freeform surface is described by $256\times 256$ control points. We employ a $3\times 3$ uniform filter for the simulated irradiance distribution to decrease the noise of the MC method. Figure 6 shows the simulated irradiance distributions of our algorithm under different settings of the number of light rays $n$, the number of receiver cells $m=m_c\times n_c$, and the source sizes. Table 1 gives the corresponding time consumption, precision of the solved intersections, and the success rates.

As shown in Table 1, for each case of the first example, the success rate $\tau$ is higher than 99.9%. $\tau$ could be reduced to below about 96% without using the parameter transformation method. In addition, the parameter transformation method improves the convergence of Newton’s method for solving the intersection points. The iterative solutions would not jump out of the allowed range of the parameters in the early stage of the iterative process. However, the intersection calculation errors can rise rapidly toward the edge of the surface ($u,v\to 0$ or 1) since the derivatives $(\mathrm{d}\alpha /\mathrm{d}u,\mathrm{d}\beta /\mathrm{d}v)$ of the transformation functions tend to 0. This problem could be avoided by extending the surface so that the active area is bounded away from the edge of the surface.

The second example is about the ray tracing of a more complex double-freeform lens with the maximum dimensions of $15.33\mathrm{mm}\times 16.47\mathrm{mm}\times 10.04\mathrm{mm}$ [see Fig. 7(a)], whose refractive index is also 1.4932. Each of the freeform surfaces contains $1024\times 1024$ control points. This freeform lens aims at generating a portrait of Johann Carl Friedrich Gauss^[42] with the size of $240\mathrm{mm}\times 240\mathrm{mm}$ at $z=100\mathrm{mm}$ from a point-like source whose size is $0.001\mathrm{mm}\times 0.001\mathrm{mm}$. We perform the PRT using ${10}^7$ rays. The simulated irradiance distribution is shown in Fig. 7(b). The maximum and average values of the distances from the calculated intersection points on the exit surface to the light rays are $2.33\times {10}^{-7}\mathrm{mm}$ and $5.16\times {10}^{-8}\mathrm{mm}$, respectively. Our program only takes $\sim 11.52\mathrm{s}$ to perform the simulation, and $\tau$ is about 99.951%.

We analyze the time complexity of our PRT program for the first example with respect to the number of rays, the number of receiver cells, the size of the light source, and the number of control points. We use the control variable method to practically test the time consumption of the program. The results are presented in Fig. 8. The time complexity for the number of rays is $O(n)$. Because $n$ is much larger than the number of GPU cores, parallel computing is still done in batches on the hardware. As shown in Fig. 8(a), the simulation time increases linearly as $n$ increases. For the number of receiver cells, the total time complexity is approximately $O(1)$. As shown in Fig. 8(b), as $m$ increases, the simulation time remains almost unchanged. The reason is that most of the time is spent on the intersection-seeking process. Theoretically, the time complexity of the ray energy statistical process is $O(m)$. However, since $n$ is typically much larger than $m$ to ensure low-level Monte Carlo noise, the impact of $m$ on the overall simulation time is negligible. The total time complexity for the source size is about $O(1)$, as illustrated in Fig. 8(c). The source size only affects the range of the starting points of the randomly sampled rays. We adopt the 2-order analytic form of the NURBS basis functions, which can avoid many recursive operations in base function modeling. Thus, the total time complexity for the control points of the lens surfaces is $O(1)$, as demonstrated in Fig. 8(d).

We have proposed a PRT method for freeform lenses with complex NURBS surfaces. In this method, we employ the inverse transform sampling strategy to uniformly sample rays in space. We also adopt the analytical form of the B-spline basis function of 2 degrees to achieve a fast surface interpolation. More importantly, we proposed a parameter transform method to reduce the initial solution requirement of Newton’s method for solving intersections, which thus increases the success rate of parallel ray tracing. Two examples verify that the proposed PRT method is efficient and effective for quick irradiance evaluations of complex freeform irradiance tailoring lenses.

The proposed PRT method is also applicable to freeform surfaces with other expressions, including $XY$ polynomials, Zernike polynomials, and radial base functions. Our method can help designers quickly evaluate their design results, especially for freeform surfaces with high degrees of freedom. As a relatively fast irradiance evaluation technique, our method can also be applied to the optimal design of freeform lenses. The proposed method is currently restricted to those cases in which a single light ray does not have multiple intersections with a surface. However, such a restriction is fulfilled by most of the freeform lens designs for illumination and beam-shaping applications. Future work may include the development of non-sequential PRT algorithms and the application of freeform illumination lens design using the automatic differentiation technique in the framework of MindSpore.