Level-set-based inverse lithography for photomask synthesis

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 OCIS codes: (110.3960) Microlithography; (110.5220) Photolithography; (110.1758) Computational imaging. 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). 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). 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). 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–14760 (2008). 6. E. Y. Lam and A. K. Wong, “Computation lithography: virtual reality and virtual virtuality,” Opt. Express 17(15), 12259–12268 (2009). 7. A. K. Wong and E. Y. Lam, “The nebulous hotspot and algorithm variability,” Proc. SPIE 7275, 727509 (2009). 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). 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). 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). 12. Y. Granik, “Solving inverse problems of optical microlithography,” Proc. SPIE 5754, 506–526 (2004). 13. Y. Granik, K. Sakajiri, and S. Shang, “On objectives and algorithms of inverse methods in microlithography,” Proc. SPIE 6349, 63494R (2006). 14. A. Poonawala and P. Milanfar, “Prewarping techniques in imaging: applications in nanotechnology and biotechnology,” Proc. SPIE 5674, 114–127 (2005). 15. A. Poonawala and P. Milanfar, “Mask design for optical microlithography: an inverse imaging problem,” IEEE Trans. Image Process. 16(3), 774–788 (2007). 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). (C) 2009 OSA 21 December 2009 / Vol. 17, No. 26 / OPTICS EXPRESS 23690 #119341 $15.00 USD Received 30 Oct 2009; revised 4 Dec 2009; accepted 4 Dec 2009; published 10 Dec 2009 17. X. Ma and G. R. Arce, “Generalized inverse lithography methods for phase-shifting mask design,” Opt. Express 15(23), 15066–15079 (2007). 18. X. Ma and G. R. Arce, “PSM design for inverse lithography with partially coherent illumination,” Opt. Express 16(24), 20126–20141 (2008). 19. V. Singh, B. Hu, K. Toh, S. Bollepalli, S. Wagner, and Y. Borodovsky, “Making a trillion pixels dance,” Proc. SPIE 6924, 69240S (2008). 20. A. Poonawala, Y. Borodovsky, and P. Milanfar, “ILT for double exposure lithography with conventional and novel materials,” Proc. SPIE 6520, 65202Q (2007). 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). 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). 24. S. Osher and R. P. Fedkiw, “Level set methods: an overview and some recent results,” J. Comput. Phys. 169(2), 463–502 (2001). 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). 26. D. Adalsteinsson and J. A. Sethian, “A unified level set approach to etching, deposition and lithography II: threedimensional simulations,” J. Comput. Phys. 122(2), 348–366 (1995). 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). 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). 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). 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). 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). 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 Handbook 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). 38. C. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes,” J. Comput. Phys. 77(2), 439–471 (1988). 39. C. Shu and S. Osher, “Efficient implementation of essentially non-oscillatory shock-capturing schemes II,” J. Comput. Phys. 83(1), 32–78 (1989). 40. M. Minoux, Mathematical Programming: Theory and Algorithms (Wiley, New York, 1986).

Frequencies of a two-density inhomogeneous drum," J. Comput

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][2][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 timeconsuming, 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][17][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.

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.

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 T {·} 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Û(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 T {·} 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] sig(U(x)) = 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.

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 R. This functional takes as its argument an image U and returns a scalar that quantifies the quality of the pattern. Smaller values of R correspond to more desirable images. R 1 , R 2 , and R 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): R 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.
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 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 (12) 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).

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 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 λ ∇ · ∇U |∇U| . Preconditioning Eq. (15) for speed improvement is also possible [29, 32].

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 T (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 T (U) at U, will produce a small variation δ φ that reduces F(U), which leads to ∂ φ ∂t = −|∇φ |α(x,t).  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 O (N 2 log N), assuming an N × N image.

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 firstorder 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.

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) Fig. 1. Simulation of lithographic imaging with different mask patterns computed using different spatial schemes. The first column denotes the input U(x), the second column I aerial (x), and the third column I(x). Row (a) uses the target circuit pattern I 0 as input. Row (b) uses the pattern derived by our level set algorithm, with first-order temporal accuracy and ENO1 spatial accuracy. Row (c), (d) and (e) are similar to (b), but with ENO2, ENO3 and WENO spatial accuracy, respectively.  Fig. 1(c) using ENO2) and Pattern Error (pixel difference) in Fig. 1 Experiment Fig. 1(a) Fig. 1(b) Fig. 1(c) Fig. 1(d) Fig. 1 Fig. 2. Simulation of lithographic imaging with different mask patterns computed using different temporal schemes. The first column denotes the input U(x), the second column I aerial (x), and the third column I(x). Row (a) uses the target circuit pattern I 0 as input. Row (b) uses the pattern derived by our level set algorithm, with ENO2 spatial and first-order temporal accuracy . Row (c) and (d) are similar to (b), but with second-order and third-order temporal accuracy, respectively.  Fig. 2(b) using first-order temporal accuracy) and Pattern Error (pixel difference) in Fig. 2 Experiment Fig. 2(a) Fig. 2(b) Fig. 2(c) Fig. 2 Fig. 3(b) and Fig. 3(d) using gradient method normalized against that in Fig. 3(c) and Fig. 1(c) using the proposed method, respectively) and Pattern Error (pixel difference) in Fig. 3 Experiment Fig. 3(a) Fig. 3(b) Fig. 3(c) Fig. 3 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.

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.