## Abstract

Freeform optics generating specific irradiance distributions have been used in various applications for some time now. While most freeform optics design algorithms assume point sources or perfectly collimated light, the search for algorithms for non-idealized light sources with finite spatial as well as angular extent is still ongoing. In this work, such an approach is presented where the resulting irradiance distribution of a freeform optical surface is calculated as a superposition of pinhole images generated by points on the optical surface. To compute the required arrangement of the pinhole images for a prescribed irradiance pattern, the expectation maximization algorithm from statistics is applied. The result is then combined with a ray-targeting approach for finding the shape of the corresponding freeform optical surface. At its current state, the approach is applicable to Gaussian input irradiances, single-sided freeform optics and for the paraxial case. An example freeform optical surface for laser material processing is shown and discussed demonstrating the performance and the limitations of the presented approach.

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

## 1. Introduction

For some time now, freeform optics have been used in several application areas to generate specific prescribed irradiance distributions for a user-defined illumination [1]. By adapting a high number of degrees of freedom, the freeform optics’ surfaces are shaped to deflect the light of an incoming light source in such a way that a desired irradiance distribution forms on a given target surface. This has been exploited in many applications such as street lighting, car headlamps or laser materials processing [2–5].

There exists a number of different design algorithms that numerically derive the necessary shape of the optics’ surfaces, which mainly assume the incoming light to be either collimated or emitted from a point-like source. This is especially a valid approximation if the freeform optics’ size and its distance to the light source are large compared to the extent of the light source [6]. In a more scientific way, the quantity of étendue is introduced which measures the extent of the light source in phase space [7]. For a light bundle with constant luminance, the étendue is essentially defined as the integral of the product of the light’s extent and its propagation direction. For the light sources described above, the étendue becomes zero as either the extent is zero (point source) or the propagation direction is zero (collimated light). It is necessary to mention, that there are other possible representations of such idealized light sources with étendue being zero, which e.g. emit the rays in a fan.

The broad range of existing freeform optics design algorithms for these light sources includes the Supporting Quadric method [8,9], methods based on the solution of a Monge-Ampère type differential equation [10,11] or ray-targeting methods [12,13]. For the latter methods, a so-called ray-mapping is computed which specifies for each ray of the light source the target position it shall be deflected to so that in total the prescribed irradiance distribution is achieved. Then, the freeform surfaces are reconstructed from this information using e.g. an optimization algorithm. To obtain a smooth freeform optical surface, a so-called integrability condition for the resulting normal field must be met.

In the case of light sources that cannot be approximated with étendue-zero light, the design of freeform optics is much more challenging. This is due to the fact that for these light sources, a single point on the freeform surface is reached by several rays that each ideally requires a slightly different surface curvature which results in design conflicts and imperfect illumination [14]. Here, only a few successful design strategies have been described. In some cases, where the deviation from an étendue-zero light source is not too significant, the feedback algorithm [14,15] iteratively provides designs for a point source or collimated light, whose performance is then compared to the required irradiance distribution. Based on the results, the target irradiance distribution is adapted, and the design recomputed. Another class of algorithms uses direct optimization of the freeform optical surface in combination with an efficient approach of calculating the generated irradiance distribution [15–17]. Such approaches are for example based on the edge-ray principle [16] or on numerical integration of an irradiance integral [15,17]. These calculations are numerical challenging and often lead to irregular oscillating surfaces. In a recent publication [18], an approach for homogeneous illumination of polygonal targets is presented that first calculates the surface boundary and afterwards interpolates the surface interior to account for a homogeneous illumination. However, results are only presented for homogeneous light sources with well-defined edges and homogeneous target distributions. Finally, in the SMS3D method [19] several incoming wavefronts that might e.g. be emitted from the light source’s corner points are coupled to outgoing wavefronts from which the freeform surfaces can be reconstructed. In a special case of this approach, the wavefronts are coupled using a pinhole approximation [20]. In this case, a point on the freeform surface is approximated as a pinhole and the freeform optical surface is designed in such a way that the light exiting all the pinholes in a given direction equals the specified intensity value. In its current formulation, this method can only be applied for Lambertian light sources with well-defined edges. Additionally, the prescribed illuminance is described through an intensity (in W / sr) instead of an irradiance (in W / m^{2}) distribution.

Expectation maximization (EM) [21] is an iterative statistical approach for the estimation of model parameters in the existence of unknown latent variables. An exemplary and common application of this is estimating the parameters of Gaussian mixture models (GMM) [22]. Given a number of random samples drawn from a GMM, the application of EM yields the model parameters that most likely explain the sample distribution, i.e. a specification of its mean positions, the standard deviations or correlation matrices and the weights of the GMM component distributions. The unknown latent variables in this case are, for each given sample, the assignment of the Gaussian function that created it. In optics, EM has mostly been implemented for image reconstruction [23,24]. In a recent approach, EM has been applied to describe the four-dimensional phase space incident on a micro-structured freeform optics which enables the efficient sampling of rays for the subsequent ray-tracing [25].

In this work, EM is used for the design of freeform optics for extended light sources with a Gaussian radiant exitance distribution and a finite opening angle. It is used to compute a mapping that specifies the required target positions for a given set of rays emanating from the light source, corresponding to the mapping used in étendue-zero design algorithms. Similarly to [20], this approach is based on the pinhole approximation but the rays for the mapping are the rays from the center of the light source to the pinholes. In analogy to point source algorithms, the mapping is then used to reconstruct the freeform optical surface using ray-targeting.

Section 2 provides details of the simulation approach for which results are presented in Section 3. The algorithm uses many approximations that are discussed in Section 4. Section 5 concludes with a short summary of the work and an outlook.

## 2. Simulation method

The general idea exploited in this paper is the computation of the irradiance distribution generated by a freeform surface using the pinhole approximation (see Fig. 1). In this, the freeform surface is approximated by an adequate number of small pinholes (e. g. slightly larger than the number of degrees of freedom). Each pinhole generates an image of the light source whose position on a target plane is specified by the surface normal and by Snell’s law. The resulting irradiance distribution is then given by the superposition of all single source images.

To use this for the design of freeform optics, a two-step algorithm is suggested. At first, the required arrangement of the different source images for a prescribed irradiance pattern is computed using the expectation maximization algorithm. Afterwards, the result is used as input for a ray-targeting method using the source’s central rays.

#### 2.1 Step 1: Computation of the required image arrangement using EM

For the computation of the required image arrangement, the source is assumed to have a Gaussian radiant exitance pattern (as in Fig. 1 a)) and a constant angular emission pattern with a maximum emission angle. This means that the source irradiance follows a two-dimensional Gaussian distribution where the standard deviation ${\sigma _S}$ is related to the source’s radius ${r_L}$ (defined by that radial distance in which the relative irradiance – integrated over 2π – has dropped to $1/{e^2}$) through ${r_L} = 2 \cdot {\sigma _S}$. This choice is motivated by the application of laser beam shaping for high power materials processing. The assumption of constant angular emission on the other hand is simplified to some degree for depiction of the general idea.

The source is assumed further to be incoherent corresponding to high-power diode laser characteristics so that coherence and interference effects can be neglected.

Furthermore, a paraxial situation with small angles between the rays and the optical as well as the target surface is assumed. Then, the source images on the target plane can be supposed to be identical with equal weight and Gaussian with a constant standard deviation ${\sigma _T}$ which can be obtained from the input Gaussian’s standard deviation ${\sigma _S}$ and the geometrical setup ($n$: index of refraction, $\alpha :$ maximum source emission angle, ${r_{FF}}$: radius of the (round) freeform surface, ${d_{LT}}$: distance between freeform surface and target plane):

This relationship can be obtained from the fact that the beam parameter product stays constant upon refraction (see also [17]).

Under these assumptions, the expectation maximization algorithm is applied in the formulation of the solution of the (restricted) Gaussian mixture model problem (see e.g. [22]). Thus, each Gaussian with index *j* and mean ${\mu _j} = ({\mu _{j,x}},{\mu _{j,y}})$ is of the form

Here, ${\sigma _T}$ is a constant.

At first, *n* random 2D sample points are drawn with a probability distribution that is given through the (normalized) prescribed irradiance pattern. Thus, the density of samples corresponds to the brightness of the irradiance pattern and, in brighter areas, a larger number of samples are drawn than in darker areas. The expectation maximization now enables the fitting of *k* Gaussian distributions with a constant standard deviation ${\sigma _T}$ to the specified 2D point cloud. To this end, the EM algorithm utilizes a modification of the log-likelihood function that exhibits an identical lower bound behavior:

Herein, ${\gamma _{ij}}$ denotes the probability that the $i$-th sample is drawn from the $j$-th Gaussian function.

The starting point of the algorithm is then an initial distribution of the Gaussian functions. Here, a specific grid is chosen where the number of Gaussians equals the number of pinholes on the freeform surface. Furthermore, the grid of pinholes and the grid of initial Gaussians are associated in such a way, that there is a 1-to-1 relationship with the same neighboring characteristics (see Fig. 1 b)). This means that the initial mean positions of the Gaussians created by neighboring pinholes are neighbors themselves with the same connectivity.

The adaptation of this initial configuration is then performed in an iterative manner in two steps, according to the EM algorithm. In the first step of each iteration *m*, which is called the expectation step, the probability ${\gamma _{ij}}$ is estimated via

The second step, the maximization step, performs a maximum-likelihood fit of each Gaussian to the points associated with it. Using the representation from Eq. (2), the values for the mean positions in the next iteration are calculated analytically by maximizing $Q$

These two steps – associating the random points with one of the Gaussians and calculating the new mean positions – are performed iteratively until the achieved update in the maximum likelihood function lies below a certain threshold. The number of necessary iterations for a required threshold depends on the size of the pinhole images as the expectation maximization converges faster for larger standard deviations ${\sigma _T}$. For the examples in Section 3, the number of iterations ranges from ∼50 (for the examples in Fig. 3 and Fig. 6) to 252 (for the example in Fig. 5).

This approach yields the mean (i.e. center positions) of a number of 2D Gaussian functions whose superposition approximates the prescribed irradiance function.

It is worth to mention that the EM always yields a result but that only a convergence to a stationary value of the likelihood function in Eq. (3) is guaranteed [26]. This means that in our method, the GMM parameters may converge to a local maximum, the global maximum or saddle points. To partially remedy this local convergence behavior, it is possible to start with multiple initial parameter sets. In practical cases however, the solution of the EM step has been shown satisfactory without any adaptation of the initial guess from Fig. 1 b).

#### 2.2 Step 2: Ray-targeting for freeform optics reconstruction

The approach in the previous section yields an arrangement of a set of Gaussian distributions that approximate a prescribed irradiance pattern. From this, the required freeform optical surface can be reconstructed using ray-targeting schemes. As each of the pinholes is associated with one of the Gaussian functions, the target positions of the central rays (i.e. that ray incident on the pinhole coming from the center of the source) are directly given by the mean positions of the Gaussians and can be used in a conventional ray-targeting approach.

For this, the freeform optical surface is represented by a B-Spline surface whose degrees of freedom are given by the z-coordinates of the control points, where the z-axis is the main propagation direction. The positions of the pinholes are then fixed at equal distances in real space.

In the ray-targeting step, the least-squares sum for the intersection points of the central rays with the target ${t_j}$

is minimized using a standard minimization algorithm to fix the control points of the B-spline surface and reconstruct the corresponding freeform optical surface.At this point, it is important to highlight the difference to ray-targeting in a point-source algorithm. In the presented approach, the information of the extended source is included in the mapping itself as the finite extent and expected irradiance distribution of the source images is accounted for in the EM step. Thus, the mapping that is used for this ray-targeting step differs significantly from a point-source algorithm’s mapping for the same setup.

## 3. Simulation results

The presented algorithm is evaluated by designing a freeform optical surface for laser material processing that is intended to create an application-specific irradiance distribution as shown in Fig. 2(a) (The irradiance distributions shown throughout this paper are arbitrarily normalized in such a way that the source total power amounts to 1 W). The irradiance distribution has been calculated to realize an optimized temperature profile for the process of laser hardening [27,28]. For simulation purposes, the light source is assumed to be immersed within the optical material of the freeform lens which has an index of refraction of 1.5. For the incoming beam a Gaussian radiant exitance distribution with a radius of 4 mm and a finite divergence angle of 0.03333 rad are assumed which corresponds to a laser beam with 200 mm rad beam parameter product. The optical surface is round and has a radius of 10 mm from which follows a distance between light source and optical surface of 300 mm to accumulate all incoming radiation. Furthermore, the optical surface is designed to create the irradiance distribution on a target at a distance of 100 mm. From Eq. (1) follows that the standard deviation of the identical Gaussians on the target plane is 2.01 mm.

For this situation, the design algorithm of the previous section is applied using 2131 degrees of freedom for the freeform surface and the irradiance pattern is simulated assuming 2732 distinct pinholes. The number of pinholes is slightly larger than the number of degrees of freedom to obtain a well-defined problem. Figure 3 displays results of different steps within the design process. Figure 3(a) shows the distribution of the random points drawn with the irradiance distribution from Fig. 2(a) used as probability distribution. In Fig. 3(b), the resulting mean positions after the expectation maximization step are shown. Figure 3(c) additionally shows the irradiance distribution that is expected using this arrangement of the Gaussian images before the freeform optics design has been performed. Additionally, the image of a single Gaussian pinhole image is shown in the inset to visualize its extent. Figure 3(d) then presents the design of the resulting freeform optics which is very smooth. However, there are different zones or facets visible which originate in the systematic of lines and clusters in Fig. 3(b). These in turn arise as the irradiance distribution is a tiling of many source images, similar to results shown in [18]. Figure 3(e) displays the irradiance distribution that is obtained with this freeform optics as a superposition of the Gaussian images in the pinhole approximation after the ray-targeting has been performed. Herein, only the central rays’ position on the target plane and the Gaussian images as shown in the inset in Fig. 3(c) are used for the calculation. Finally, Fig. 3(f) presents an independent ray-tracing result of the freeform optics with a total power of 1 W.

The most obvious conclusion from these figures is that the resulting irradiances in Fig. 3(c), (e) and (f) deviate significantly from the prescribed irradiance pattern from Fig. 2(a). The edges are less sharp and the overall distribution is smoothed. These deviations are a natural consequence of the used extended source. As the image in Fig. 3(c) is composed of Gaussian distributions of significant extent (see inset), it is physically impossible to significantly decrease the difference to the prescribed irradiance in Fig. 2(a). However, the overall shape of the three boundary lines is strongly visible and especially the deviation between the resulting irradiance pattern of the EM step (Fig. 3(c)) and the ray-tracing result (Fig. 3(f)) is very low.

For comparison, a freeform optical surface has been designed using an algorithm for étendue-zero light sources [12] which assumes a point source with a corresponding emission pattern and opening angle. In a ray-tracing simulation, this freeform optical surface is irradiated with the real laser beam and the resulting irradiance distribution on the target plane is shown in Fig. 4(a). When compared to Fig. 3(f), it is obvious that the approach from Section 2 yields an irradiance pattern much closer to the prescribed one than using the étendue-zero approach. The edges as well as the peaks on the lower corners are much more prominent for the surface from the novel algorithm. This is quantified in Fig. 4(b) by comparing two cross-sections through the irradiance patterns along the x-axis at y = 0 mm. While the peak-to-valley ratio in the central area (marked by the range between the vertical dashed lines) is only 0.13 W/cm^{2} for the point source algorithm, the new algorithm creates an irradiance distribution with a peak-to-valley ratio of 0.39 W/cm^{2}. Furthermore, the edges at the lines’ ends are steeper for the presented algorithm (0.77 W/cm^{2} from x = -6 mm to x = -4 mm for the new algorithm vs. 0.41 W/cm^{2} from x = -6 mm to x = -4 mm for the point source algorithm). Both points indicate that with the new presented approach steeper edges can be obtained.

A special characteristic of the presented algorithms becomes visible when the radius for the situation in Fig. 2(b) is decreased while keeping the other parameters constant. In this case, the extent of the pinhole images is much smaller. Figure 5 shows the resulting irradiance pattern for a beam radius as small as 1 mm. The irradiance distribution from Fig. 2(a) is approximately obtained. A minor difference is that the edges are not perfectly straight (see discussion below for an interpretation). Thus, the algorithm also works for sources with small extent.

An important drawback of the presented approach becomes clear when a less paraxial situation is used for the formation of the prescribed irradiance distribution. To this end, the distance between the freeform lens and the target plane is reduced and the target size adapted accordingly to assure a comparable pinhole image size. Figure 6 shows ray-tracing results for freeform optical surfaces creating the same irradiance distribution in different distances. It becomes obvious that with decreasing distance the performance decreases and the irradiance distributions become blurred. This is a consequence of the non-paraxial nature of these geometric setups and the resulting distortion of the Gaussians on the target plane which are not incorporated in the design algorithm from Section 2.

## 4. Discussion

Although the previous section demonstrates a successful example of the applicability of the described approach, it also presents its limitations.

Most importantly, the expectation maximization step does not ensure that the calculated image positions lead to valid freeform optical surfaces. In other words, the integrability condition might be violated. The success of the design in Fig. 3 arises from the paraxial geometry, from the choice for the initial positions of the Gaussian distributions in the EM step as well as from the fact that the EM seems to keep neighboring characteristics for this application example constant. The last observation is not strictly guaranteed by the EM formulation used in this article. Therefore, the unfulfilled integrability condition might be an explanation for the difference between the distributions in Fig. 3(c) and (e) as well as the distortion of the lines in Fig. 5. Thus, for more complex application cases, a more rigorous investigation of the integrability condition by adapting the EM step will be necessary.

A restriction to paraxial cases arises as the pinhole images are assumed to be identical. This is a strong restriction as it also includes the assumption that the angular emission pattern of the light source is homogeneous. Especially in non-paraxial situations, the source images on the target plane will vary with position which depends on the one hand on the position of the pinhole on the freeform optical surface and on the other hand on the induced deflection. In the presented examples from Fig. 6, this is more restricting than the possibly violated integrability condition.

Additionally, the EM converges to a local minimum, which only yields an approximately satisfying approximation of the prescribed irradiance distribution by Gaussian functions. There seem to be still deviations that can be traced back to the local minimization that is obtained with EM and possibly also to the finite number of Gaussians. i.e. pinholes. One consequence is for example that there are minor irregularities on the freeform optical surfaces itself which become visible when investigating the surface curvature in a CAD program. Running the EM with different initial positions of the Gaussians substantiates that. However, for now, the correspondence between prescribed and resulting irradiance distribution is satisfying enough. This especially holds as the major difference between Fig. 3 and Fig. 2(a) arises due to the nature of the extended source and is independent of the EM result.

Furthermore, in the examples above, only a single freeform surface is designed. Studying two-sided freeform lenses with a significant influence of the entrance surface requires a thorough investigation of the image propagation.

Finally, the presented approach is only validated for Gaussian input distributions. This is especially convenient as the second EM-step yields an analytic expression for the new mean positions of the Gaussian distribution. However, by using other mixture model techniques from statistics, it seems possible to formulate adequate expressions for a broader range of input distributions e.g. by using super-Gaussian distributions or a description of the light source as a superposition of Gaussian functions.

## 5. Summary

A promising idea for the design of freeform optics for extended sources has been presented that combines the expectation maximization algorithm from statistics in the formulation to solve a Gaussian mixture model with a ray-targeting approach. The performance of the algorithm has been shown using an irradiance distribution for laser material processing. Here, sharper features are obtained than using a point source algorithm for the studied geometry with an increase in peak-to-valley values up to a factor 3. However, due to the relevant assumptions made, there are strong limitations to the given approach. These result especially in non-paraxial cases in irradiance distributions that differ significantly from the prescribed ones. These deviations probably are due to the violation of the integrability condition and to the fact that for non-paraxial cases the Gaussians cannot be assumed to be identical. The adaptation of the algorithm to non-paraxial cases will be the content of future work.

## Funding

Deutsche Forschungsgemeinschaft (EXC-2023 Internet of Production – 390621612).

## Acknowledgments

This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2023 Internet of Production – 390621612.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **R. Wu, Z. Feng, Z. Zheng, R. Liang, P. Benítez, J. C. Miñano, and F. Duerr, “Design of freeform illumination optics,” Laser Photonics Rev. **12**(7), 1700310 (2018). [CrossRef]

**2. **Y. C. Lo, K. T. Huang, X. H. Lee, and C. C. Sun, “Optical design of a Butterfly lens for a street light based on a double-cluster LED,” Microelectron. Reliab. **52**(5), 889–893 (2012). [CrossRef]

**3. **A. Timinger and H. Ries, “Street-lighting with LEDs,” Proc. SPIE **7103**, 71030H (2008). [CrossRef]

**4. **M. Berens, A. Bruneton, A. Bäuerle, M. Traub, R. Wester, J. Stollenwerk, and P. Loosen, “Multiple intensity distributions from a single optical element,” Proc. SPIE **8834**, 88340M (2013). [CrossRef]

**5. **F. Klocke, M. Schulz, and S. Gräfe, “Optimization of the laser hardening process by adapting the intensity distribution to generate a top-hat temperature distribution using freeform optics,” Coatings **7**(6), 77 (2017). [CrossRef]

**6. **I. Moreno and C. C. Sun, “LED array: where does far-field begin?” Proc. SPIE **7058**, 70580R (2008). [CrossRef]

**7. **J. Chaves, * Introduction to nonimaging optics*, 1

^{st}ed. (CRC, 2008).

**8. **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]

**9. **V. Oliker, “Controlling light with freeform multifocal lens designed with supporting quadric method (SQM),” Opt. Express **25**(4), A58–A72 (2017). [CrossRef]

**10. **Z. Feng, B. D. Froese, R. Liang, D. Cheng, and Y. Wang, “Simplified freeform optics design for complicated laser beam shaping,” Appl. Opt. **56**(33), 9308–9314 (2017). [CrossRef]

**11. **J. H. M. ten Thije Boonkkamp, C. R. Prins, W. L. IJzerman, and T. W. Tukker, “The Monge-Ampere equation for freeform optics,”. In Freeform Optics (2015), pp. FTh3B-4.

**12. **A. Bäuerle, A. Bruneton, R. Wester, J. Stollenwerk, and P. Loosen, “Algorithm for irradiance tailoring using multiple freeform optical surfaces,” Opt. Express **20**(13), 14477–14485 (2012). [CrossRef]

**13. **Z. Liu, P. Liu, and F. Yu, “Parametric optimization method for the design of high-efficiency freeform illumination system with a LED source,” Chin. Opt. Lett. **10**(11), 112201 (2012). [CrossRef]

**14. **Y. Luo, Z. Feng, Y. Han, and H. Li, “Design of compact and smooth freeform optical system with uniform illuminance for LED source,” Opt. Express **18**(9), 9055–9063 (2010). [CrossRef]

**15. **R. Wester, G. Müller, A. Völl, M. Berens, J. Stollenwerk, and P. Loosen, “Designing optical freeform surfaces for extended sources,” Opt. Express **22**(S2), A552–A560 (2014). [CrossRef]

**16. **A. Hirst, J. Muschaweck, and P. Benítez, “Fast, deterministic computation of irradiance values using a single extended source in 3D,” Opt. Express **26**(14), A651–A656 (2018). [CrossRef]

**17. **A. Völl, R. Wester, M. Berens, P. Buske, J. Stollenwerk, and P. Loosen, “Accounting for laser beam characteristics in the design of freeform optics for laser material processing,” Adv. Opt. Technol. **8**(3-4), 279–287 (2019). [CrossRef]

**18. **D. A. Birch and M. Brand, “Design of freeforms to uniformly illuminate polygonal targets from extended sources via edge ray mapping,” Appl. Opt. **59**(22), 6490 (2020). [CrossRef]

**19. **P. Gimenez-Benitez, J. C. Miñano, J. Blen, R. M. Arroyo, J. Chaves, O. Dross, M. Hernandez, and W. Falicoff, “Simultaneous multiple surface optical design method in three dimensions,” Opt. Eng. **43**(7), 1489–1503 (2004). [CrossRef]

**20. **S. Sorgato, J. Chaves, H. Thienpont, and F. Duerr, “Design of illumination optics with extended sources based on wavefront tailoring,” Optica **6**(8), 966–971 (2019). [CrossRef]

**21. **G. J. McLachlan and T. Krishnan, * The EM algorithm and extensions*, 2

^{nd}ed, (John Wiley & Sons, 2007).

**22. **M. R. Gupta and Y. Chen, * Theory and use of the EM algorithm*, (Now Publishers Inc, 2011), Chap. 3.

**23. **K. T. Lay and A. K. Katsaggelos, “Image identification and restoration based on the expectation-maximization algorithm,” Opt. Eng. **29**(5), 436–446 (1990). [CrossRef]

**24. **N. Cao, A. Nehorai, and M. Jacob, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express **15**(21), 13695–13708 (2007). [CrossRef]

**25. **M. Berens, A. Völl, R. Wester, J. Stollenwerk, and P. Loosen, “Density estimation in optical phase space for optimizing micro-optical elements on freeform surfaces,” Proc. SPIE **10694**, 106940N (2018). [CrossRef]

**26. **C. F. Jeff Wu, “On the Convergence Properties of the EM Algorithm,” Ann. Statist. **11**(1), 95–103 (1983). [CrossRef]

**27. **A. Völl, J. Stollenwerk, and P. Loosen, “Computing specific intensity distributions for laser material processing by solving an inverse heat conduction problem,” Proc. SPIE **9741**, 974105 (2016). [CrossRef]

**28. **D. Burger, “Beitrag zur Optimierung des Laserhärtens,” Ph.D. dissertation (University Stuttgart, 1988).