## Abstract

A new optical component of transmissive faceted structure that allows the reshaping of collimated white light beam is depicted and realized. After being illuminated, each facet element in the structure can slightly tilt along its own axes to deflect the incident light. The calculation of the tilt angles is made by a numerical approach, which is the inverse solution of the analytical relationships between the two-dimensional tilt angles of the facet element and the local coordinates of the outgoing ray on the target plan. For a predefined illumination pattern on a fixed screen, the numerical computation of the 2D tilt angle matrix of the faceted structure is described. In the meantime, a Monte Carlo ray tracing program is created to calculate and verify the irradiance map, all compared with well-recognized commercial illumination software programs. Afterwards, the component is fabricated using an innovative additive technology, inspired by 3D printing and proposed by Luximprint.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The shaping of white light beam with high quality is important in non-imaging optics. Various kinds of optical elements have been proposed to reshape light with broad spectrum. Lightpipe is used for example in the micro display and is focused on the efficiency and uniformity of the light on the output plan [1–7]. In the areas of solar energy, non-imaging Fresnel lens have been designed as the concentrators [8–12]. Freeform optics for beam shaping involves a series of calculation and optimization methods [13–17]. But the fabrication of the surface profile of the free-form is still very complicated and expensive. Some other solutions are also possible in the field of faceted reflector. But the overall geometrical structure is often rotational symmetric, and arbitrary intensity profile has not been reported yet [18]. In computer graphics domain, Weyrich et al. create the idea of “microfacet height field” by series of algorithms to obtain a custom reflectance [19]. Independently, we see some faceted structures in micro dimensions called mirror cells arrays fabricated by two photon polymerizations in the recent literature [20,21].

We propose another original approach based on geometrical optics to allow white light beam shaping and to achieve arbitrary intensity profile on a fixed plan. The optical configuration contains a white light source, the transmissive faceted structure, and the detector. The faceted structure is composed of a matrix of square facets distributed in two directions. The geometrical dimensions of the facets are large with respect to the wavelength in this work.

The analytical mapping from the tilt angles of the facet element to the ray coordinates on the target plan is calculated. With the reverse solution of the analytical mapping, the 2D tilt angles matrix of the faceted structure can be computed numerically, when a desired illumination pattern is demanded. Then, using a house-made Monte Carlo raytracing program [22], the irradiance map can be calculated directly in Matlab (local illumination). Afterward, the calculated faceted model is transferred to the optical software Zemax OpticStudio 16 or LightTools 8.2 from Synopsis (global illumination), which confirm the light redistribution results on the detector. The local illumination result in Matlab and the global illumination result in illumination software programs are numerically compared with each other. The realization of the designed component with an innovative technology of 3D Printing proposed by Luximprint [23] is also demonstrated.

## 2. Optical configuration

The optical configuration is composed of a source of visible light (*S*), the transmissive faceted structure (*FS*) and the detector (*D*), given in Fig. 1. In our design, the incident beam from the source (*S*) is collimated and uniform. The transmissive faceted structure employs a matrix of *N _{Fx}* ×

*N*facets.

_{Fy}*FS*has a flat front surface and a faceted back one. The incident rays, which are normal to the front surface, are refracted by the output facets whose angles are modified along the axis

*x*and

*y*to get the desired illumination on

*D*.

Here, a complete coordinate system is defined: the global coordinates are (*x*, *y*, *z*), the local coordinates on the detector are (*X*, *Y*). *O* (0, 0, 0) is the center of the faceted structure. *D _{o}* is the center of the detector, with the global coordinates (0, 0,

*d*) and the local coordinates noted (0, 0).

_{FD}*d*is the distance between the rear base plan of

_{FD}*FS*and

*D*.

The vectors $\overrightarrow{s}$, ${\overrightarrow{n}}_{F}$ and ${\overrightarrow{n}}_{D}$ are normal to the plan of *S*, the back face of the facet element *F _{m,n}* and the plane

*D*. Here, the 2D local view of

*FS*in the

*x*-

*y*plan is presented. The 2D local view of

*D*opposite to ${\overrightarrow{n}}_{D}$ is also presented.

*F*denotes one facet element, (

_{m,n}*m*,

*n*) is the facet index:

*m*and

*n*is an integer in the range of [1,

*N*] and [1,

_{Fx}*N*]. The size of the facet element is

_{Fy}*p*and

_{x}*p*along the axis

_{y}*x*and axis

*y*. The size of

*D*is

*l*×

_{D}*l*mm

_{D}^{2}(or

*N*×

_{D}*N*pixels).

_{D}The ray mapping is important for a successful beam shaping. In Fig. 1, *F* (*m*, *n*, *u, v*) denotes the intersection point of the ray (*m*, *n*, *u*, *v*) on the back face of the facet *F _{m,n}*,

*D*(

*m*,

*n*,

*u, v*) denotes the intersection point of the ray (

*m*,

*n*,

*u*,

*v*) on the detector

*D*. (

*u*,

*v*) denotes the local index of the ray on the facet

*F*: the values of

_{m,n}*u*and

*v*are within [-1, 1]. Analytical equations are needed to solve the mathematical mapping between the ray intersection on the faceted structure $\overrightarrow{F}\left(m,n,u,v\right)$ in global coordinates (

*x*,

*y*,

*z*) and the ray intersection on the detector$\overrightarrow{D}\left(m,n,u,v\right)$in the local coordinates (

*X*,

*Y*). When

*F*tilts with (

_{m,n}*α*,

_{m,n}*β*) along the global axis

_{m,n}*x*and axis

*y*, a sketch to declare the position of

*F*(

*m*,

*n*,

*u, v*) is given in Fig. 2.

In Fig. 2, *F* (*m*, *n*, 0, 0) is the center of the facet *F _{m,n}*, whose position remains unchanged during the angular modification. In the

*y*-

*z*plan [Fig. 2(a)], the facet tilts with

*α*along the global axis

_{m,n}*x*(

*α*>0). In the

_{m,n}*x*-

*z*plan, [Fig. 2(b)], the facet tilts with

*β*along the global axis

_{m,n}*y*(

*β*>0). In the 3D space, when

_{m,n}*F*tilts with (

_{m,n}*α*) along the global axis

_{m,n}, β_{m,n}*x*and

*y*, the global coordinates $\overrightarrow{F}\left(m,n,u,v\right)$ are given by:

## 3. Analytical equations

The unit normal vector ${\overrightarrow{n}}_{F}$ of the facet element *F _{m,n}* is given as:

The unit vector of the refracted ray $\overrightarrow{t}$ is derived from the Snell’s law [16], and is given by:

*n*(resp.

_{1}*n*) denotes the refractive index of the faceted structure (resp. surrounding material). In our case, the surrounding material is air (

_{2}*n*= 1).

_{2}The unit vector of the parallel incident ray $\overrightarrow{s}$ is given by:

where*θ*(resp.

_{i}*θ*) is the incident (resp. refracted) angle.

_{t}The incident angle *θ _{i}* is calculated according to the dot product of ${\overrightarrow{n}}_{F}$ and $\overrightarrow{s}$, and is given by:

The refracted angle *θ _{t}* is calculated according to the Snell’s law in the trigonometric form and is given by:

*b*is given by:

To describe the trajectory of one refracted ray coming from the faceted structure to the detector, the mathematical line equation for one refracted ray in the 3D space is then given by Eq. (9):

In Eq. (9), *k* is a parameter for the line function,$\overrightarrow{F}\left(m,n,u,v\right)$ is already given by Eq. (1). The function of the plan of the detector is given by:

The intersection between the refracted ray and the plan of the detector is the ray incident point on the detector $\overrightarrow{D}\left(m,n,u,v\right)$. By solving the simultaneous set of Eq. (9) and Eq. (10), the value of the parameter *k* is calculated for $\overrightarrow{D}\left(m,n,u,v\right)$in global coordinates (*x*, *y*, *z*), as followed:

In Fig. 1, the local coordinates$\left(D\left(m,n,u,v\right)\cdot \overrightarrow{X},D\left(m,n,u,v\right)\cdot \overrightarrow{Y}\right)$ on the detector equals the global coordinates $\left(D\left(m,n,u,v\right)\cdot \overrightarrow{x},D\left(m,n,u,v\right)\cdot \overrightarrow{y}\right)$in the optical system. In conclusion, when the rear face of the facet element *F _{m,n}* tilts with (

*α*) along the axis

_{m,n}, β_{m,n}*x*and the axis

*y*, the local coordinates of the refracted ray (

*m*,

*n*,

*u, v*) on detector is given by:

## 4. Calculation of the local illumination result in Matlab

Section 3 gives the mathematical mapping (*f _{x}*,

*f*) between the 2D tilt angles (

_{y}*α*) and the local coordinates $\left(D\left(m,n,u,v\right)\cdot \overrightarrow{X},D\left(m,n,u,v\right)\cdot \overrightarrow{Y}\right)$ on the detector. With the reverse solution $\left({f}_{\alpha}^{-1},{f}_{\beta}^{-1}\right)$, the numerical computation of the local illumination result in Matlab is now explained. Here, we give a sketch to declare the general steps in the computation of the irradiance map in Matlab, and the parallel simulation in Zemax and LightTools.

_{m,n}, β_{m,n}In Fig. 3, the first step $\left({f}_{\alpha}^{-1},{f}_{\beta}^{-1}\right)$ is the numerical calculation of the tilt angles (*α _{m,n}, β_{m,n}*), depending on the desired illumination pattern, which is explained in section 4.1. The second step (

*f*,

_{x}*f*) is the numerical calculation of the local illumination result in Matlab, which is explained in section 4.2. In parallel, from the reverse solution $\left({f}_{\alpha}^{-1},{f}_{\beta}^{-1}\right)$, the calculated faceted model is also transferred into Zemax and LightTools. The non-sequential raytracing is progressed in the two-optical illumination software. So, the global illumination result in Zemax and LightTools is compared with the local illumination result in Matlab.

_{y}In Matlab, the ray is launched randomly, and is refracted only once on the surface. The local illumination only considers the ray path in the direct forward sequence without the real physical effects as reflection losses or diffraction [22]. In Zemax and LightTools, the non-sequential raytracing, considers multiple reflections, inter-reflection, retro-reflection, split, etc.

#### 4.1 Calculation of the tilt angles

After refraction from *FS*, the incident beam is divided into *N _{Fx}* by

*N*sub-beams. Due to the facet surface flatness considered here, each sub-beam (

_{Fy}*m*,

*n*) has one predefined direction. The desired illumination pattern is prescribed by the local coordinates of the center ray (

*m*,

*n,*0, 0) of each sub-beam. Here, the local coordinates are termed as $\left({X}_{D}\left(m,n\right),{Y}_{D}\left(m,n\right)\right)$. The 2D tilt angles (

*α*) can be calculated in Matlab with the numerical reverse solution $\left({f}_{\alpha}^{-1},{f}_{\beta}^{-1}\right)$:

_{m,n}, β_{m,n}The analytical approach in section 3 is developed to calculate the ray path after the inclined facets. Equation (12) calculates from the angles of the facets to the coordinates of the ray on the target plane, which are non-linear unordinary functions. So, inversing the equations (going from the desired ray coordinates on the target plane to retrieve the angles) is not possible with a simple analytical approach. Instead, we chose a numerical method to deal with the inverse problem in Eq. (13), which uses a specific function in Matlab called Vpasolve, inside the symbolic Toolbox.

#### 4.2 Monte Carlo raytracing program

From the first step of the numerical reverse solution $\left({f}_{\alpha}^{-1},{f}_{\beta}^{-1}\right)$, the tilt angle matrix is calculated. Afterwards, the second step is developed to calculate the local illumination result in Matlab.

A Monte-Carlo method is employed using the rand function inside Matlab. *N _{ray}* is the number of rays impinging on

*FS*with a random position (

*x*,

_{F}*y*), as follows:

_{F}*N*] is a 1 ×

_{ray}*N*vector which contain random values in the interval [0, 1]. In sequence, vectors

_{ray}*x*[1,

_{F}*N*] and

_{ray}*y*[1,

_{F}*N*] give the

_{ray}*x*and

*y*coordinates of

*N*rays on

_{ray}*FS*.

The difference between Eq. (14) and Eq. (1) is that, for Eq. (1), it calculates the ray coordinates on the structure, supposing the facet index (*m*, *n*), the local index (*u*, *v*) and the tilt angle (*α*, *β*) are known in advance. Relying on Eq. (1), the analytical relationship between the tilt angle (*α*, *β*) and the ray local coordinates on the detector is calculated and given in Eq. (12).

Equation (14) supposes the ray random positions on the faceted structure (*x _{F}*,

*y*) are known in advance. Afterwards, the facet index (

_{F}*m*,

*n*) and the local index (

*u*,

*v*) are derived with Eq. (15) and Eq. (16). Then, the ray local coordinates (

*X*,

_{D}*Y*) on the detector are calculated with the forward solution (

_{D}*f*,

_{x}*f*) in Eq. (12).

_{y}In Eq. (15), vectors *m*[1, *N _{ray}*] and

*n*[1,

*N*] denote the facet index, on which

_{ray}*N*rays incident, are given by:

_{ray}In Eq. (16), vectors *u*[1, *N _{ray}*] and

*v*[1,

*N*] calculate the ray local index on the facet element

_{ray}*F*, are given by:

_{m,n}Based on Eq. (15) and (16), the ray local coordinates on the detector are calculated with the forward solution in Eq. (12). To display the image, the local coordinates (*X _{D}*,

*Y*), whose value is within [-

_{D}*l*/2,

_{D}*l*/2] are rounded to the pixel index (

_{D}*X*,

_{p}*Y*), which are given by:

_{p}*X*,

_{p}*Y*) are integers and their values are within the interval between 1 and

_{p}*N*.

_{D}In Eq. (17), vectors *X _{p}*[1,

*N*] and

_{ray}*Y*[1,

_{p}*N*] record the pixel index of each ray on the detector. With the two vectors, the number of incident rays on every pixel of the detector is counted and recorded into a matrix with the dimension

_{ray}*N*×

_{D}*N*. In the incoherent mode, the optical power on pixel

_{D}*p*of the detector is determined by the number of rays incident on the pixel. The incoherent irradiance value on each pixel is the incident power over unit area of it. The matrix recording the number of rays on every pixel of the detector is converted to the matrix recording the incoherent irradiance value on every pixel, which is imaged to display the irradiance map in Matlab.

## 5. Design results

The local illumination result in Matlab and the global illumination result in Zemax and LightTools are presented in this section. The desired illumination pattern is a four-symmetric spaced squares (S4) on a plane, as shown in Fig. 4.

*FS* is made of 6 × 6 facets: the size of the facet element (*p _{x}* ×

*p*) is here 5 mm × 5 mm, to be large in comparison with the wavelength. The distance (

_{y}*d*) between the rear base plan of the faceted structure and the detector is 50 mm. The material (refractive index and transmission coefficient) of the faceted structure was defined according to the provided information from Luximprint, whose transparent 3D printing technology is adopted to realize the fabrication. The refractive index was used at the wavelength of 604 nm which corresponds to the primary wavelength of the source (Ocean Optics HL 2000) employed in our experiment. The raytracing configuration in Zemax and the 3D scheme of

_{FD}*FS*are shown in Fig. 5.

In the design (Table 1), the maximum and minimum feature height (*h _{max}*,

*h*) of the faceted structure is critical for the fabrication to follow. It needs to consider the

_{min}*z*coordinates of the four vertices (

*z*1,

*z*2,

*z*3,

*z*4) of each facet element

*F*. As the tilt angles of each facet is not the same, every facet element has different feature height. For the faceted structure matrix (

_{m,n}*N*×

_{Fx}*N*),

_{Fy}*h*is the maximum feature height of the

_{max}*N*×

_{Fx}*N*facets,

_{Fy}*h*is the minimum feature height of the

_{min}*N*×

_{Fx}*N*facets. In Fig. 6, the maximum feature height (

_{Fy}*h*) of the faceted structure is then 2.73 mm. According to Luximprint technology, the length of PMMA substrate is 2 mm.

_{max}Ocean Optics HL 2000 is applied as the source for the simulation: the spectrum is given in Fig. 7.

In Fig. 7, the peak intensity of Ocean Optics HL 2000 is at 604 nm. The spectrum has a strong output in the VIS-NIR band. With white light illumination, the 2D incoherent irradiance map is given in Fig. 8 with three different types of simulation.

The 2D irradiance matrix obtained with LightTools is transferred to Matlab. The 3D irradiance map and 1D irradiance diagram (*x* = 10mm, *y* = 10 mm) are plotted.

In Fig. 9, the 3D and 1D irradiance clearly reveal the high uniformity on the target plan. The high uniformity is obtained thanks to the collimated and uniform incident beam. The 2D angular modulations of the facets within such a small objective distance (50 mm) do not alter this property.

To evaluate the quality of the design, we calculate quality factors, which consist of the transmission (*T*), the efficiency (*η*), and the correlation coefficient (*C*). The transmission (*T* = *P _{D}* /

*P*) gives the optical power gained on the detector (

_{i}*P*) to the incident power on the faceted structure (

_{D}*P*). The efficiency factor (

_{i}*η*=

*P*/

_{T}*P*) measures the optical power in the target area (

_{D}*P*) to the total power on the detector (

_{T}*P*). It describes the ability of faceted structure to reshape and concentrate the light inside the targeting area. Correlation coefficient factor [24–27] is used here to assess the “difference” between the local illumination result in Matlab and the global illumination result in Zemax or LightTools. The function of the correlation coefficient (

_{D}*C*), uses four sub-factors:

_{MZ}or C_{ML}*I*(

_{M}*p*) is the incoherent irradiance per pixel from the local illumination result.

*I*(

_{Z}*p*) (resp.

*I*(

_{L}*p*)) is the incoherent irradiance per pixel from the global illumination result in Zemax or LightTools.

*I*(

_{M}*p*) and

*I*(

_{Z}*p*) (resp.

*I*(

_{L}*p*)) are normalized before calculation, so that the correlation coefficient (

*C*or

_{MZ}*C*) is in the interval [0, 1].

_{ML}In Table 2, the main reason of the difference between the local and the global illumination is that the real physical effects like interreflection, retroreflection, are not considered in the local illumination. The correlation coefficient between the two distinct types of illumination results is above 95%, and the efficiency factors are very similar in the two cases. Both the correlation coefficient and the efficiency factor confirm the quality of the numerical method, at the same time verify the very low influence of the facet edges and the shadowing effects on the distribution of the refracted beam on the detector. In particular, the low shadowing effect is maintained in our various design examples, which is not just limited to this case.

For example, for the illumination pattern of the four squares with a new arrangement in the blow Fig. 10 (a), the correlation coefficient between the local illumination result in Matlab and the global illumination result in LightTools is 96.62%; the efficiency is about 99.84% with our program, and about 99.92% in LightTools. For the illumination pattern of the rectangular hole in Fig. 10 (b), the correlation coefficient between the local illumination result in Matlab and the global illumination result in LightTools is 95.24%; the efficiency is about 99.99% with our program, and about 99.93% in LightTools.

In the above Fig. 10. (a), the size of the facet element (*p _{x}* ×

*p*) is 5 mm × 5 mm, the number of facets (

_{y}*N*×

_{Fx}*N*) is 8 × 8; on the detector, the size of each square (

_{Fy}*L*×

_{X}*L*) is 5 mm × 5 mm, the distance between each square is 40 mm. In Fig. 10. (b), the size of the facet element (

_{Y}*p*×

_{x}*p*) is 1 mm × 1 mm, the number of facets (

_{y}*N*×

_{Fx}*N*) is 20 × 20; on the detector, the size of the inner hole is 20 mm.

_{Fy}For different design examples, the self-shadowing and the edges of the facets may cause severe artefacts only when the angles of the facets are large compared to the working distance. In our case, whatever the desired image, the interval of the possible angles is controlled to the minimum in the design process: this means that the self-shadowing effect is restricted to be as small as possible. Moreover, for the presented design in Fig. 8, both the chosen interval of the angles which is [-15°, 15°] and the facet element dimensions are a compromise between the shadowing effects and the fabrication specifications coming from additive manufacturing Luximprint. For the illumination pattern in Fig. 10 (a), the interval of the angles is [-17°, 17°]; for the illumination pattern in Fig. 10 (b), the interval of the angles is [-5°, 5°].

## 6. Polychromatic effect

As our target is the beam shaping of white light, the change of the refractive index due to the large spectral bandwidth will affect the design result. The spectral sensitivity of the refracted ray is calculated here to evaluate this polychromatic impact. The spectrum of Ocean optics HL 2000 is [0.26, 1] µm. The discrete refractive index information provided by Luximprint is on a smaller interval. In Fig. 11, for each sampled wavelength, the vertical axis is the maximum of ray local coordinate deviation $\left(\Delta X,\Delta Y\right)$on detector due to the change of refractive index $\Delta {n}_{1}\left(\lambda \right)={n}_{1}\left(\lambda \right)-{n}_{1}\left(604\right)$, where the refractive index at 604 nm is employed as the basis.

From Fig. 11, within [0.26, 0.41] µm bandwidth, the maximum deviation is larger than one pixel (0.5mm/pixel). For Ocean Optics HL 2000 in Fig. 7, the intensity weight in [0.26, 0.41] µm is lower than 10 percent. So, the polychromatic effect within [0.26, 0.41] µm can be ignorable, which is also outside the visible region. For spectrum within [0.41, 0.87] µm, the ray coordinate deviation is lower than one pixel, which confirms the material dispersion has no impact on the design result.

Here, the chromatic dispersion effect is examined only for one fixed material, which is linked to the additive manufacturing technology. However, as we have stated previously, to limit the shadowing effect to a minimum level, the interval of the angles is controlled during the design process. At the same time, the prism dispersion effect is limited to a very low level. Chromatic effect will only be sensible if we have large interval angles. With the control of the interval of the angles, the possible material is not restricted to a specific one: it works for a variety of dispersive materials.

Except for the source polychromatic effect analysis, the impact of the source divergence and the variation of the working distance Δ*d _{FD}* on the quality of the illumination pattern are also evaluated. With the maintain of the quality decrease of the illumination pattern within 20 percent, the faceted structure can tolerant 2 degrees divergence of the incident beam, and 10 mm variation of the working distance. The quality decrease is measured by the four criteria: efficiency, uniformity, root means square error and correlation coefficient between the calculation result in Matlab and the simulation result in Zemax or LightTools.

## 7. Fabrication and experimental results

The fabrication of the transmissive faceted structure employed the transparent 3D printing technology from Luximprint. This optical 3D printing technology is inspired by the digital additive manufacturing process to fabricate typical array-typed component. Then, it is possible to use standard .step, .stl or .iges files to transfer the design information, which is under the requirements of the manufacturing process.

In the experimental setup, the employed source is Ocean Optics HL 2000 and is connected to an optical fiber P600-2-VIS-NIR, with a core diameter of 600 microns. Besides, the collimation of the output light is realized with a simple achromatic lens, the focal length *f* = 100 mm. In Fig. 12, the experiment setup in one dimension is given.

In Fig. 12, a square diaphragm manufactured by 3D Printing is put after the collimator to adjust the size of the beam. Benefits from the cooperation with Luximprint, a smooth surface profile fabrication technology was chosen. A tracing paper is put after the sample, which allows the image to be recorded by a digital camera. The main advantage of the tracing paper is the possibility to observe directly the irradiance map without parallax or distortion. The distance between the sample and the tracing paper *d _{FD}* is 50 mm. The used camera is a DSLR Canon Camera 1000D with a standard EF-S zoom objective 18-55 mm. The experiment results before and after the tracing paper are given in Fig. 13, and the efficiency

*η*in the target area is around 65%. Efficiency

*η*is calculated in Matlab platform, which uses a specific multicolor image processing.

## 8. Conclusion

This article demonstrates the design and the fabrication of a transparent faceted structure to do beam shaping of a white light source. The optical concept is based on a matrix of square facets inclined in two directions. Firstly, numerical solutions to calculate the tilt angles and the raytracing method to obtain the local illumination pattern on the target plan were explained in detail. Moreover, the global illumination results in Zemax OpticStudio 16 and LightTools 8.2 were compared with the local illumination result, which shows us the validity of our approach. The fabrication was accomplished with an innovative additive technology proposed by Luximprint with an optical transparent material: the dimensions of the facet element are 5 mm × 5 mm, and the component contains 6 × 6 facets. The experimental result of the irradiance map was recorded with a very simple optical setup and an estimation of the efficiency was also given using an image signal processing method to detect the target area. The optical concept of the faceted structure is proposed here in transmissive case. The reflective case is also conceivable with our method to take account different parameters like the off-axis angle of the source, the off-axis angle of the illumination screen, etc…

## Acknowledgments

The authors would like to acknowledge the financial support from University of Strasbourg for the fabrication of the element, the cooperation program with INSA Strasbourg, and the China Scholarship Council for the PhD Grant of Lihong Liu.

## References

**1. **N. Bokor and N. Davidson, “Anamorphic, adiabatic beam shaping of diffuse light using a tapered reflective tube,” Opt. Commun. **201**(4–6), 243–249 (2002). [CrossRef]

**2. **F. Fournier, W. J. Cassarly, and J. P. Rolland, “Method to improve spatial uniformity with lightpipes,” Opt. Lett. **33**(11), 1165–1167 (2008). [CrossRef] [PubMed]

**3. **J. F. Van Derlofske and T. A. Hough, “Analytical model of flux propagation in light-pipe systems,” Opt. Eng. **43**(7), 1503–1510 (2004). [CrossRef]

**4. **F. Fournier and J. Rolland, “Optimization of freeform lightpipes for light-emitting-diode projectors,” Appl. Opt. **47**(7), 957–966 (2008). [CrossRef] [PubMed]

**5. **I. Moreno, “Output irradiance of tapered lightpipes,” J. Opt. Soc. Am. A **27**(9), 1985–1993 (2010). [CrossRef] [PubMed]

**6. **C. M. Cheng and J. L. Chern, “Illuminance formation and color difference of mixed-color light emitting diodes in a rectangular light pipe: an analytical approach,” Appl. Opt. **47**(3), 431–441 (2008). [CrossRef] [PubMed]

**7. **C. H. Chen, C. C. Chen, and W. C. Liang, “Light pipe line beam shaper,” Opt. Rev. **14**(4), 231–235 (2007). [CrossRef]

**8. **J. J. O’Gallagher and R. Winston, “Nonimaging solar concentrator with near uniform irradiance for photovoltaic arrays,” Proc. SPIE **4446**, 52–56 (2002).

**9. **E. M. Kritchman, A. A. Friesem, and G. Yekutieli, “Efficient Fresnel lens for solar concentration,” Sol. Energy **22**(2), 119–123 (1979). [CrossRef]

**10. **R. Leutz, A. Suzuki, A. Akisawa, and T. Kashiwagi, “Design of a nonimaging Fresnel lens for solar concentrators,” Sol. Energy **65**(6), 379–387 (1999). [CrossRef]

**11. **W. T. Xie, Y. J. Dai, R. Z. Wang, and K. Sumathy, “Concentrated solar energy applications using Fresnel lenses: A review,” Renew. Sustain. Energy Rev. **15**(6), 2588–2606 (2011). [CrossRef]

**12. **K. Ryu, J. G. Rhee, K. M. Park, and J. Kim, “Concept and design of modular Fresnel lenses for concentration solar PV system,” Sol. Energy **80**(12), 1580–1587 (2006). [CrossRef]

**13. **J. C. Miñano, P. Benitez, and A. Santamaria, “Free-form optics for illumination,” Opt. Rev. **16**(2), 99–102 (2009). [CrossRef]

**14. **J. C. Miñano, P. Benítez, J. Blen, and A. Santamaría, “High-efficiency free-form condenser overcoming rotational symmetry limitations,” Opt. Express **16**(25), 20193–20205 (2008). [CrossRef] [PubMed]

**15. **D. Michaelis, P. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. **36**(6), 918–920 (2011). [CrossRef] [PubMed]

**16. **H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A **19**(3), 590–595 (2002). [CrossRef] [PubMed]

**17. **R. A. Hicks, “Controlling a ray bundle with a free-form reflector,” Opt. Lett. **33**(15), 1672–1674 (2008). [CrossRef] [PubMed]

**18. **S. R. David, C. T. Walker, and W. J. Cassarly, “Faceted reflector design for uniform illumination,” Proc. SPIE **3482**, 437–446 (1998). [CrossRef]

**19. **T. Weyrich, P. Peers, W. Matusik, and S. Rusinkiewicz, “Fabricating microgeometry for custom surface reflectance,” ACM Transactions on Graphics (Proc. SIGGRAPH) **28**, 321–326 (2009). [CrossRef]

**20. **D. Asoubar, C. Hellmann, H. Schweitzer, M. Kuhn, and F. Wyrowski, “Customized homogenization and shaping of LED light by micro cells arrays,” Proc. SPIE **9383**, 93831B (2015). [CrossRef]

**21. **E. Harnisch, M. Russew, J. Klein, N. König, H. Crailsheim, and R. Schmitt, “Optimization of hybrid polymer materials for 2PP and fabrication of individually designed hybrid microoptical elements thereof,” Opt. Mater. Express **5**(2), 456–461 (2015). [CrossRef]

**22. **G. Patow and X. Pueyo, “A survey of inverse surface design from light transport behavior specification,” Comput. Graph. Forum **24**(4), 773–789 (2005). [CrossRef]

**23. **Luximprint, Additive Fabrication Services for Functional and Decorative Optical Plastics, http://www.luximprint.com.

**24. **M. Flury, P. Gérard, Y. Takakura, P. Twardworski, and J. Fontaine, “Investigation of M2 factor influence for paraxial computer-generated hologram reconstruction using a statistical method,” Opt. Commun. **248**(4–6), 347–357 (2005). [CrossRef]

**25. **N. Bokor and Z. Papp, “Monte Carlo method in computer holography,” Opt. Eng. **36**(4), 1014–1020 (1997). [CrossRef]

**26. **N. Bokor and Z. Papp, “Computational method for testing computer-generated holograms,” Opt. Eng. **35**(10), 2810–2815 (1996). [CrossRef]

**27. **B. K. Jennison, J. P. Allebach, and D. W. Sweeney, “Iterative approaches to computer-generated holography,” Opt. Eng. **28**(6), 286629 (1989). [CrossRef]