## Abstract

Inverse lithography technology (ILT) treats photomask design for microlithography as an inverse mathematical problem. We show how the inverse lithography problem can be addressed as an obstacle reconstruction problem or an extended nonlinear image restoration problem, and then solved by a level set time-dependent model with finite difference schemes. We present explicit detailed formulation of the problem together with the first-order temporal and second-order spatial accurate discretization scheme. Experimental results show the superiority of the proposed level set-based ILT over the mainstream gradient methods.

©2009 Optical Society of America

## 1. Introduction

#### 1.1. Inverse lithography

Optical microlithography, a critical step in the integrated circuit (IC) fabrication process, is increasingly challenging with a continuous shrinkage of the circuit minimum feature size. The bandlimited imaging system acts as a lowpass filter and introduces significant distortions caused by diffraction. Several resolution enhancement techniques (RETs) [1–3] have been developed to date; taking advantage of the increasing computational power and optimization algorithms, inverse lithography has shown great promise to meet the challenges in future technology nodes.

Inverse lithography synthesizes an “optimal” mask given a desired image, and is unconstrained by the topology of the original physical design [4, 5]. It has been pointed out that while the specifications of inverse lithography make engineering sense, the current performance still falls short, making this a “virtual virtuality” as it is not yet applicable to real-world manufacturing [6, 7].

Nevertheless, various mathematical techniques have been put forward to tackle inverse lithography. One of the earliest work was by Liu and Zakhor, who applied the branch and bound algorithm in [8] and what they called the “bacteria” algorithm in [9] for an optimal binary and phase shifting mask (PSM) design. These methods, however, are very time-consuming, which limits their applicability. In the mid-1990s, iterative methods were applied to generate binary masks [10] and Pati and Kailath used Gercberg-Saxton algorithms to generate PSMs [11]. In recent years, Granik classified methods for solving inverse mask problems on linear, quadratic and non-linear formulations [12], based on which a model-based technique for SRAF insertion is developed [13]; Poonawala and Milanfar designed a model-based optical proximity correction (OPC) system and developed a faster optimization algorithm using steepest-descent [14, 15]. Similarly, in [5, 16–18], iterative methods for PSMs are also presented, treating the problem as an inverse problem and solving it by an optimization process. This has also been expanded to include pixelated mask technology [19], double exposure lithography [20], manufacturability enhancement [21] and designs specifically taking into account robustness to variations, such as in focal length [22]. Meanwhile, level set approach has also been actively explored as a feasible alternative [4, 23]. However, the exact formulations and technical discussions are not provided in sufficient details in the open literature, making it hard to assess the merits and drawbacks of this approach. In this paper, we provide a complete level set formulation of the problem, together with a computationally efficient numerical solution to the proposed level set-based formulation.

#### 1.2. Level set method

Level set as a mathematical technique is pioneered by Osher [24] and Sethian. Optical lithography is in fact one of the earlier applications considered by the latter [25–28]. However, these papers aimed to compute a lithographic development profile for a given evolving front. In other words, this is “forward lithography” in contrast to the inverse problem we treat in this paper. Meanwhile, level set methods have been investigated for various image processing problems [29], including inverse imaging. For example, level set has been used to tackle inverse problems involving obstacles [30, 31], and in nonlinear deblurring and noise removal where a time-dependent model is developed, using the total variation (TV) norm as the regularization functional [32].

We address the inverse lithography problem as an extension of level set motion for nonlinear image restoration, or alternatively, as an obstacle reconstruction problem. This is different from conventional deblurring methods in image processing, as we have to take into account the resist effects which impose a nonlinear threshold-like effect on the circuit pattern. We also present the formulation with sufficient details so others can follow, analyze, and improve, with the goal of fostering more open discussion on ways to apply level set to optical lithography. Furthermore, we develop a novel approach to solve the level set formulation with high fidelity and robustness, by incorporating the time-dependent model with available difference schemes.

## 2. The Constrained level set time-dependent model formulation

#### 2.1. Mathematical model

For simplicity in both mathematics and implementation, we concentrate on coherent sources and binary masks, but our methodology can be readily applied to partially coherent illumination and PSMs. The projection lithography imaging process can be described as a forward model

where the boldface **x** denotes spatial coordinates (*x,y*), and *𝓣* {·} maps the input intensity function *U*(**x**) to the output intensity function *I*(**x**). Due to the lowpass nature of the optical lithography imaging system, *I*(**x**) is typically a blurred version of *U*(**x**). Suppose we denote the desired circuit pattern as *I*
_{0}(**x**). The objective of inverse lithography is to find a predistorted input intensity function *U*̂(**x**) which minimizes its distance with the desired output, i.e.,

in which *d*(·, ·) is an appropriately defined distance metric, such as the *ℓ*
_{2} norm.

The lithography process *𝓣*{·} in Eq. (1) can be divided into two functional blocks, namely the projection optics effects (aerial image formation) and resist effects. If the optical system is coherent, the aerial image intensity distribution *I*
_{aerial}(**x**) is related to the input intensity function *U*(**x**) and the point spread function (PSF) *H*(x) by

where * denotes convolution, and *H*(**x**) is taken as the inverse Fourier transform of a disc function [33]. The resist effects can be approximated using a logarithmic sigmoid function [15]

with *a* being the steepness of the sigmoid and *t _{r}* being the threshold. Putting together, we can write the image formation equation for a coherent imaging system as

In what follows, we will drop the arguments **x** whenever there is no ambiguity.

#### 2.2. Optimization framework

From Eqs. (3) and (5), we can see that inverse lithography involves deconvolution and is generally ill-posed. One possibility is to consider a variational formulation of the model that regularizes the problem by incorporating the formulation with a functional *𝓡*. This functional takes as its argument an image *U* and returns a scalar that quantifies the quality of the pattern. Smaller values of *𝓡* correspond to more desirable images. *𝓡*
_{1}, *𝓡*
_{2}, and *𝓡*
_{3} below are common examples of the regularization (note that ∇ denotes the gradient, Δ denotes the Laplacian [32], and Ω is an open and bounded domain of R^{2} which contains *U*):

*𝓡*
_{3} above is called the total variation (TV) regularization. Quadratic penalty terms and complexity penalty terms can also be used [15]. Since TV norm does not penalize discontinuities in *U*, thus allowing the recovery of its edges, we opt to use it in this work. The inverse lithography

problem can be written as

$$\mathrm{subject}\phantom{\rule{.2em}{0ex}}\mathrm{to}\phantom{\rule{1em}{0ex}}{\int}_{\Omega}{\left(\mathrm{sig}\left({\mid H*U\mid}^{2}\right)-{I}_{0}\right)}^{2}\text{d}\mathbf{x}=\epsilon ,$$

where *ε* is the parameter controlling the approximation of the aerial image and the hard thresholding to the desired circuit pattern, with a smaller *ε* indicating better approximation. The Lagrangian of the above expression is

$\lambda {\int}_{\Omega}\mid \nabla U\mid \text{d}\mathbf{x}+\frac{1}{2}({\int}_{\Omega}{(\mathrm{sig}\left({\mid H*U\mid}^{2}\right)-{I}_{0})}^{2}\text{d}\mathbf{x}-\mid \Omega \mid \epsilon ),$,

where |Ω| is the area of the bounded domain, or number of pixels in *U*, and its Euler-Lagrange equations, with homogeneous Neumann boundary conditions for *U*, are

where *α*(**x**) is defined as

$$=\alpha \{H*\left[({I}_{0}-\mathrm{sig}(H*U))\odot \mathrm{sig}(H*U)\odot \left(1-(H*U)\right)\odot (H*U)\right]\},$$

and ⊙ denotes element-by-element multiplication. The derivation of *α*(**x**) in Eq. (12) requires a few steps and can be found in [5] and [15]. If we assume *λ* to be known, the equivalent unconstrained problem of Eq. (9) can be written as

and its Euler-Lagrange equation is the same as Eq. (10).

#### 2.3. Time-dependent scheme

A stable explicit time-dependent scheme can be applied to solve Eq. (10), where the choice of parameters is almost user-independent and the programming is comparatively simple, yielding

with *t* denoting time and *α*(**x**) also updated to a time-dependent function. The solution procedure is a parabolic equation with time as an evolution parameter. Equation (14) moves each level curve of *U* normal to itself, with velocity equal to the curvature of the level surface divided by the magnitude of the gradient of *U*. The constraints are included in *α*(*x, t*) to prevent distortion and to obtain a nontrivial steady state.

A disadvantage of this approach, however, is that it converges slowly to steady state, since the parabolic term is close to singular for small gradients. To cope with this, we use an improved time-dependent model which accelerates the movement of the level curves of *U* and regularizes the parabolic term in a nonlinear way, by multiplying the right hand side of Eq. (14) with the magnitude of the gradient [24, 32]. Thus, we have instead

This can be viewed as a convection-diffusion equation, with morphological convection and anisotropic diffusion [32]. It has a clearer physical meaning than Eq. (14) because it moves each level curve of *U* normal to itself with velocity equal to -*α*(*x, t*), combined with a curvaturedriven flow characterized by velocity $\lambda \nabla \xb7\left(\frac{\nabla U}{\mid \nabla U\mid}\right)$. Preconditioning Eq. (15) for speed improvement is also possible [29, 32].

#### 2.4. Obstacle reconstruction

We can also treat the inverse lithography problem as an obstacle reconstruction problem [30, 31], which is an inverse problem involving obstacles where the desired unknown is a region consisting of several subregions. We give *U* a level set description by introducing an unknown function *ϕ*(**x**), which is related to *U* by defining

where *U*
_{int}=1 and *U*
_{ext}=0 if we are dealing with binary masks. Now we can define the inverse lithography problem as finding *ϕ*(x) such that *𝓣*(*U*)≈*I*
_{0}. Solving this with a least squares fit to the approximation is equivalent to seeking the minimizer of

The boundary of the subregions in *U* where *U*=*U*
_{int} is governed by the zero level set of *ϕ*, namely, *ϕ*(**x**)=0. By taking the variation of the equation *ϕ*(**x**)=0, we find

where *δϕ* and *δ*
**x** represent the small variation of *ϕ* and **x**, respectively. If *ϕ* depends on both **x** and time *t*, we associate the evolution of the subregions in *U* with *ϕ*(**x**, *t*). The minimum requirement for *δϕ*(**x**, *t*) is that *F*(*U*) be a decreasing function of time, i.e., the resulting variation of *F*(*U*) or *δF*(*U*) is negative. Assuming that each point moves perpendicular to the surface with normal velocity *α*(**x**, *t*), *δ*
**x** satisfies

According to [30], defining *α*(*x, t*) as

where *J*(*U*) is the Jacobian of *𝓣*(*U*) at *U*, will produce a small variation *δϕ* that reduces *F*(*U*), which leads to

Note that *α*(**x**, *t*) is exactly the same as Eq. (15). The above equation is different from it only without the convection part, and Eq. (15) is easily derived from Eq. (20) by adding a normal speed which is λ times the curvature of the level surface.

As such, Eq. (20) is more general in the sense that adding regularization terms other than TV norm will result in different normal speeds. Meanwhile, adding TV norm tends to reduce noise in the image by minimizing the length of the contours. Moreover, reducing the TV norm favors cutting the corners in *U*, which is undesirable in photomasks [34]. Therefore, we will use Eq. (20) to solve the inverse lithography problem in this paper. Besides, when enforcing mask manufacturability rules, the inverse lithography solution becomes a constrained optimization if we regularize the inverse problem, and using the above equation will be more convenient in terms of describing the geometrical constraints, for example, regularization terms including, but not limited to, mean curvature regularization (TV norm) and average mean curvature affect the geometry of the contour, the former minimizes the length of the contours while the latter minimizes the length of the boundaries keeping the total area of the enclosed regions fixed, giving preference to smoother contours and contours enclosing larger regions.

We also remark that the complexity of computing *α*(**x**, *t*) is governed by the convolution operation. Therefore, the algorithmic complexity is *𝒪*(*N*
^{2} log*N*), assuming an *N*×*N* image.

#### 2.5. Numerical schemes

Equation (20) is a partial differential equation (PDE). Once *ϕ* and *α* are defined at every grid point on the Cartesian grid, the first-order derivatives in space and time of Eq. (20) can be approximated using finite difference techniques. An appropriate difference scheme should be designed for our level set approach. It is natural to use a higher order accurate difference scheme, for example, the 2–4 implicit scheme in [35], if we want to achieve a better input mask. However, more accurate difference scheme requires more computation, and therefore a tradeoff should be made between accuracy and computation time.

It is possible to evaluate the spatial derivatives of *ϕ* using a first-order accurate forward difference and a backward difference, or a second-order accurate central difference, together with upwind differencing [36] which chooses an approximation to the spatial derivatives based on the sign of *α*. However, we can improve the first-order accurate upwind scheme significantly by a more accurate approximation of the spatial derivative of *ϕ*. Harten et al. develop the essentially nonoscillatory (ENO) polynomial interpolation method that uses the smoothest possible polynomial interpolation to find *ϕ*, and then differentiates it to get the spatial derivatives [37]. The implementation is further refined in [38, 39]. Frequently used ENO schemes include first-, second-, and third-order accurate essentially nonoscillatory (ENO1, ENO2, ENO3) methods, as well as fifth-order accurate weighted ENO (WENO) method that takes a convex combination of ENO approximations. Interested readers may refer to [24, 29, 32, 36] for more details.

For the temporal discretization in Eq. (20), a rather simple first-order accurate method is the forward Euler method. However, to achieve a higher-order temporal discretization, we can use the total variation diminishing (TVD) Runge-Kutta (RK) method [38], where the basic first-order accurate TVD RK scheme reduces to the forward Euler method. In this paper, we aim to find a consistent finite difference approximation to Eq. (20) making a tradeoff between accuracy and computation time by comparing experimental results using combinations of different temporal and spatial difference schemes.

## 3. Experimental Results

We apply the inverse lithography technique outlined above in designing various circuit patterns, and compare the results that use different temporal and spatial accuracies. We use the same imaging system parameters: *λ*=193nm, *NA*=0.85, resolution=10nm/pixel, thresholding *t _{r}*=0.3, and therefore the same PSF

*H*(

**x**) for various experiments. The optimization approach stops after 50 iterations.

Consider Fig. 1. Each row has three sub-images, where the first column is the input mask pattern *U*(**x**), the second column is the corresponding aerial image *I*
_{aerial}(**x**), and the third column is the binary output circuit pattern *I*(**x**). The various rows denote different input patterns. In (a), we use the target circuit pattern as input, while in (b)–(e), we design the input mask using the target circuit pattern as input to our level set algorithms. All of them use a first-order temporal accuracy, but the difference is that (b) uses an ENO1 spatial accuracy, (c) uses ENO2, (d) uses ENO3 and (e) uses WENO. All images are of size 170×140 pixels, while *H*(**x**) is a jinc function of size 51×51.

As we would expect, using the desired circuit pattern as the input to the lithographic system produces unacceptable circuit patterns due to the blurring in the imaging system. Using the mask computed by our level set algorithm produces a circuit pattern closer to the intended one, as shown in Fig. 1(b)–1(e). However, the pattern in (b) is still unacceptable because of the merging patterns in the lower right circled area. Improving the spatial accuracy to ENO2, ENO3 and WENO, as are the cases in Fig. 1(c)–1(e), produce better masks and the resultant output patterns are closer to the intended circuit. We can also examine the pattern error numerically, by counting how many pixels are different between the intended circuit pattern and the one obtained after the simulated lithographic system. This is given in Table 1. Computation time in Fig. 1(c) which uses spatial accuracy ENO2 is set as 1.00, against which other computation times are normalized. It can be seen that increasing spatial accuracy from ENO1 to ENO2 brings a significant decrease in pattern error, which agrees with the visual analysis earlier, while further spatial accuracy increase from ENO2 to ENO3 and WENO does not lead to a considerable decrease in pattern error at the cost of increased computation time. Figure 1 and Table 1 is typical of such simulation experiments with different input pattern masks producing similar visual observation, pattern error and computation time analysis, suggesting that better input mask is designed using the proposed level set approach and ENO2 spatial differential scheme as demonstrated in the numerical experiments.

Another set of simulation is given in Fig. 2, this time with a regular array of rectangular shapes representing the contact layer with size 140×140. Illustration of sub-images in various rows and columns is the same as that in Fig. 1. In Fig. 2(b)–2(d) where the target circuit pattern is used as the input to the proposed level set approach, all using ENO2 spatial accuracy, but different temporal schemes, i.e., first-order temporal accuracy using forward Euler methods in (b), second- and third-order temporal accuracy using TVD RK methods in (c) and (d). The pattern error of Fig. 2 is presented in Table 2. In Table 2, computation times in Fig. 2(c) and 2(d) are normalized against that in Fig. 2(b). From Fig. 2 and Table 2, the masks computed using the level set approach improve the closeness of the circuit pattern to the intended one. We can see that more computation time comes along with higher order temporal discretization, but without much improvement in pattern error. These observations in Fig. 2 are applicable to simulation experiments with different input pattern masks which suggest first-order temporal accuracy is the right choice in the proposed level set approach. Combining the first-order temporal accuracy with the aforementioned ENO2 spatial accuracy, we have decided the difference schemes of the proposed level set approach.

Gradient methods play an important role in solving the inverse lithography problem in Eq. (2) [5, 15], which can be either steepest descent or conjugate gradient. In Fig. 3, the restored input pattern masks using conjugate gradient method [40] and level set methods are presented. Table 3 presents the pattern error and computation time of Fig. 3, where the computation time in Fig. 3(b) is normalized against that in Fig. 3(c) while the computation time in Fig. 3(d) is normalized against Fig. 1(c). Comparing the two sets of restoration performances in Fig. 3(d) and Fig. 3(b) using gradient method with Fig. 1(c) and Fig. 3(c) using the proposed level set approach, together with normalized computation time and pattern error in Table 3, it is concluded that the proposed level set methods achieve a considerably better performance in pattern error than gradient methods, with almost equal computation time. It is also noted the restored input pattern masks tend to eliminate small unwanted block objects in the input pattern masks generated by gradient methods which is of great importance in mask manufacturing.

It should be noted that upon mask manufacturing rules, often geometric constraints are enforced. With its ability to handle topological merging and breaking when the topology is unknown, the level set approach to inverse lithography is in a much better position to describe the geometric constraints than conventional inverse lithography. Therefore, the potential of level set approach to inverse lithography upon mask mask rule checking (MRC) and user-defined requirement is substantial.

## 4. Conclusion

We have formulated a level set approach for photomask design in optical microlithography as an obstacle reconstruction problem or an extended nonlinear image restoration problem. Detailed level set formulations for inverse lithography have been explicitly shown which are not found elsewhere in the literature. The design of a predistorted mask to counteract the optical distortion is considered as an inverse mathematical problem and falls into the level set optimization framework by describing the inverse problem as a constrained time-dependent model PDE which is solved by finite difference schemes. Experiments of solving the PDE with different orders of temporal and spatial accuracy are conducted, suggesting first-order temporal accuracy and second-order essentially nonoscillatory spatial accuracy are appropriate difference schemes in the proposed level set method. Comparison with existing gradient methods shows the superiority of the proposed level set approach in pattern error improvement and elimination of small unwanted block objects in mask manufacturing.

## Acknowledgment

This work was supported in part by the Research Grants Council of the Hong Kong Special Administrative Region, China under Projects HKU 7139/06E, 7174/07E and 7134/08E.

## References and links

**1. **A. K.-K. Wong, *Resolution Enhancement Techniques in Optical Lithography* (SPIE Press, Bellingham, WA, 2001).

**2. **F. Schellenberg, “Resolution enhancement technology: the past, the present, and extensions for the future,” Proc. SPIE **5377**, 1–20 (2004).
[CrossRef]

**3. **L. W. Liebmann, S. M. Mansfield, A. K. Wong, M. A. Lavin, W. C. Leipold, and T. G. Dunham, “TCAD development for lithography resolution enhancement,” IBM J. Res. Develop **45(5)**, 651–665 (2001).
[CrossRef]

**4. **L. Pang, Y. Liu, and D. Abrams, “Inverse lithography technology (ILT): a natural solution for model-based SRAF at 45nm and 32nm,” Proc. SPIE **6607**, 660739 (2007).
[CrossRef]

**5. **S. H. Chan, A. K. Wong, and E. Y. Lam, “Initialization for robust inverse synthesis of phase-shifting masks in optical projection lithography,” Opt. Express **16(19)**, 14746–14761(2008).
[CrossRef]

**6. **E. Y. Lam and A. K. Wong, “Computation lithography: virtual reality and virtual virtuality,” Opt. Express **17(15)**, 12259–12268 (2009).
[CrossRef]

**7. **A. K. Wong and E. Y. Lam, “The nebulous hotspot and algorithm variability,” Proc. SPIE **7275**, 727509 (2009).
[CrossRef]

**8. **Y. Liu and A. Zakhor, “Optimal binary image design for optical lithography,” in *Proc. SPIE* **1264**, 401–412 (1990).

**9. **Y. Liu and A. Zakhor, “Binary and phase-shifting image design for optical lithography,” Proc. SPIE **1463**, 382–399 (1991).
[CrossRef]

**10. **S. Sherif, B. Saleh, and R. De Leone, “Binary image synthesis using mixed linear integer programming,” IEEE Trans. Image Process. **4(9)**, 1252–1257 (1995).
[CrossRef]

**11. **Y. C. Pati and T. Kailath, “Phase-shifting masks for microlithography: automated design and mask requirements,” J. Opt. Soc. Am. A **11(9)**, 2438–2452 (1994).
[CrossRef]

**12. **Y. Granik, “Solving inverse problems of optical microlithography,” *Proc. SPIE* **5754**, 506–526 (2004).
[CrossRef]

**13. **Y. Granik, K. Sakajiri, and S. Shang, “On objectives and algorithms of inverse methods in microlithography,” *Proc. SPIE* **6349**, 63494R (2006).
[CrossRef]

**14. **A. Poonawala and P. Milanfar, “Prewarping techniques in imaging: applications in nanotechnology and biotechnology,” Proc. SPIE **5674**, 114–127 (2005).
[CrossRef]

**15. **A. Poonawala and P. Milanfar, “Mask design for optical microlithography: an inverse imaging problem,” IEEE Trans. Image Process. **16(3)**, 774–788 (2007).
[CrossRef]

**16. **S. H. Chan and E. Y. Lam, “Inverse image problem of designing phase shifting masks in optical lithography,” in *Proceedings of IEEE International Conference on Image Processing*, pp. 1832–1835 (2008).

**17. **X. Ma and G. R. Arce, “Generalized inverse lithography methods for phase-shifting mask design,” Opt. Express **15(23)**, 15066–15079 (2007).
[CrossRef]

**18. **X. Ma and G. R. Arce, “PSM design for inverse lithography with partially coherent illumination,” Opt. Express **16(24)**, 20126–20141 (2008).
[CrossRef]

**19. **V. Singh, B. Hu, K. Toh, S. Bollepalli, S. Wagner, and Y. Borodovsky, “Making a trillion pixels dance,” *Proc. SPIE* **6924**, 69240S (2008).
[CrossRef]

**20. **A. Poonawala, Y. Borodovsky, and P. Milanfar, “ILT for double exposure lithography with conventional and novel materials,” *Proc. SPIE* **6520**, 65202Q (2007).
[CrossRef]

**21. **N. Jia, A. K. Wong, and E. Y. Lam, “Regularization of inverse photomask synthesis to enhance manufacturability,” *Proc. SPIE* **7520**, 752032 (2009).

**22. **N. Jia, A. K. Wong, and E. Y. Lam, “Robust photomask design with defocus variation using inverse synthesis,” *Proc. SPIE* **7140**, 71401W (2008).
[CrossRef]

**23. **L. Pang, G. Dai, T. Cecil, T. Dam, Y. Cui, P. Hu, D. Chen, K. Baik, and D. Peng, “Validation of inverse lithography technology (ILT) and its adaptive SRAF at advanced technology nodes,” *Proc. SPIE* **6924**, 69240T (2008).
[CrossRef]

**24. **S. Osher and R. P. Fedkiw, “Level set methods: an overview and some recent results,” J. Comput. Phys. **169(2)**, 463–502 (2001).
[CrossRef]

**25. **D. Adalsteinsson and J. A. Sethian, “A unified level set approach to etching, deposition and lithography I: algorithms and two-dimensional simulations,” J. Comput. Phys. **120(1)**, 128–144 (1995).
[CrossRef]

**26. **D. Adalsteinsson and J. A. Sethian, “A unified level set approach to etching, deposition and lithography II: three-dimensional simulations,” J. Comput. Phys. **122(2)**, 348–366 (1995).
[CrossRef]

**27. **D. Adalsteinsson and J. A. Sethian, “A unified level set approach to etching, deposition and lithography III: complex simulations and multiple effects,” J. Comput. Phys. **138(1)**, 193–223 (1997).
[CrossRef]

**28. **J. A. Sethian and D. Adalsteinsson, “An overview of level set methods for etching, deposition, and lithography development,” IEEE Trans. Semicond. Manuf. **10**, 167–184 (1997).
[CrossRef]

**29. **S. Osher and N. Paragios, *Geometric Level Set Methods in Imaging, Vision, and Graphics* (Springer Verlag New York, NJ, USA, 2003).

**30. **F. Santosa, “A level-set approach for inverse problems involving obstacles,” ESAIM Contröle Optim. Calc. Var. **1**, 17–33 (1996).
[CrossRef]

**31. **S. Osher and F. Santosa, “Level set methods for optimization problems involving geometry and constraints I. Frequencies of a two-density inhomogeneous drum,” J. Comput. Phys. **171(1)**, 272–288 (2001).
[CrossRef]

**32. **A. Marquina and S. Osher, “Explicit algorithms for a new time dependent model based on level set motion for nonlinear deblurring and noise removal,” SIAM J. Sci. Comp. **22**, 387–405 (2000).
[CrossRef]

**33. **A. K.-K. Wong, *Optical Imaging in Projection Microlithography* (SPIE Press, Bellingham, WA, 2005).

**34. **T. Chan, S. Esedoglu, F. Park, and A. Yip, “Recent developments in total variation image restoration,” in *Hand-book of Mathematical Models of Computer Vision*, pp. 17–32 (Springer Verlag, 2005).

**35. **Y. Shen, N. Wong, and E. Y. Lam, “Interconnect thermal simulation with higher order spatial accuracy,” in *Proceedings of IEEE Asia Pacific Conference on Circuits and Systems*, pp. 566–569 (2008).

**36. **S. Osher and R. Fedkiw, *Level Set Methods and Dynamic Implicit Surfaces* (Springer, 2002).

**37. **A. Harten, B. Engquist, S. Osher, and S. Chakravarthy, “Uniformly high order accurate essentially non-oscillatory schemes III,” J. Comput. Phys. **71**, 231–303 (1987).
[CrossRef]

**38. **C. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes,” J. Comput. Phys. **77(2)**, 439–471 (1988).
[CrossRef]

**39. **C. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes II,” J. Comput. Phys. **83(1)**, 32–78 (1989).
[CrossRef]

**40. **M. Minoux, *Mathematical Programming: Theory and Algorithms* (Wiley, New York, 1986).