1. Field of the Invention
The present invention relates to a method for image de-blurring, based on estimating an image's Lipschitz exponent a using a direct computational technology rather than an iterative technique. More particularly, the method of the present invention is related to de-blurring of an image based on singular integrals and Fast Fourier Transform (FFT) algorithms.
2. Description of the Related Art
Most images fix, y) are not differentiable functions of the space variables x and y. Rather, they exhibit edges, singularities, localized sharp features, and various other kinds of important fine-scale details or texture. For such non-smooth imagery, prior art has generally formulated the ill-posed image-deblurring problem incorrectly. This often leads to flawed reconstructions, where vital small-scale information has been smoothed out, or where unexpected noise-induced graininess obscures fine detail.
Digital image acquisition plays an ever-increasing role in science, technology, and medicine, and image deblurring is becoming an increasingly important image processing activity. In addition, with the widespread use of digital cameras and camera cell-phones, there is growing general interest in the possibilities of post-processed digital image enhancement. In another direction, wavefront coding is a revolutionary new idea in Optics that calls for deliberately designing an imperfect lens, see, e.g., D. MacKenzie, Novel Imaging Systems Rely On Focus-Free Optics, SIAM News, Volume 36#6, July-August 2003, the contents of which are hereby incorporated by reference in their entirety. That lens produces a blurred image, but one where the depth of field has been significantly increased. Mathematical deconvolution is subsequently applied to the blurred image to remove the designed blur. This results in a superior photograph where distant and close-in objects are equally well focused. Carl Zeiss Inc., is reportedly set to manufacture such a lens. In these applications, fast computational throughput for large size imagery is very desirable. Some deconvolution methods involve computationally intensive nonlinear iterative procedures, typically requiring hours of computing time. Direct (non-iterative) deconvolution methods, that can process 1000×1000 images in less than a minute of computing time, are considered real-time methods, and are highly sought after.
Image deblurring is a difficult ill-posed mathematical problem, requiring for its correct solution prior knowledge and specification of the smoothness characteristics in the unknown exact sharp image f(x, y). However, most commonly occurring images are not differentiable functions of the space variables x and y. Rather, these images display edges, localized sharp features, and various other fine-scale details or texture. For such non-smooth imagery, prior art has generally stabilized the ill-posed deblurring problem by prescribing L2 bounds for the sharp image f(x, y). The L2-Tikhonov-Miller method, is the best-known example of that approach, see, e.g., K. Miller, Least Squares Methods For Ill-Posed Problems With A Prescribed Bound, SIAM J. Math. Anal., 1 (1970), pp 52-74; R. L. Lagendijk and J. Biemond, Iterative Identification and Restoration of Images, Kluwer Academic Publisher, Norwell, Mass., 1991, the contents of both of which are hereby incorporated in their entirety. In another direction, considerable research during the last ten years has been based on the assumption that images belong to BV(R2), the space of functions of bounded variation. This has led to nonlinear partial differential equation (PDE) deblurring procedures, where bounds are prescribed on the total variation or TV seminorm
The Marquina-Osher TV algorithm is one of the most widely used PDE deblurring methods, see, e.g., Marquina-Osher, Explicit Algorithms For A New Time Dependent Model Based On Level Set Motions For Nonlinear Deblurring And Noise Removal, SIAM J. Sci. Comput., 22 (2000), pp. 387-405, the contents of which are hereby incorporated in their entirety. However, each of these two general deblurring approaches is fundamentally flawed theoretically, and that flaw often translates into poor quality reconstructions. Thus, prescribing L2 bounds insufficiently constrains the Tikhonov-Miller solution, which is typically found to be contaminated by noise. Also, as was recently proved, most natural images are not of bounded variation, see Y. Gousseau and J. M. Morel, Are Natural Images of Bounded Variation?, SIAM J. Math Anal., 33 (2001), pp. 634-648, the contents of which are hereby incorporated by reference in their entirety. As a result, Marquina-Osher TV deblurring often leads to unacceptable loss of fine-scale information.
Correct characterization and calibration of the lack of smoothness of images is crucial in image deblurring, as well as in other image processing tasks. As functions of x and y, most images are significantly better behaved than the most general L2 functions, while being significantly less smooth than functions of bounded variation. For this reason, both the L2-Tikhonov-Miller and TV-Marquina-Osher methods are incorrectly formulated.
Thus, a satisfactory solution to the deblurring problem is needed. The present invention solves the deblurring problem by providing a method for singular integral deblurring of images. Very recently it has become apparent that so-called Lipschitz spaces Λ(α, p, q), are the correct framework for image deblurring, see A. S. Carasso, Singular Integrals, Image Smoothness, and The Recovery of Texture In Image Deblurring, NISTIR #7005, June 2003, National Institute of Standards and Technology, Gaithersburg, Md. 20899 and A. S. Carasso, Singular Integrals, Image Smoothness, and The Recovery of Texture In Image Deblurring, SIAM J. Appl. Math., 64(2004), pp. 1749-1774. Such spaces describe functions with fractional derivatives, and can easily accommodate non-smooth images. The Lp Lipschitz exponent α, 0<α<1, measures the fine-scale content of an image, provided that image is relatively noise free. Heavily textured images have low values for a, while a large value of a indicates that an image is relatively smooth. Estimating an image's Lipschitz exponent a is a delicate problem.
In the method of the present invention, image deblurring is treated as an entirely separate issue from image smoothness characterization. The present invention is based on a new and fundamental reformulation of the image deblurring problem in which the Lipschitz space characterization of the unknown desired sharp image f(x, y) has been explicitly incorporated. This Lipschitz space characterization pertains to the space Λ(α, 2, ≈), and is expressed in terms of the singular integral (SI) method of the present invention, which is fully developed and discussed in A Direct Procedure For Classifying Image Smoothness, Based On Singular Integral Operators and Fast Fourier Transform Algorithms, Aug. 15, 2003, which is hereby incorporated by reference in its entirety. An entirely new energy functional is provided for the deblurring problem that is used to define the deblurred image. In this functional, Fourier analysis and singular integral operators play essential roles in specifying and enforcing the required Lipschitz space information.
The present invention is directed to a method for solving the image deblurring problem where one is given a noisy blurred image g(x, y), and a known shift-invariant point spread function p(x, y). Each of these is assumed to be a 2J×2J array of non-negative numbers, with the array p(x, y) summing to unity. This problem has the mathematical form
p(x, y)f(x, y)=g(x, y),
where denotes convolution. We may write
g(x, y)=ge(x, y)+n(x, y), (2)
where ge(x, y) is the hypothetical exact blurred image that would have been recorded in the absence of any noise, and n(x, y) presumed small, represents the cumulative effects of all noise processes and other errors affecting final acquisition of the digitized array g(x, y). The noise may be multiplicative. Neither ge(x, y) nor n(x, y) are known, and only their sum g(x, y) is known. Denoting the unknown exact sharp image by fe(x, y), we have
p(x, y)fe(x, y)=ge(x, y), 0≦fe(x, y)≦MAX. (3)
Note that blurring is an averaging process that reduces maximum image intensity. Hence, we necessarily have 0≦ge(x, y)<MAX in Eq. (3), if MAX represents the actual maximum pixel value in the sharp image fe(x, y). The fundamental difficulty in deblurring is due to the fact that g(x, y) is given, not ge(x, y). Thus, one must solve Eq.(1), rather than Eq.(3). Because of the notorious ill-conditioned nature of that problem, special precautions are necessary to find a solution f(x, y) in Eq.(1) that is a good approximation to the correct solution fe(x, y) in Eq.(3). A-priori information about fe(x, y) and n(x, y) is an essential element in the successful solution of Eq.(1).
For any 2J×2J array i(x, y), define the discrete L2 norm of i(x, y) as follows
Denote the forward discrete Fourier transforms of the 2J×2J arrays p(x, y) and g(x, y), by {circumflex over (p)}(ξ, η) and ĝ(ξ,η), respectively, where ξ, η are integers and −J≦ξ, η≦J. These discrete Fourier transforms can be computed using Fast Fourier Transform (FFT) algorithms.
This invention is a direct non-iterative deblurring method, capable of restoring images of size 1024×1024 in seconds of CPU time. The method is therefore useable in real-time applications in such fields as Astronomy, Electron Microscopy, and Medicine. The prior art Lucy-Richardson, Van Cittert, Landweber, Maximum Entropy, and related methods, are examples of iterative methods, see, respectively, {L. B. Lucy, An Iterative Technique For The Rectification of Observed Distributions, Astronomical Journal, 79 (1974), pp. 745-754; R. L. Lagendijk and J. Biemond, Iterative Identification and Restoration of Images, Kluwer Academic Publishers, Norwell, Mass., 1991; L. Landweber, An Iteration Formula For Fredholm Integral Equations Of The First Kind, Am. J. Math., 73 (1951), pp. 615-624; and E. S. Meinel, Maximum Entropy Image Restoration: Lagrangian and Recursive Techniques, J. Opt. Soc. Amer. Series A, 5 (1988), pp. 25.29} the respective contents of which are hereby incorporated in their entirety. These methods typically require several thousand iterations to resolve fine detail, and necessitate several hours of CPU time for large size imagery. Such methods are not useable in real-time applications. Accordingly, iterative methods are not considered relevant prior art for the present invention, and are not discussed further. However, a survey and evaluation of the merits and shortcomings of these iterative methods is given in A. S. Carasso, Linear and Non-linear Image Deblurring: A Documented Study, SIAM J. Numer. Anal., 36(1999), pp. 1659-1689, the contents of which are hereby incorporated by reference in their entirety.
The prior art true Wiener filtering method requires a-priori knowledge of the exact power spectra of both the unknown desired sharp image fe(x, y) and the noise component n(x, y), i.e., the quantities |{circumflex over (f)}e(ξ,η)|, and |{circumflex over (n)}(ξ,η)|, for −J≦ξ, η≦J. This represents a total of 8J2 numbers. That information is almost never available in practice, and true Wiener filtering is unrealizable, except in controlled numerical experiments. However, the method is of considerable theoretical interest because it provides the mathematically optimal deblurred image, i.e., that image whose departure from fe(x, y) is the smallest in the L2 norm. Examples of True Wiener filtering are included in each of
The prior art Tikhonov-Miller method requires a-priori knowledge of positive constants E and M, such that
∥n∥2≦ε, ∥f∥2≦M, ε<<M. (5)
The Tikhonov-Miller method is considered to be a realizable approximate version of the unrealizable true Wiener method, requiring only two a-priori numbers, versus 8J2. The resulting deblurred image may be far from optimal, however. In fact, the quality of the Tikhonov-Miller image is usually much inferior to that of the true Wiener image. An example of Tikhonov-Miller deblurring is included in
The TV method is currently very much in vogue. The method is computationally intensive, and is typically implemented by numerically solving a nonlinear partial differential equation of anisotropic diffusion type. The great strength of this method is its ability to filter out noise while preserving edges. However, a very major disadvantage of the TV method is its tendency to eliminate vital small-scale information, unbeknownst to the user. In fact, the method is based on a mathematically false premise. Most natural images are not of bounded variation, and the true value of the constant KTV in Eq. (6) is infinite. Any implementation of TV deblurring must use a finite value for KTV. This can lead to unacceptable smoothing out of texture and other fine detail. Examples of TV deblurring are included in
The Singular Integral or SI method of the present invention is a direct, non-iterative method requiring only four a-priori numbers. Nevertheless, the SI method consistently delivers deblurred imagery that is quite close in quality to that found in unrealizable, true Wiener filtering. In particular, the SI method avoids the noise-induced graininess that is typical in the prior art Tikhonov-Miller method, as well as the unacceptable loss of texture that is generally found in the prior art Total Variation (TV) method.
No role is played in the method of the present invention by any of multi-resolution analysis, wavelet theory, wavelet expansions, decay of wavelet coefficients, image compression, or image quantization, and no knowledge whatsoever of any of these topics is required to implement the present approach computationally. In fact, very modest computational effort is actually necessary to implement the Lipschitz space deblurring method of the present invention that is based on new energy functional π(f, λ). The desired deblurred images is the unique function fSI(x, y) that minimizes π(f, λ) over all f∈L2. Commonly available FFT algorithms can be used to construct fSI(x, y) directly and this results in the fast direct deblurring method of the present invention. Note that use of the energy functional π(f, λ). requires specifying a value of λ, where 0<λ≦1. Each choice of λ translates into a specific choice of singular integral semigroup Sλt. The choice of singular integral operator becomes the one that is used both for classifying image smoothness, and for regularizing the ill-posed deblurring problem. Therefore, each choice of λ defines an instantiation of the functional π(f, λ) as the concrete, correctly formulated, Lipschitz space image deblurring method of the present invention. Because of the multiplicity of distinct values of λ that can be used, an entire class of image deblurring methods is actualy subsumed under the umbrella of the functional π(f, λ), through the family Sλt. The present invention therefore includes the infinitely many distinct methods defined by the family Sλt, 0<λ≦1.0.
In experiments conducted using the method of the present invention, the particular choice λ=0.5 was made, corresponding to the Poisson Singular Integral. However, similar results can be obtained with other choices of λ.
It is to be understood by persons of ordinary skill in the art that the following descriptions are provided for purposes of illustration and not for limitation. An artisan understands that there are many variations that lie within the spirit of the invention and the scope of the appended claims. Unnecessary detail of known functions and operations may be omitted from the current description so as not to obscure the present invention.
The present invention is directed to a two-part image deblurring method having a classification module for classifying image smoothness and an image deblurring module, that are both based on singular integral operators and Fast Fourier Transform algorithms. Given a square image f(x, y), the present invention includes a first method for determining a Lipschitz exponent a of the square image, for a particular user-selected LP norm. This determination implies that the given image belongs to a specific Lipschitz space Λ(α, p, ∞). In addition, for a specific user-selected choice of singular integral operator Sλt, the image classifying method determines a positive constant C such that
∥Sλtf−f∥p≦C∥f∥ptαa/2λ0<t≦0.1. (7)
When p=2, the information contained in Eq.(7) is the basis for a singular integral image deblurring method associated with the operator Sλt.
The first module is a classification module comprising the steps of:
1. Given a noise-free 2J×2J image f(x, y), applying commonly available Fast Fourier Transform Algorithms to form the forward Fourier Transform of f(x, y), denoted by {circumflex over (f)}(ξ, η), where ξ and η are integers and −J≦ξ, η≦J.
2. For a preselected fixed value of λ, where 0.5≦λ≦1, and a preselected fixed finite sequence {tn}n=1N tending to zero, forming the function ĝλ(ξ, η, tn) for each tn, where
ĝλ(ξ, η, tn)=e−t
The Poisson singular integral, denoted by Ut, corresponds to the choice λ=0.5 in Eq.(8), while the Gaussian singular integral, denoted by Gt, corresponds to λ=1.0. Other values of λ correspond to more general isotropic Levy stable singular integrals denoted by Stλ. Values of λ<0.5 are generally less useful in this method. For a wide class of 512×512 and 1024×1024 8-bit images, the sequence tn=(0.95)n, n=1, . . . , 300, has been found adequate.
3. For each tn, using FFT algorithms to form the inverse Fourier transform of ĝλ(ξ, η, tn) denoted by gλ(x, y, tn). In the function gλ(x, y, tn), all negative values are reset to the value zero. It is not necessary to display the image corresponding to gλ(x, y, tn). However, that image would be a blurred version of the original image f(x, y). For each tn let hn(x, y) denote the difference between f(x, y) and gλ(x, y, tn),
hn(x, y)=f(x, y)−gλ(x, y, tn) (9)
4. For any 2J×2J image i(x, y), and any integer p≧1, defining the discrete Lp norm of i(x, y), as follows
where the values p=1 and p=2 are of particular interest.
5. For a predetermined value of p, calculating the positive number ∥f∥p, for the given image f(x, y).
6. For each tn, calculating the finite sequence {μ(tn)}n=1N where and where hn(x, y) is defined in Eq.(9). Note that steps 2,-6 can be performed simultaneously, for numerous distinct tn values, on multiprocessor computers. Such parallel computation of μ(tn) becomes advantageous for large size imagery.
7. Plotting in μ(tn) v. ln (tn) for n=1, . . . , N.
An example of the resulting curve is illustrated in
In
In
In fact, an infinite variety of singular kernels can be applied to estimate image Lipschitz exponents a. General radially symmetric Levy stable singular operators Sλt, include Gaussian kernels (λ=1.0), and Poisson kernels (λ=0.5), as special cases. Referring now to
8. Characteristic of the solid curves in
9. Using commonly available Linear Least Squares algorithms to find the straight lines Γ and Σ, as shown in
10. For the preselected specific value of p, the returned values for σ, b, and δ in step 9, yield the Lipschitz exponent α through
α=2σλ (12)
as well as the important inequality
∥Sλtf−f∥p≦C∥f∥ptσ, 0<t≦0.1, (13)
where
C=exp(b+δ). (14)
This implies that f(x, y) belongs to Λ(α, p, ≈). When p=2, the values of C and σ in Eq.(13) represent important a-priori information about the image f(x, y) that can be used to reconstruct f(x, y) from a noisy blurred version g(x, y).
NOTE: The choice of sequence {tn}, and the preferred intervals for the least squares fits mentioned above, have been found to work well for a wide class of 8-bit images. However, it should be expected and understood by those skilled in the art that some modification of these parameters may be necessary to accommodate other types of images.
Each row of images in
corresponding (C, α) values in the L2 norm, for the Poisson singular integral Ut. Images of similar objects tend to have approximately equal values for (C, a). Hence, given a blurred Elizabeth Taylor image, in a preferred embodiment of the present invention, one can use (C, α) values corresponding to a sharp Marilyn Monroe image as a-priori information in a deblurring method according to the present invention. More generally, extensive databases of (C, α) values for sharp, noise-free images of similar objects, are a powerful tool in applications such as diagnostic medical imaging, surveillance, environmental monitoring, and several areas of non-destructive evaluation.
The second module is a deblurring module comprising the steps of
1. The first two a-priori numbers in the Singular Integral method are the positive constants E and M in Eq.(5), just as in the Tikhonov-Miller method. The SI method assumes that the unknown exact sharp image fe(x, y) belongs to the Lipschitz space Λ(α, 2, ≈). The first step is choosing singular integral operator SA and positive constants C and α with α<1, such that
∥Stλfe−fe∥2≦C∥f∥2tα/2λ0<t≦t*. (15)
Typical choices are t*=0.1, and λ=0.5, corresponding to the Poisson Singular Integral (PSI) method, or λ=1.0, corresponding to the Gaussian Singular Integral (GSI) method. Since infinitely many distinct choices of λ are possible, the SI method actually encompasses an entire class consisting of infinitely many distinct methods. However, with fixed λ and t*, the required a-priori information in the SI method consists of the four constants {ε, M, C, α}. Examples of pairs (C, α) for diverse images when λ=0.5 and t*=0.1, are given in Tables 1 and 2. The values of C and α are confined to fairly narrow ranges for a wide class of images, and plausible guesses for (C, α) are generally possible.
2. For integer ξ, η, with −J≦ξ, η≦J, let ρ=(ξ2+η2)λ, With t*>0 as in Eq. (15), defining {circumflex over (z)}(ξ, η, λ, t*) by
3. Next, with the given pair (C, α), let δ=α/(2λ), and defining the positive constant Γt* by
4. Using the given pair (ε, M) in Eq. (5), and using the discrete forward Fourier transforms ĝ(ξ, η) and {circumflex over (p)}(ξ, η) of g(x, y) and p(x, y) respectively, forming the Fourier domain function {circumflex over (f)}SI (ξ, η) defined by
In Eq. (18), {circumflex over ({overscore (p)})}(ξ, η) denotes the complex conjugate of {circumflex over (p)}(ξ, η).
5. Using Fast Fourier Transform routines, forming the Inverse Fourier Transform of {circumflex over (f)}SI (ξ, η) in Eq. (18). Reset any negative pixel value to the value zero. With the constant MAX as in Eq. (3), reset any pixel value exceeding MAX to the value MAX. The resulting function, denoted by fSI(x, y) is defined to be the Singular Integral deblurred image.
6. Using commonly available software visualization tools, displaying fSI(x, y) on a computer screen or other imaging device. Examples of successful SI restorations are shown in
7. Creating multiple trial reconstructions. The exact a-priori Lipschitz space information (C, a) may not be available for a given image. However, as shown in TABLES 1 and 2, values of C and a are confined to a fairly narrow range for a diverse class of images. Moreover, similar objects often have approximately equal (C, a) values for fixed A and t*. Useful databases of (C, a) values for various classes of objects can therefore be compiled. In
Starting from an informed, plausible, initial guess for (C, a), it is feasible and helpful to explore computationally the effects of varying the values of C and α. This is achieved by simply varying the constant Γt* in Eq. (18). With fixed λ and t*, each choice of pair (C, α) defines a new constant Γt* in Eq. (17). With the stored pre-computed arrays {circumflex over (p)}(ξ,η), ĝ(ξ, η), and {circumflex over (z)}(ξ,η,λ,t*), multiple Fourier domain functions {circumflex over (f)}SI(λ, η) can be generated in Eq. (18), corresponding to multiple constants Γt*. Parallel computation, if available, may be used to generate and invert the multiple Fourier domain functions {circumflex over (f)}SI(ξ, η). Thus, steps 4 and 5 of the SI method, which constitute the heart of the SI method, can be done in parallel. However, it should be noted that computation of these multiple SI reconstructions is also quite easily accomplished on sequential machines. Visual inspection of the multiplicity of reconstructions, together with some prior knowledge of expected salient features in the exact solution fe(x, y), can guide the user towards selecting the best Singular Integral restoration.
|{circumflex over (n)}(ξ, η)|, |{circumflex over (f)}e(ξ, η)|.
While the preferred embodiments of the present invention have been illustrated and described, it will be understood by those skilled in the art that the image processor as described herein is illustrative, and various changes and modifications may be made to the algorithms and equivalents may be substituted for elements thereof, without departing from the true scope of the present invention. In addition, many modifications may be made to adapt the teachings of the present invention to a particular situation without departing from its central scope. Therefore, it is intended that the present invention not be limited to the particular embodiments disclosed as the best mode contemplated for carrying out the present invention, but that the present invention include all embodiments falling within the scope of the appended claims.
This application claims priority from U.S. Application Ser. No. 60/574,787 and from U.S. Application Ser. No. 60/574,788, both filed May 27, 2004, the contents of both of which are hereby incorporated by reference in their entirety.
Number | Date | Country | |
---|---|---|---|
60574787 | May 2004 | US | |
60574788 | May 2004 | US |