5012
OPTICS LETTERS / Vol. 38, No. 23 / December 1, 2013
Target flux estimation by calculating intersections between neighboring conic reflector patches Cristina Canavesi,1,* William J. Cassarly,2 and Jannick P. Rolland1 1
The Institute of Optics, University of Rochester, 275 Hutchison Rd, Rochester, New York 14627, USA 2 Synopsys, Inc., 1100 Hunt Club Drive, Wooster, Ohio 44691, USA *Corresponding author: [email protected] Received August 1, 2013; revised October 15, 2013; accepted October 24, 2013; posted October 25, 2013 (Doc. ID 194936); published November 21, 2013
We propose a fast algorithm to estimate the flux collected by conic reflector patches, based on the calculation of intersections between neighboring patches. The algorithm can be employed in conjunction with the supporting ellipsoids algorithm for freeform reflector design and is shown to be orders of magnitude faster and more scalable than the commonly used Monte Carlo ray tracing approach. © 2013 Optical Society of America OCIS codes: (080.4298) Nonimaging optics; (080.4228) Nonspherical mirror surfaces; (220.2945) Illumination design; (080.1753) Computation methods. http://dx.doi.org/10.1364/OL.38.005012
Reflectors are commonly used in a variety of illumination applications, including street lighting, automotive headlamps, light projection engines, and medical lighting. In general, if the setup does not have rotational or translational symmetry, a freeform optical surface governed by a nonlinear second-order partial differential equation of the Monge–Ampère is required. No analytic calculation is available to describe the surface, and standard numerical integration techniques often prove to be unstable. Usually a point source problem is solved first, from which the surface for the finite source is then derived. Several iterative design methods have been proposed to derive the reflective surface to direct light from the given point source to the desired target, which is discretized in a finite number of target points or directions for the near- and far-field problems, respectively. Ries and Muschaweck use a multigrid approach, but few details are provided on its actual implementation [1]. Two other methods, the linear programming and supporting ellipsoids algorithms, use the reflective property of conics by building the reflector from patches of conics, each of which reflects a set of rays from the source to a different target point. All the conic patches share a common focus at the origin, where the point source is located. The linear programming approach, proposed independently by Wang and co-workers [2,3], addresses the far-field problem with a variational approach using paraboloids as the conics that make up the reflector but is severely limited by computational complexity as the resolution of the target discretization increases [4,5]. Additionally, it was shown that the variational approach cannot be used for the near-field problem [6]. The Oliker supporting ellipsoids method is versatile and can be used to produce arbitrary target distributions in either the near or far field [7]. A major limitation of the supporting ellipsoid method is the time needed for the flux evaluation, which is usually done by Monte Carlo ray tracing [8]. The motivation for the research presented here is to improve the computation time of the supporting ellipsoids method by reducing the time needed to perform the flux evaluation. In the following, we refer to the near-field problem, but the results can easily be extended to the far-field case 0146-9592/13/235012-04$15.00/0
with no additional complexity by replacing the expression of an ellipsoid with that of a paraboloid. In the first step of the supporting ellipsoids method as described in [9], the ellipsoids are initialized so that all the flux is collected by one ellipsoid, called the reference ellipsoid; when all ellipsoids collect some flux, they are said to be supporting. The ellipsoids are iteratively scaled, and the target flux is evaluated to determine whether the exit condition is met. The merit function is usually defined as the squared sum of the differences between desired and actual normalized targets. Two solutions exist for any given problem: convergent or divergent geometry [8]. A convergent reflector is one for which the ellipsoids reflect rays from the point source toward the opposite side of the target with respect to the source, and rays reflected from different targets cross; a divergent reflector directs rays on the same side of the target, and rays from different ellipsoids do not cross. The selection between the convergent or divergent case must be made before the initialization of the ellipsoids. The target evaluation is repeated several times in each iteration cycle of the supporting ellipsoids algorithm. Typically, Monte Carlo ray tracing is used; rays are generated randomly according to the source distribution, and the number of rays that each ellipsoid collects is counted [10]. Assuming equal weighting for all rays, the flux collected by an ellipsoid is proportional to the number of rays it collects. In the general brute force approach, rays are traced from the point source to all ellipsoids; the distance to each ellipsoid is evaluated, and the ellipsoid for which the distance is smaller (for the convergent case) or greater (for the divergent case) is considered to collect that ray. The variance decreases linearly with the number of samples; thus the error of the estimate is given by the inverse of the square root of the number of rays. To halve the error in the estimation, four times as many rays have to be traced [10]. Let us consider a target made up of N x × N y points. To evaluate the target distribution with an accuracy of 1∕sqrtn, it is necessary to trace nN x N y rays (n rays per ellipsoid). Since all rays are traced to all ellipsoids in the brute force approach, a total of nN x N y N x N y ray tracing operations are needed to evaluate the target © 2013 Optical Society of America
December 1, 2013 / Vol. 38, No. 23 / OPTICS LETTERS
distribution. The computation time increases linearly with the number of rays per ellipsoid, n, and quadratically with the total number of ellipsoids, N x N y . The flux collected by each ellipsoid patch is proportional to its area; if the perimeter of each patch is identified, the area can be calculated and used to estimate the collected flux [11]. Finding the perimeter of each conic patch consists of finding the vertices between three or more patches, and the boundary curves between pairs of conics, as shown in Fig. 1 for a reflector made of 4 × 4 ellipsoids. The problem of finding intersections between quadric surfaces has been explored for several decades in the field of computer graphics; recently Dupont et al. made a significant contribution to the field by developing a method to calculate the intersection curves between two generic quadrics [12]. This method can be applied to the supporting ellipsoids problem; however, since it deals with quadric surfaces arbitrarily located in space, the framework used is overly complicated for the specific case considered in the supporting ellipsoids method— that of a collection of ellipsoids sharing a common focus at the origin—resulting in an unfavorable computation time. Thus, we derive an intersection calculation method that has a significantly faster computation than the Dupont method for the specific case considered. We express a prolate ellipsoid of revolution (prolate spheroid) in its polar form as ρ
f ; 1 − eS · T
(1)
where ρ is the distance to the ellipsoid along a given ray direction, f and e are the focal parameter and ellipticity of the ellipsoid, and S and T are unit vectors that describe the source ray direction and target point, respectively. One focus of the ellipsoid is at the origin, and the second focus is at the target point. Since all ellipsoids share a common focus at the origin, the problem of finding the intersections between two ellipsoids is equivalent to finding the values of S for which Eq. (1) produces the same value of ρ. In general, the intersection occurs along a curve that in the geometry
Fig. 1. Reflector surface obtained as a solution of the supporting ellipsoids method with flux evaluation performed by Monte Carlo ray tracing for a 4 × 4 discrete target, as rendered in LightTools. The vertices and boundaries of one ellipsoid patch are highlighted.
5013
produced within the supporting ellipsoids method has the shape of an ellipse. Kochengin and Oliker have in fact shown how the intersection between two ellipsoids with a common focus forms a circular disk when projected on the unit sphere [13]. Conversely, the intersection between three ellipsoids occurs not along a curve, but at individual points in space. Finding the intersection point between three ellipsoids or the intersection curve between two ellipsoids that fall within the desired portion of the reflector is a mechanical exercise of algebra; an example of the intersections calculated for two and three ellipsoids is shown in Fig. 2. Once the intersection curves and points that describe the border and vertices of an ellipsoid are known, the area of the ellipsoid patch can be calculated in cosine space to provide an estimation of the collected flux [11]. The computation can be further simplified by only calculating the three-ellipsoid intersection points to estimate the area of each ellipsoid patch. The intersection curve describing the boundary between two neighboring ellipsoids is an ellipse, but the portion of the ellipse that makes up the border of interest is nearly linear, and can be approximated with a straight line. Thus, once the intersection points between three neighboring ellipsoids are known, it is unnecessary to also calculate the intersection curves between pairs of neighboring ellipsoids; the points can be connected with straight lines to represent the boundary, and the area in the projected x–y plane is calculated as the area of a polygon. Let us consider a target made of a grid of N x × N y points. For each quadruplet of ellipsoids, four intersection points between three ellipsoids can be calculated. Thus, the total number of three-ellipsoid intersection operations to be performed is 4N x − 1N y − 1 at most. The underlying assumption made to calculate the intersections is that the ellipsoids that make up the reflector are arranged in a manner that closely mirrors the disposition of the target points; thus, neighboring reflector patches direct light to neighboring target points. With this assumption, it is possible to identify the triplets of ellipsoids that intersect at each vertex point. Similarly, it is possible to derive the intersection curve that constitutes the boundary between two ellipsoids. The flux estimation with Monte Carlo ray tracing and with the three-ellipsoid intersection calculation was implemented in Visual Basic on an Intel i-7
Fig. 2. (a) Two-ellipsoid and (b) three-ellipsoid intersections calculated with the intersection method. The intersection curves (ellipses) between two ellipsoids are shown with a solid line; the dots in (b) represent the two intersection points that occur in space.
5014
OPTICS LETTERS / Vol. 38, No. 23 / December 1, 2013
2.7 GHz processor in a Windows 7 virtual machine with 2 GB of RAM and two processor cores. A single Monte Carlo ray tracing operation (one ray and one ellipsoid) is approximately 15 times faster than a single threeellipsoid intersection calculation. However, since the number of Monte Carlo ray tracing calculations (with a chosen n number of rays per ellipsoid) needed to evaluate the target distribution, nN x N y N x N y , is significantly higher than the 4N x − 1N y − 1 intersection calculations required, and it increases quadratically with the number of targets, the intersection method is proven to be several orders of magnitude faster, and more scalable than the Monte Carlo approach for all cases of practical interest (n greater than 1, and N x , N y greater than 4). The relative computation time versus the number of target points is shown in Fig. 3 using the intersection method and Monte Carlo ray tracing (with 10 or 100 rays per ellipsoid), normalized to the computation time for the intersection method with 1024 target points. The intersection method is based on the assumption that the arrangement of neighboring ellipsoids is known; a further assumption made is that the ellipsoids that are along the border of the reflector are also known. As a result, the intersection method cannot be used from the beginning of the supporting ellipsoid algorithm described in [7]; instead, we start using Monte Carlo ray tracing to evaluate the target distribution, and when the assumption is valid we switch to flux estimation by the intersection method. As an example, let us consider a uniform square target of width 10 located at z 100, discretized in a grid of 10 × 10 points, and a point source placed at the origin emitting a Lambertian distribution into a 40° square in the tangent plane. The reflector obtained after 100 iterations of the supporting ellipsoids method run with target distribution evaluation by Monte Carlo ray tracing with 100 rays per ellipsoid is shown in Fig. 4(a). Figure 4 depicts the projection in the x–y plane of the reflector surface; an average of 100 Monte Carlo rays per ellipsoid are depicted, with different colors representing different ellipsoids, and the boundaries of each ellipsoid calculated with the intersection method are also shown. Since in this example the supporting ellipsoids method was stopped before converging to the solution after 100 iterations, the ellipsoid at the bottom left (the reference ellipsoid) notably collects more flux than all the other
Fig. 3. Relative flux estimation time for the intersection method (cross) and Monte Carlo ray tracing with 10 rays/ ellipsoid (circle) or 100 rays/ellipsoid (square), normalized to the computation time for 1024 target points with the intersection method.
Fig. 4. Projection in the x–y plane of the rays hitting the 10 × 10 reflector obtained with the supporting ellipsoids algorithm (a) with Monte Carlo flux estimation for 100 iterations, and (b) switching to the intersection method until the final solution is obtained. Rays hitting different ellipsoids have different colors for visualization purposes; the vertices and boundaries calculated with the intersection method are shown in black.
ellipsoids, as expected. Once the reflector of Fig. 4(a) was obtained, the supporting ellipsoids method was continued using the intersection method for the target distribution evaluation. The final reflector solution obtained is shown in Fig. 4(b). The estimated flux for the converged solution of Fig. 4(b) is shown in Fig. 5. The gray scale image depicts 10 × 10 bins with values equal to the normalized flux calculated with Monte Carlo ray tracing and the intersection method. For the example considered (100 target points and 1000 rays per ellipsoid), 10 million Monte Carlo ray tracing operations are needed to perform one target evaluation, compared to 324 three-ellipsoid calculations with the intersection method. The Monte Carlo ray trace, done with 1000 rays per ellipsoid, took 365 ms, while the intersection method took 0.795 ms (459 times faster than the Monte Carlo ray trace). The statistical error produced using 1000 rays in the Monte Carlo ray tracing is approximately 3%; on the other hand, the area estimation obtained with the intersection method has a normalized standard deviation below 0.05%. In order to obtain such a level of accuracy with the Monte Carlo approach, 4 million rays per ellipsoid would have to be traced; extrapolating the computation time from the values
Fig. 5. Map of the normalized flux collected by the 10 × 10 reflector of Fig. 4(b) as calculated by (a) Monte Carlo ray tracing with 1000 rays/ellipsoid and (b) intersection method. The scale of the images is set to display the entire range of values of the Monte Carlo ray tracing result.
December 1, 2013 / Vol. 38, No. 23 / OPTICS LETTERS
Fig. 6. (a) Projection in the x–y plane of a reflector that produces a target with 45° twist, and (b) normalized flux estimated with the intersection method (the scale for the visualization is set at 5%).
obtained with 1000 rays per ellipsoid, tracing four million rays per ellipsoid would take ∼14; 600 s, or in excess of 4 h. Figure 6 shows an example of a reflector that produces a square uniform 10 × 10 target distribution rotated 45° about the z axis with 0.13% normalized standard deviation; the intersection method is still capable of correctly estimating the flux even when the patches differ considerably in shape. We have implemented a flux estimation method for conic reflector patches based on calculation of the intersection points between triplets of ellipsoids. A major advantage of the intersection method is that, unlike Monte Carlo ray tracing, it does not face a trade-off between accuracy and computation speed. Additionally, the complexity of the intersection method increases only linearly with the resolution of the discrete target, thus proving to be more scalable than the quadratic dependence on the number of target points obtained with Monte Carlo ray tracing. The intersection method can be used in conjunction with the supporting ellipsoids algorithm to produce a reflector for a prescribed target illumination. In the current implementation, due to the assumptions made regarding the a priori knowledge of neighboring ellipsoids from the target point distribution, the intersection
5015
method still requires that a brute force Monte Carlo ray tracing be done for a few iterations of the supporting ellipsoids algorithm, until a stable starting point is achieved. The assumptions made in the current implementation of the intersection method for flux estimation can be relaxed to address a wider range of cases. Future work will address the applicability of the flux estimation method to various target configurations of interest, quantification of its impact on the convergence of the supporting ellipsoids method, and possible use in conjunction with other reflector design methods. The authors acknowledge Synopsys, Inc. for the educational license of LightTools. This work was funded under a fellowship from Synopsys, Inc., the National Science Foundation (IIP-1236999, EECS-1002179, IIP1338877), and the NYSTAR Foundation (C050070). References 1. H. Ries and J. Muschaweck, J. Opt. Soc. Am. A 19, 590 (2002). 2. X.-J. Wang, Calc. Var. Partial Differential Equation 20, 329 (2004). 3. T. Glimm and V. Oliker, J. Math. Sci. 117, 4096 (2003). 4. V. Oliker, Proc. SPIE 5942, 594207 (2005). 5. C. Canavesi, W. J. Cassarly, and J. P. Rolland, Opt. Express 20, 4050 (2012). 6. T. Graf and V. I. Oliker, Inverse Probl. 28, 025001 (2012). 7. V. I. Oliker, in Trends in Nonlinear Analysis, M. Kirkilionis, S. Kromker, R. Rannacher, and F. Tomi, eds. (SpringerVerlag, 2003), p. 191. 8. F. R. Fournier, W. J. Cassarly, and J. P. Rolland, Proc. SPIE 7423, 742302 (2009). 9. S. A. Kochengin and V. O. Oliker, Comput. Visualization Sci. 6, 15 (2003). 10. P. Shirley and R. K. Morley, Realistic Ray Tracing, 2nd ed. (A K Peters, 2003). 11. D. G. Koch, Proc. SPIE 1780, 226 (1992). 12. L. Dupont, D. Lazard, S. Lazard, and S. Petitjean, J. Symb. Comput. 43, 168 (2008). 13. S. A. Kochengin and V. I. Oliker, Numer. Math. 79, 553 (1998).
OPTICS LETTERS / Vol. 38, No. 23 / December 1, 2013
Target flux estimation by calculating intersections between neighboring conic reflector patches Cristina Canavesi,1,* William J. Cassarly,2 and Jannick P. Rolland1 1
The Institute of Optics, University of Rochester, 275 Hutchison Rd, Rochester, New York 14627, USA 2 Synopsys, Inc., 1100 Hunt Club Drive, Wooster, Ohio 44691, USA *Corresponding author: [email protected] Received August 1, 2013; revised October 15, 2013; accepted October 24, 2013; posted October 25, 2013 (Doc. ID 194936); published November 21, 2013
We propose a fast algorithm to estimate the flux collected by conic reflector patches, based on the calculation of intersections between neighboring patches. The algorithm can be employed in conjunction with the supporting ellipsoids algorithm for freeform reflector design and is shown to be orders of magnitude faster and more scalable than the commonly used Monte Carlo ray tracing approach. © 2013 Optical Society of America OCIS codes: (080.4298) Nonimaging optics; (080.4228) Nonspherical mirror surfaces; (220.2945) Illumination design; (080.1753) Computation methods. http://dx.doi.org/10.1364/OL.38.005012
Reflectors are commonly used in a variety of illumination applications, including street lighting, automotive headlamps, light projection engines, and medical lighting. In general, if the setup does not have rotational or translational symmetry, a freeform optical surface governed by a nonlinear second-order partial differential equation of the Monge–Ampère is required. No analytic calculation is available to describe the surface, and standard numerical integration techniques often prove to be unstable. Usually a point source problem is solved first, from which the surface for the finite source is then derived. Several iterative design methods have been proposed to derive the reflective surface to direct light from the given point source to the desired target, which is discretized in a finite number of target points or directions for the near- and far-field problems, respectively. Ries and Muschaweck use a multigrid approach, but few details are provided on its actual implementation [1]. Two other methods, the linear programming and supporting ellipsoids algorithms, use the reflective property of conics by building the reflector from patches of conics, each of which reflects a set of rays from the source to a different target point. All the conic patches share a common focus at the origin, where the point source is located. The linear programming approach, proposed independently by Wang and co-workers [2,3], addresses the far-field problem with a variational approach using paraboloids as the conics that make up the reflector but is severely limited by computational complexity as the resolution of the target discretization increases [4,5]. Additionally, it was shown that the variational approach cannot be used for the near-field problem [6]. The Oliker supporting ellipsoids method is versatile and can be used to produce arbitrary target distributions in either the near or far field [7]. A major limitation of the supporting ellipsoid method is the time needed for the flux evaluation, which is usually done by Monte Carlo ray tracing [8]. The motivation for the research presented here is to improve the computation time of the supporting ellipsoids method by reducing the time needed to perform the flux evaluation. In the following, we refer to the near-field problem, but the results can easily be extended to the far-field case 0146-9592/13/235012-04$15.00/0
with no additional complexity by replacing the expression of an ellipsoid with that of a paraboloid. In the first step of the supporting ellipsoids method as described in [9], the ellipsoids are initialized so that all the flux is collected by one ellipsoid, called the reference ellipsoid; when all ellipsoids collect some flux, they are said to be supporting. The ellipsoids are iteratively scaled, and the target flux is evaluated to determine whether the exit condition is met. The merit function is usually defined as the squared sum of the differences between desired and actual normalized targets. Two solutions exist for any given problem: convergent or divergent geometry [8]. A convergent reflector is one for which the ellipsoids reflect rays from the point source toward the opposite side of the target with respect to the source, and rays reflected from different targets cross; a divergent reflector directs rays on the same side of the target, and rays from different ellipsoids do not cross. The selection between the convergent or divergent case must be made before the initialization of the ellipsoids. The target evaluation is repeated several times in each iteration cycle of the supporting ellipsoids algorithm. Typically, Monte Carlo ray tracing is used; rays are generated randomly according to the source distribution, and the number of rays that each ellipsoid collects is counted [10]. Assuming equal weighting for all rays, the flux collected by an ellipsoid is proportional to the number of rays it collects. In the general brute force approach, rays are traced from the point source to all ellipsoids; the distance to each ellipsoid is evaluated, and the ellipsoid for which the distance is smaller (for the convergent case) or greater (for the divergent case) is considered to collect that ray. The variance decreases linearly with the number of samples; thus the error of the estimate is given by the inverse of the square root of the number of rays. To halve the error in the estimation, four times as many rays have to be traced [10]. Let us consider a target made up of N x × N y points. To evaluate the target distribution with an accuracy of 1∕sqrtn, it is necessary to trace nN x N y rays (n rays per ellipsoid). Since all rays are traced to all ellipsoids in the brute force approach, a total of nN x N y N x N y ray tracing operations are needed to evaluate the target © 2013 Optical Society of America
December 1, 2013 / Vol. 38, No. 23 / OPTICS LETTERS
distribution. The computation time increases linearly with the number of rays per ellipsoid, n, and quadratically with the total number of ellipsoids, N x N y . The flux collected by each ellipsoid patch is proportional to its area; if the perimeter of each patch is identified, the area can be calculated and used to estimate the collected flux [11]. Finding the perimeter of each conic patch consists of finding the vertices between three or more patches, and the boundary curves between pairs of conics, as shown in Fig. 1 for a reflector made of 4 × 4 ellipsoids. The problem of finding intersections between quadric surfaces has been explored for several decades in the field of computer graphics; recently Dupont et al. made a significant contribution to the field by developing a method to calculate the intersection curves between two generic quadrics [12]. This method can be applied to the supporting ellipsoids problem; however, since it deals with quadric surfaces arbitrarily located in space, the framework used is overly complicated for the specific case considered in the supporting ellipsoids method— that of a collection of ellipsoids sharing a common focus at the origin—resulting in an unfavorable computation time. Thus, we derive an intersection calculation method that has a significantly faster computation than the Dupont method for the specific case considered. We express a prolate ellipsoid of revolution (prolate spheroid) in its polar form as ρ
f ; 1 − eS · T
(1)
where ρ is the distance to the ellipsoid along a given ray direction, f and e are the focal parameter and ellipticity of the ellipsoid, and S and T are unit vectors that describe the source ray direction and target point, respectively. One focus of the ellipsoid is at the origin, and the second focus is at the target point. Since all ellipsoids share a common focus at the origin, the problem of finding the intersections between two ellipsoids is equivalent to finding the values of S for which Eq. (1) produces the same value of ρ. In general, the intersection occurs along a curve that in the geometry
Fig. 1. Reflector surface obtained as a solution of the supporting ellipsoids method with flux evaluation performed by Monte Carlo ray tracing for a 4 × 4 discrete target, as rendered in LightTools. The vertices and boundaries of one ellipsoid patch are highlighted.
5013
produced within the supporting ellipsoids method has the shape of an ellipse. Kochengin and Oliker have in fact shown how the intersection between two ellipsoids with a common focus forms a circular disk when projected on the unit sphere [13]. Conversely, the intersection between three ellipsoids occurs not along a curve, but at individual points in space. Finding the intersection point between three ellipsoids or the intersection curve between two ellipsoids that fall within the desired portion of the reflector is a mechanical exercise of algebra; an example of the intersections calculated for two and three ellipsoids is shown in Fig. 2. Once the intersection curves and points that describe the border and vertices of an ellipsoid are known, the area of the ellipsoid patch can be calculated in cosine space to provide an estimation of the collected flux [11]. The computation can be further simplified by only calculating the three-ellipsoid intersection points to estimate the area of each ellipsoid patch. The intersection curve describing the boundary between two neighboring ellipsoids is an ellipse, but the portion of the ellipse that makes up the border of interest is nearly linear, and can be approximated with a straight line. Thus, once the intersection points between three neighboring ellipsoids are known, it is unnecessary to also calculate the intersection curves between pairs of neighboring ellipsoids; the points can be connected with straight lines to represent the boundary, and the area in the projected x–y plane is calculated as the area of a polygon. Let us consider a target made of a grid of N x × N y points. For each quadruplet of ellipsoids, four intersection points between three ellipsoids can be calculated. Thus, the total number of three-ellipsoid intersection operations to be performed is 4N x − 1N y − 1 at most. The underlying assumption made to calculate the intersections is that the ellipsoids that make up the reflector are arranged in a manner that closely mirrors the disposition of the target points; thus, neighboring reflector patches direct light to neighboring target points. With this assumption, it is possible to identify the triplets of ellipsoids that intersect at each vertex point. Similarly, it is possible to derive the intersection curve that constitutes the boundary between two ellipsoids. The flux estimation with Monte Carlo ray tracing and with the three-ellipsoid intersection calculation was implemented in Visual Basic on an Intel i-7
Fig. 2. (a) Two-ellipsoid and (b) three-ellipsoid intersections calculated with the intersection method. The intersection curves (ellipses) between two ellipsoids are shown with a solid line; the dots in (b) represent the two intersection points that occur in space.
5014
OPTICS LETTERS / Vol. 38, No. 23 / December 1, 2013
2.7 GHz processor in a Windows 7 virtual machine with 2 GB of RAM and two processor cores. A single Monte Carlo ray tracing operation (one ray and one ellipsoid) is approximately 15 times faster than a single threeellipsoid intersection calculation. However, since the number of Monte Carlo ray tracing calculations (with a chosen n number of rays per ellipsoid) needed to evaluate the target distribution, nN x N y N x N y , is significantly higher than the 4N x − 1N y − 1 intersection calculations required, and it increases quadratically with the number of targets, the intersection method is proven to be several orders of magnitude faster, and more scalable than the Monte Carlo approach for all cases of practical interest (n greater than 1, and N x , N y greater than 4). The relative computation time versus the number of target points is shown in Fig. 3 using the intersection method and Monte Carlo ray tracing (with 10 or 100 rays per ellipsoid), normalized to the computation time for the intersection method with 1024 target points. The intersection method is based on the assumption that the arrangement of neighboring ellipsoids is known; a further assumption made is that the ellipsoids that are along the border of the reflector are also known. As a result, the intersection method cannot be used from the beginning of the supporting ellipsoid algorithm described in [7]; instead, we start using Monte Carlo ray tracing to evaluate the target distribution, and when the assumption is valid we switch to flux estimation by the intersection method. As an example, let us consider a uniform square target of width 10 located at z 100, discretized in a grid of 10 × 10 points, and a point source placed at the origin emitting a Lambertian distribution into a 40° square in the tangent plane. The reflector obtained after 100 iterations of the supporting ellipsoids method run with target distribution evaluation by Monte Carlo ray tracing with 100 rays per ellipsoid is shown in Fig. 4(a). Figure 4 depicts the projection in the x–y plane of the reflector surface; an average of 100 Monte Carlo rays per ellipsoid are depicted, with different colors representing different ellipsoids, and the boundaries of each ellipsoid calculated with the intersection method are also shown. Since in this example the supporting ellipsoids method was stopped before converging to the solution after 100 iterations, the ellipsoid at the bottom left (the reference ellipsoid) notably collects more flux than all the other
Fig. 3. Relative flux estimation time for the intersection method (cross) and Monte Carlo ray tracing with 10 rays/ ellipsoid (circle) or 100 rays/ellipsoid (square), normalized to the computation time for 1024 target points with the intersection method.
Fig. 4. Projection in the x–y plane of the rays hitting the 10 × 10 reflector obtained with the supporting ellipsoids algorithm (a) with Monte Carlo flux estimation for 100 iterations, and (b) switching to the intersection method until the final solution is obtained. Rays hitting different ellipsoids have different colors for visualization purposes; the vertices and boundaries calculated with the intersection method are shown in black.
ellipsoids, as expected. Once the reflector of Fig. 4(a) was obtained, the supporting ellipsoids method was continued using the intersection method for the target distribution evaluation. The final reflector solution obtained is shown in Fig. 4(b). The estimated flux for the converged solution of Fig. 4(b) is shown in Fig. 5. The gray scale image depicts 10 × 10 bins with values equal to the normalized flux calculated with Monte Carlo ray tracing and the intersection method. For the example considered (100 target points and 1000 rays per ellipsoid), 10 million Monte Carlo ray tracing operations are needed to perform one target evaluation, compared to 324 three-ellipsoid calculations with the intersection method. The Monte Carlo ray trace, done with 1000 rays per ellipsoid, took 365 ms, while the intersection method took 0.795 ms (459 times faster than the Monte Carlo ray trace). The statistical error produced using 1000 rays in the Monte Carlo ray tracing is approximately 3%; on the other hand, the area estimation obtained with the intersection method has a normalized standard deviation below 0.05%. In order to obtain such a level of accuracy with the Monte Carlo approach, 4 million rays per ellipsoid would have to be traced; extrapolating the computation time from the values
Fig. 5. Map of the normalized flux collected by the 10 × 10 reflector of Fig. 4(b) as calculated by (a) Monte Carlo ray tracing with 1000 rays/ellipsoid and (b) intersection method. The scale of the images is set to display the entire range of values of the Monte Carlo ray tracing result.
December 1, 2013 / Vol. 38, No. 23 / OPTICS LETTERS
Fig. 6. (a) Projection in the x–y plane of a reflector that produces a target with 45° twist, and (b) normalized flux estimated with the intersection method (the scale for the visualization is set at 5%).
obtained with 1000 rays per ellipsoid, tracing four million rays per ellipsoid would take ∼14; 600 s, or in excess of 4 h. Figure 6 shows an example of a reflector that produces a square uniform 10 × 10 target distribution rotated 45° about the z axis with 0.13% normalized standard deviation; the intersection method is still capable of correctly estimating the flux even when the patches differ considerably in shape. We have implemented a flux estimation method for conic reflector patches based on calculation of the intersection points between triplets of ellipsoids. A major advantage of the intersection method is that, unlike Monte Carlo ray tracing, it does not face a trade-off between accuracy and computation speed. Additionally, the complexity of the intersection method increases only linearly with the resolution of the discrete target, thus proving to be more scalable than the quadratic dependence on the number of target points obtained with Monte Carlo ray tracing. The intersection method can be used in conjunction with the supporting ellipsoids algorithm to produce a reflector for a prescribed target illumination. In the current implementation, due to the assumptions made regarding the a priori knowledge of neighboring ellipsoids from the target point distribution, the intersection
5015
method still requires that a brute force Monte Carlo ray tracing be done for a few iterations of the supporting ellipsoids algorithm, until a stable starting point is achieved. The assumptions made in the current implementation of the intersection method for flux estimation can be relaxed to address a wider range of cases. Future work will address the applicability of the flux estimation method to various target configurations of interest, quantification of its impact on the convergence of the supporting ellipsoids method, and possible use in conjunction with other reflector design methods. The authors acknowledge Synopsys, Inc. for the educational license of LightTools. This work was funded under a fellowship from Synopsys, Inc., the National Science Foundation (IIP-1236999, EECS-1002179, IIP1338877), and the NYSTAR Foundation (C050070). References 1. H. Ries and J. Muschaweck, J. Opt. Soc. Am. A 19, 590 (2002). 2. X.-J. Wang, Calc. Var. Partial Differential Equation 20, 329 (2004). 3. T. Glimm and V. Oliker, J. Math. Sci. 117, 4096 (2003). 4. V. Oliker, Proc. SPIE 5942, 594207 (2005). 5. C. Canavesi, W. J. Cassarly, and J. P. Rolland, Opt. Express 20, 4050 (2012). 6. T. Graf and V. I. Oliker, Inverse Probl. 28, 025001 (2012). 7. V. I. Oliker, in Trends in Nonlinear Analysis, M. Kirkilionis, S. Kromker, R. Rannacher, and F. Tomi, eds. (SpringerVerlag, 2003), p. 191. 8. F. R. Fournier, W. J. Cassarly, and J. P. Rolland, Proc. SPIE 7423, 742302 (2009). 9. S. A. Kochengin and V. O. Oliker, Comput. Visualization Sci. 6, 15 (2003). 10. P. Shirley and R. K. Morley, Realistic Ray Tracing, 2nd ed. (A K Peters, 2003). 11. D. G. Koch, Proc. SPIE 1780, 226 (1992). 12. L. Dupont, D. Lazard, S. Lazard, and S. Petitjean, J. Symb. Comput. 43, 168 (2008). 13. S. A. Kochengin and V. I. Oliker, Numer. Math. 79, 553 (1998).
