Convolution is commonly performed on signals in many contexts, including the fields of sound, still image, video, lithography, and radio (radar) signal processing. Typically, the signals to be convolved are pattern signals. Each of the expressions “pattern” and “pattern signal” is used herein in a broad sense to denote a one-dimensional sequence or two-dimensional (or higher dimensional) array of data words (which can be, but need not be pixels). Typically, the data words comprise binary bits, and the convolution is performed in discrete fashion on the binary bits using software, digital signal processing circuitry, custom hardware, or FPGA systems (field programmable gate array based computing systems).
The term “data” herein denotes one or more signals indicative of data, and the expression “data word” herein denotes one or more signals indicative of a data word.
The motivations for implementing convolution rapidly, even when processing data indicative of very large patterns, are myriad. The present invention was motivated by the need for proximity correction in the field of lithography. In such problems, one attempts a two-dimensional convolution between data indicative of a large pattern “p” (where the pattern is a pixel array) and a diffusion kernel “d”. Often the kernel “d” is a Gaussian or a superposition of Gaussians, or is otherwise a smooth kernel. More specifically, the present invention grew out of attempts to establish a suitable “O(NN′)” algorithm (an algorithm requiring not more than on the order of NN′ multiplications and additions) for convolving a two-dimensional pattern comprising NN′ pixels, where each of N and N′ is very large) with a Gaussian kernel (or other smooth kernel) such that the convolution is exact or very close to exact.
The objective in performing proximity correction (in the field of lithography) is to generate a “raw” optical signal (or “raw” electron beam signal) which can be input to a set of reflective or refractive optics (or electron beam optics), in order to cause the output of the optics to produce a desired pattern on a mask or wafer. To determine the characteristics of a raw optical signal (or raw electron beam signal) that are needed to produce the desired pattern on the mask or wafer, a deconvolution operation is typically performed on a very large array of pixels (which determine a pattern “p”) in order to correct for the well known proximity problem. In the case of electron beam lithography, the proximity problem results from electron scattering in the substrate (mask or wafer) being written. Such scattering exposes broadened areas on the substrate to electrons (i.e., an area surrounding each pixel to be written in addition to the pixel itself), with the scattering effectively broadening the electron beam beyond the beam diameter with which the beam is incident on the substrate.
In nearly all proximity correction schemes, such a deconvolution operation includes at least one convolution step. Accordingly, in performing typical proximity correction, a very large array of pixels (determining a pattern “p”) must be convolved with a diffusion kernel. Although such a convolution is typically performed on a pattern comprising a very large array of binary pixels, this restriction is not essential in the following discussion and is not essential to implementation of the invention. Indeed, the invention can implement convolution on data indicative of any pattern “p” with a smooth convolution kernel “d” having characteristics to be described below.
For data indicative of a pattern “p” and a convolution kernel “d” we consider the cyclic convolution:                     (                  d          ⁢                      ×            C                    ⁢          p                )            n        =                  ∑                              i            +            j                    ≡                      n            ⁢                          (                              mod                ⁢                                                                   ⁢                N                            )                                          ⁢                        d          i                ⁢                  p          j                      ,where ×C denotes that the convolution operator has cyclic character, and an acyclic convolution which differs only in the indicial constraint and range:                     (                  d          ⁢                      ×            A                    ⁢          p                )            n        =                  ∑                              i            +            j                    =          n                    ⁢                        d          i                ⁢                  p          j                      ,where “×A” denotes that convolution operator has acyclic character.
For simplicity, we restrict much of the discussion herein to one-dimensional cases (in which the pattern p is an ordered set of N data values and the kernel is an ordered set of M values). Despite this, it should be appreciated that in typical embodiments of the invention, the pattern is two-dimensional (a two-dimensional array pjk of data values determines the pattern) and the summation defining the convolution (a summation which corresponds to either one of the summations set forth in the previous paragraph) is over index k as well as index j of the array pjk. In the case of a two-dimensional pattern p determined by an N by N′ array of data values, the indices n, i, j and domain lengths N in the formulae set forth in the previous paragraph are 2-vectors.
In one-dimensional cases, the result of the cyclic convolution has length N (it comprises N data values), and the result of the acyclic convolution has length M+N−1.
It is standard that a cyclic convolution d×Cp can be cast as an equivalent matrix-vector product:d×Cp≡Dp,where D is the circulant matrix of d (hereinafter the “circulant” of d), whose 1-dimensional form is defined (assuming that N is greater than 3) as:   D  =            (                                                  d              0                                                          d              1                                                          d              2                                                          d              3                                            ⋯                                              d                              N                -                1                                                                                        d                              N                -                1                                                                        d              0                                                          d              1                                                          d              2                                            ⋯                                              d                              N                -                2                                                                          ⋮                                                                                                                                                                                                                                                                            ⋮                                                              d              1                                                          d              2                                                          d              3                                                          d              4                                            ⋯                                              d              0                                          )        .  Therefore, conventional methods for cyclic convolution can be cast in the language of matrix algebra. Acyclic convolution can be obtained with similar matrix manipulations.
For the sake of simplicity, we will use the symbol × hereinbelow to denote convolution having either acyclic or cyclic character. In most of the discussion the symbol will refer to convolution having acyclic character. Those of ordinary skill in the art will recognize that, given a specified acyclic convolution, a corresponding cyclic convolution can be implemented by slight modification of the parameters (e.g., the boundary conditions and definition of the circulant of the kernel) that determine the acyclic convolution.
Above-referenced U.S. patent application Ser. No. 09/480,908 discloses a fast convolution method whose central idea (in one-dimensional embodiments) is to approximate a smooth kernel d by a polynomial spline kernel ƒ (where ƒ is a spline function ƒ(x) which is piecewise a polynomial of degree δ with L pieces f1(x)), and then to use appropriate operators that annihilate (or flatten) each polynomial of given degree (in a manner to be explained) to calculate the convolution of ƒ and p quickly. In some embodiments, the smooth kernel d is approximated by a spline kernel ƒ which is not a polynomial spline kernel, but which consists of L pieces defined over adjacent segments of its domain (in typical two-dimensional cases, the latter spline kernel is a radially symmetric function whose domain is some continuous or discrete set of values of the radial parameter). Though “spline” convolution as described in U.S. application Ser. No. 09/480,908 has features reminiscent of conventional wavelet schemes and is an O(N) algorithm (as are wavelet schemes), an advantage of “spline” convolution is that it can be performed (on data indicative of a pattern p consisting of N data values) with cN arithmetic operations (multiplications and additions), whereas conventional wavelet convolution on the same data would require bN arithmetic operations, where the factor “b” is typically (i.e., with typical error analysis) significantly larger than the factor “c.” In other words, the implied big-O constant for the spline convolution is significantly smaller than the typical such constant for conventional wavelet convolution.
Spline convolution, as described in U.S. application Ser. No. 09/480,908, is a method for performing cyclic or acyclic convolution of a pattern “p” (i.e., data indicative of a pattern “p”) with a smooth diffusion kernel d, to generate data indicative of the convolution result r=Dp, where D is the circulant of d. The pattern “p” can be one-dimensional in the sense that it is determined by a continuous (or discrete) one-dimensional domain of data values (e.g., pixels), or it can be two-dimensional in the sense that it is determined by a continuous two-dimensional domain of data values (or a two-dimensional array of discrete data values), or p can have dimension greater than two. In typical discrete implementations, the pattern p is one-dimensional in the sense that it is determined by a discrete, ordered set of data values (e.g., pixels) pi, where i varies from 0 to N−1 (where N is the signal length), or it is two-dimensional in the sense that it is determined by an array of data values pij, where i varies from 0 to N−1 and j varies from 0 to N′−1, or it has dimension greater than two (it is determined by a three or higher-dimensional set of data values). Typically, the kernel d is determined by an array of data values dij, where i varies from 0 to N−1 and j varies from 0 to N′−1 (but the kernel d can alternatively be determined by a discrete set of data values d0 through dN−1).
In some embodiments described in U.S. application Ser. No. 09/480,908, the convolution Dp is accomplished by performing the steps of:
(a) approximating the kernel d by a polynomial spline kernel ƒ (unless the kernel d is itself a polynomial spline kernel, in which case d=ƒ and step (a) is omitted);
(b) calculating q=Bp=Δδ+1Fp, where F is the circulant of kernel ƒ, and Δδ+1 is an annihilation operator (whose form generally depends on the degree δ of the polynomial segments of ƒ) which operates on circulant F in such a manner that Δδ+1F=B is sparse; and
(c) backsolving Δδ+1r=q to determiner=Fp.
In cases in which the kernel d is itself a polynomial spline kernel (so that d=ƒ and F=D), the method yields an exact result (r=Dp). Otherwise, the error inherent in the method is (ƒ−d)×p, where × denotes convolution, and thus the error is bounded easily.
In one-dimensional cases (in which the pattern to be convolved is a one-dimensional pattern of length N), Δδ+1 has the form of the N×N circulant matrix defined as follows:       Δ          δ      +      1        =      [                                        +                          (                                                                                          δ                      +                      1                                                                                                            0                                                              )                                                            -                          (                                                                                          δ                      +                      1                                                                                                            1                                                              )                                                            +                          (                                                                                          δ                      +                      1                                                                                                            2                                                              )                                                            -                          (                                                                                          δ                      +                      1                                                                                                            3                                                              )                                                ⋯                          0                                      0                                      +                          (                                                                                          δ                      +                      1                                                                                                            0                                                              )                                                            -                          (                                                                                          δ                      +                      1                                                                                                            1                                                              )                                                            +                          (                                                                                          δ                      +                      1                                                                                                            2                                                              )                                                ⋯                          0                                      ⋮                                                                                                                                                                                                                              ⋮                                                  -                          (                                                                                          δ                      +                      1                                                                                                            1                                                              )                                                            +                          (                                                                                          δ                      +                      1                                                                                                            2                                                              )                                                ⋯                                                                           0                                      +                          (                                                                                          δ                      +                      1                                                                                                            0                                                              )                                            ]  in which each entry is a binomial coefficient, and δ is the maximum degree of the spline segments of spline kernel ƒ. For example, δ=2 where the spline kernel ƒ comprises quadratic segments. In two- or higher-dimensional cases, the annihilation operators can be defined as       Δ    =                            ∂                      x            1                                n            1                          ⁢                              ∂                          x              2                                      n              2                                ⁢                                           ⁢          …                    ⁢                           ⁢              ∂                  x          d                          n          d                      ,where ∂xhnh is the n-th partial derivative with respect to the h-th of d coordinates. For example, the Laplacian       ∇    2    ⁢      =                            ∂                      x            1                    2                ⁢                  +                                           ⁢          …                    ⁢                           +              ∂                  x          d                2            will annihilate piecewise-planar functions of d-dimensional arguments ƒ=ƒ(x1, x2, . . . xd).
In the one-dimensional case, the end points of each segment (the “pivot points”) of spline kernel ƒ may be consecutive elements di and di+1 of kernel d, and step (a) can be implemented by performing curve fitting to select each segment of the spline kernel as one which adequately matches a corresponding segment of the kernel d. In some implementations, appropriate boundary conditions are satisfied at each pivot point, such as by derivative-matching or satisfying some other smoothness criterion at the pivot points.
In some implementations described in application Ser. No. 09/480,908, step (c) includes a preliminary “ignition” step in which a small number of the lowest components of r=Fp are computed by exact multiplication of p by a few rows of F, and then a step of determining the rest of the components of r using a natural recurrence relation determined by the spline kernel and the operator Δδ+1. For example, in the one-dimensional case, the lowest components of r are r0, r1, . . . , rδ, where “δ” is the maximum degree of the spline segments of spline kernel ƒ (for example r0, r1, and r2 where the spline kernel comprises quadratic segments), and these (δ+1) components are determined by exact multiplication of p by (δ+1) rows of F. The (δ+1) components can alternatively be determined in other ways. Then, the rest of the components “rδ” are determined using a natural recurrence relation determined by the operator Δδ+1. The “ignition” operation which generates the components r0, r1, . . . , rδ, can be accomplished with O(N) computations. The recurrence relation calculation can also be accomplished with O(N) computations.
In other embodiments, the method disclosed in U.S. application Ser. No. 09/480,908 for performing the convolution r=Dp includes the steps of:
(a) approximating the kernel d by a polynomial spline kernel ƒ (unless the kernel d is itself a polynomial spline kernel, in which case d=ƒ and step (a) is omitted);
(b) calculating q=Bp=ΔδFp, where F is the circulant of kernel ƒ and Δδ is a flattening operator (whose form generally depends on the degree δ of the polynomial segments of F, and which operates on circulant F such that B=ΔδF is almost everywhere a locally constant matrix); and
(c) backsolving Δδr=q to determine r=Fp. In one-dimensional cases (in which p has length N), Δδ has the form of the N×N circulant matrix:       Δ    δ    =      [                                        +                          (                                                                    δ                                                                                        0                                                              )                                                            -                          (                                                                    δ                                                                                        1                                                              )                                                            +                          (                                                                    δ                                                                                        2                                                              )                                                            -                          (                                                                    δ                                                                                        3                                                              )                                                ⋯                          0                                      0                                      +                          (                                                                    δ                                                                                        0                                                              )                                                            -                          (                                                                    δ                                                                                        1                                                              )                                                            +                          (                                                                    δ                                                                                        2                                                              )                                                ⋯                          0                                      ⋮                                                                                                                                                                                                                              ⋮                                                  -                          (                                                                    δ                                                                                        1                                                              )                                                            +                          (                                                                    δ                                                                                        2                                                              )                                                ⋯                                                                           0                                      +                          (                                                                    δ                                                                                        0                                                              )                                            ]  in which each entry is a binomial coefficient, and δ is the maximum degree of the spline segments of spline kernel ƒ. In higher-dimensional cases, the flattening operator Δ67  is defined similarly.
In other embodiments disclosed in U.S. application Ser. No. 09/480,908, the convolution Dp (where D is the circulant of smooth kernel d) includes the steps of:
(a) approximating the kernel d by a spline kernel ƒ which is not a polynomial spline kernel (unless the kernel d is itself such a spline kernel, other than a polynomial spline kernel, in which case d=ƒ and step (a) is omitted);
(b) calculating q=Bp=AFp, where F is the circulant of kernel ƒ and A is an annihilation or flattening operator, where A operates on circulant F in such a manner that AF=B is sparse when A is an annihilation operator, and A operates on circulant F in such a manner that AF=B is almost everywhere a locally constant matrix when A is a flattening operator; and
(c) back-solving Ar=q to determine r=Fp.
To better appreciate the advantages of the present invention over conventional convolution methods, we next explain two types of conventional convolution methods: Fourier-based convolution and wavelet-based convolution.
As is well known, Fourier-based convolution relies on the elegant fact that if F is a Fourier matrix, sayFjk=e−2πijk/Nthen the transformation FDF−1 of the circulant is diagonal, whence we compute:Dp=F−1(FDF−1)Fp,where the far-right operation Fp is the usual Fourier transform, the operation by the parenthetical part is (by virtue of diagonality) dyadic multiplication, and the final operation F−1 is the inverse Fourier transform. For arbitrary D one requires actually three Fourier transforms, because the creation of the diagonal matrix FDF−1 requires one transform. However, if D is fixed, and transformed on a one-time basis, then subsequent convolutions Dp only require two transforms each, as is well known. The complexity then of Fourier-based cyclic convolution is thus O(N log N) operations (i.e., on the order of N log N multiplications and additions) for convolving a pattern p of length N (a pattern determined by N data values), because of the 2 or 3 FFTs (Fast Fourier Transforms) required. It should be noted that the Fourier method is an exact method (up to round-off errors depending on the FFT precision).
Another class of conventional convolution methods consists of wavelet convolution methods, which, by their nature, are generally inexact. The idea underlying such methods is elegant and runs as follows in the matrix-algebraic paradigm. Assume that, given an N-by-N circulant D, it is possible to find a matrix W (this is typically a compact wavelet transform) which has the properties:    (1) W is unitary (i.e. W−1 is the adjoint of W);    (2) W is sparse; and    (3) WDW−1 is sparse,where “sparse” in the present context denotes simply that any matrix-vector product Wx, for arbitrary x, involves reduced complexity O(N), rather than say O(N2).
With the assumed properties, we can calculate:Dp=W−1(WDW−1)Wpby way of three sparse-matrix-vector multiplications, noting that unitarity implies the sparseness of W−1. Therefore the wavelet-based convolution complexity is O(N) for convolving a pattern p determined by N data values, except that it is generally impossible to find, for given circulant D, a matrix W that gives both sparsity properties rigorously. Typically, if the convolution kernel d is sufficiently smooth, then a wavelet operator W (which is sparse) an be found such that within some acceptable approximation error the property (3) above holds. Above-noted properties (1) and (2) are common at least for the family of compact wavelets (it is property (3) that is usually approximate).
An advantage of “spline” convolution (in accordance with the teaching of U.S. application Ser. No. 09/480,908) over conventional wavelet convolution is that it can be performed (on data indicative of a pattern p comprising N data values) with cN arithmetic operations, whereas conventional wavelet convolution on the same data would require bN arithmetic operations, where (assuming typical error budgets) the factor “b” is significantly larger than the factor “c.” Among the other important advantages of the “spline” convolution method of application Ser. No. 09/480,908 (over conventional convolution methods) are the following: spline convolution is exact with respect to the spline kernel f, whereas wavelet convolution schemes are approximate by design (and error analysis for wavelet convolution is difficult to implement); the signal lengths for signals to be convolved by spline convolution are unrestricted (i.e., they need not be powers of two as in some conventional methods, and indeed they need not have any special form); and spline convolution allows acyclic convolution without padding with zeroes.
Separated-spline convolution in accordance with the present invention (like spline convolution in accordance with U.S. application Ser. No. 09/480,908) is an O(N) method for convolving a pattern p determined by N data values. Separated-spline convolution in accordance with the present invention has an advantage over spline convolution in accordance with U.S. application Ser. No. 09/480,908 in that separated-spline convolution in accordance with the invention can be performed (on data indicative of a two- or higher-dimensional patterns consisting of N data values) with dN arithmetic operations (multiplications and additions), whereas spline convolution in accordance with U.S. application Ser. No. 09/480,908 on the same data would require cN arithmetic operations, where the factor “c” is larger, and typically significantly larger, than the factor “d.” In other words, the implied big-O constant for separated-spline convolution according to the present invention is significantly smaller than the typical implied big-O constant for spline convolution as described in U.S. patent application Ser. No. 09/480,908.