WO2015047538A1

WO2015047538A1 - Psf estimation using multiple smaller images

Info

Publication number: WO2015047538A1
Application number: PCT/US2014/049188
Authority: WO
Inventors: Radka Tezaur
Original assignee: Nikon Corporation
Priority date: 2013-09-30
Filing date: 2014-07-31
Publication date: 2015-04-02

Abstract

A method for estimating an area point spread function (230) of an image area (225) includes the steps of: dividing the image area (225) into a plurality of area blocks (227); and estimating the area point spread function (230) of the image area (225) using a PSF cost function that sums the fidelity terms of at least two of the area blocks (227). The image area (225) can be an entire photograph (14), an image block (222) or a group of photographs (850).

Description

United States Provisional Patent Application

for

PSF ESTIMATION USING MULTIPLE SMALLER IMAGES

of

Radka Tezaur

RELATED APPLICATION

[0001] This application claims priority on U.S. Provisional Application Serial No. 61 /884,649, filed September 30, 2013, and entitled "PSF ESTIMATION USING MULTIPLE SMALLER IMAGES." As far as permitted, the contents of U.S. Provisional Application Serial Number 61 /884,649 are incorporated herein by reference.

BACKGROUND

[0002] Cameras are commonly used to capture an image of a scene that includes one or more objects. Unfortunately, some of the images can be blurred. For example, movement of the camera and/or movement of the objects in the scene during the exposure time of the camera can cause the image to be blurred. Further, if the camera is not properly focused when the image is captured, the image can be blurred. Additionally, there are other types of blur than just motion and defocus blur. For example, blur can be caused by lens aberrations in the image apparatus.

[0003] When blur is sufficiently spatially uniform, a blurred captured image can be modeled as the convolution of a latent sharp image with some point spread function ("PSF") plus noise,

B = K * L + N. Equation (1 ).

[0004] In Equation (1 ) and elsewhere in this document, (i) "B" represents a blurry image, (ii) "L" represents a latent sharp image, (iii) "K" represents a PSF kernel, and (iv) "N" represents noise (including quantization errors, compression artifacts, etc.). [0005] A non-blind PSF problem seeks to recover the PSF kernel K when the latent sharp image L is known. Alternatively, a blind problem seeks to recover both the PSF kernel K and the latent sharp image L. Both the blind and non- blind problems are complex to solve, computationally expensive, and slow. Moreover, many blurry images include areas that further complicate the problem of recovering the PSF kernel. Further, in many blurry images, the PSF varies in the image.

[0006] One common approach to solving a PSF estimation and deconvolution problem includes reformulating it as an optimization problem in which suitable cost functions are minimized. A relatively common type of a cost function used for PSF estimation is a regularized least squares cost function. Typically, a regularized least squares cost function consists of (i) one or more fidelity terms, which make the minimum conform to equation (1 ) modeling of the blurring process, and (ii) one or more regularization terms, which make the solution more stable and help to enforce prior information about the solution, such as sparseness.

[0007] An example of such a PSF cost function is c(K) = L * K - B + γ D * K Equation (2)

[0008] In Equation (2), and elsewhere in this document, (i) c(K) \s the Point Spread Function estimation cost function, (ii) Ddenotes a regularization operator,

(iii) L* K -B is a fidelity term that forces the solution to satisfy the blurring model in Equation (1 ) with a noise term that is small, (iv) ||

is a regularization term that helps to infuse prior information about arrays that can be considered a PSF kernel, (v) γ is a regularization weight that is a selected constant that helps to achieve the proper balance between the fidelity and the regularization terms (vi) the subscript p denotes the norm or pseudo-norm for the fidelity term(s) that results from assuming a certain probabilistic prior distribution modeling, (vii) the superscript p denotes the power for the fidelity term(s), (viii) the subscript q denotes the norm for the regularization term(s) that results from assuming a certain probabilistic prior distribution modeling, and (ix) the superscript q denotes the power for the regularization term(s).

[0009] When the noise derivatives are assumed to follow a Guassian prior distribution, the power p for the fidelity term(s) is equal to two (p=2); and when the image derivatives are assumed to follow a Guassian prior distribution, the power q for the regularization term(s) is equal to two (q=2). This can be referred to as a "Gaussian prior" or a "2-norm".

[0010] Alternatively, when the noise derivatives are assumed to follow a Laplacian prior distribution, the power p for the fidelity term(s) is equal to one (p=1 ); and when the image derivatives are assumed to follow a Laplacian prior distribution, the power q for the regularization term(s) is equal to one (q=1 ). This can be referred to as a "Laplacian prior" or a "1 -norm".

[0011] Still alternatively, when the noise derivatives are assumed to follow a hyper-Laplacian prior distribution, the power p for the fidelity term(s) is less than one (p<1 ); and when the image derivatives are assumed to follow a hyper- Laplacian prior distribution, the power q for the regularization term(s) is less than one (q<1 ). This can be referred to as a "hyper-Laplacian prior" or a "pseudonorm".

[0012] For the example, if a 2-norm is used in the fidelity and regularization terms of Equation (2), then a non-exclusive example of a closed form formula for the cost function of Equation 2 minimum is as follows:

F (L )F (L ) + yF (D )F (D ) Equation (3)

[0013] In Equation (3), the regularization term is necessary to compensate for zeros in F(L).

[0014] Unfortunately, existing PSF cost functions are not entirely adequate. SUMMARY

[0015] The present invention is directed to a system and method for estimating an area point spread function of an image area. In one embodiment, the method includes the steps of: dividing the image area into a plurality of area blocks; and estimating the area point spread function of the image area using a PSF cost function that sums the fidelity terms of at least two of the area blocks. This also allows for possibility to leave out areas blocks in the calculations. Further, this allows for the use of a smaller PSF array (as the PSF needs to be extended only to the size of the area block) and improves numerical stability. As a result thereof, the calculations necessary to determine the area point spread function are simplified, the computational power required to determine the area point spread function is reduced, and the calculation speed is increased.

[0016] As provided herein, the image area can be an entire photograph that is divided into a plurality of image blocks, and each image block can be one of the area blocks that is used in the PSF cost function.

[0017] Alternatively, the image area can be a portion of a photograph that is divided into a plurality of image blocks, and each image block can be one of the area blocks that is used in the PSF cost function.

[0018] Still alternatively, the photograph can be divided into a plurality of image blocks, the image area can be one of the image blocks, and the image block can be divided into the plurality of area blocks that are used in the PSF cost function.

[0019] In yet another embodiment, the image area is a group of photographs that are each blurred with approximately the same area point spread function. In this embodiment, the group of photographs can be divided into the plurality of area blocks used in the PSF cost function. In this example, each of the photographs can be one of the area blocks.

[0020] Additionally, the method can include the steps of estimating a feature value for at least some of the area blocks, and determining which area blocks have the best feature values. In this embodiment, the area blocks with the best feature values can be used for estimating the area point spread function. In this embodiment, the feature value can be evaluated by evaluating a power spectrum, and wherein the area blocks which have the largest power spectrum are assigned the best feature values. Further, in this embodiment, a predetermined number of the area blocks with the best feature values can be used in the estimation of the area point spread function.

[0021] In still another embodiment, the present invention is directed to a system for estimating an area point spread function of an image area, the system including a control system having a processor that (i) divides the image area into a plurality of area blocks; and (ii) estimates the area point spread function of the image area using a PSF cost function that sums the fidelity terms of at least two of the area blocks. For example, the system can be a camera or computer that performs the processes described herein.

BRIEF DESCRIPTION OF THE DRAWINGS

[0022] The novel features of this invention, as well as the invention itself, both as to its structure and its operation, will be best understood from the accompanying drawings, taken in conjunction with the accompanying description, in which similar reference characters refer to similar parts, and in which:

[0023] Figure 1 is a simplified illustration of an image apparatus, a computer having features of the present invention, and a captured image;

[0024] Figure 2A illustrates the captured image of Figure 1 divided into a plurality of image blocks;

[0025] Figure 2B illustrates the image blocks labeled and with a feature value;

[0026] Figure 2C is a simplified illustration of an estimated, area point spread function for the captured image of Figure 2A;

[0027] Figure 3A illustrates the captured image of Figure 1 divided into a plurality of image blocks, and each image block divided into a plurality of mircoblocks;

[0028] Figure 3B illustrates one of the image blocks, with the image block divided into a plurality of mircoblocks, labeled, and with the feature values listed;

[0029] Figure 3C illustrates three microblocks; [0030] Figure 3D illustrates three other microblocks;

[0031] Figure 4 is a simplified illustration of the point spread functions calculated with default microblocks for a number of default image blocks from an image;

[0032] Figure 5 is a simplified illustration of the point spread functions calculated with adaptive microblocks for a number of default image blocks from an image;

[0033] Figure 6 is a simplified illustration of the point spread functions calculated with adaptive microblocks for a number of adaptive image blocks from an image;

[0034] Figure 7 is a simplified illustration of an overall image point spread function for the captured image of Figure 2A;

[0035] Figure 8A is a simplified illustration of a group of photographs having features of the present invention; and

[0036] Figure 8B is a simplified illustration of an estimated, area point spread function for the group of photographs of Figure 8A.

DESCRIPTION

[0037] Figure 1 is a simplified illustration of an image apparatus 10 (e.g. a digital camera), a computer 12 (illustrated as a box), and a blurred, captured image 14 (also referred to as a photograph). The present invention is directed to one or more unique algorithms that are used to estimate a point spread function of an image area (e.g. the image 14) with improved accuracy. If the point spread function is more accurate, the resulting deblurred image will be improved. Moreover, the formulas provided herein can be solved relatively easily. This will speed up the deblurring process and will allow for image deblurring with less powerful processors (e.g. in camera processing).

[0038] As provided herein, the image apparatus 10 can include a capturing system 15 (e.g. a semiconductor device that records light electronically) that captures the image 14, and a control system 16 that uses one or more of the algorithms for estimating the point spread functions for in camera processing. Alternatively, the computer 12 can include a control system 18 that uses one or more of the algorithms for estimating the point spread functions. In either event, the control system 16, 18 can provide a deblurred latent sharp image (not shown) from the blurred image 14. Each control system 16, 18 can include one or more processors and circuits. Further, either of the control systems 16, 18 can include software that utilizes one or more methods and formulas provided herein.

As an overview, as provided herein, a blurred image area can be divided into multiple smaller area blocks, and the new formulas provided herein allow for the simultaneous use of some or all of the smaller area blocks to determine an area point spread function for the image area. As provided herein, the new formulas (e.g. cost functions which sum of the fidelity terms from different image areas) can be treated by a closed form formula that involves the sum of power spectra of the image areas. These new formulas improve the accuracy and speed of estimating the area point spread function for the blurry image area.

[0039] In one embodiment, the image area can be the entire image 14 ("photograph"), and the image 14 can be divided into a plurality of separate area blocks. In this embodiment, the formulas provided herein can be used to accurately estimate the point spread function for the entire photograph using a plurality of the area blocks. This embodiment is described in further detail with reference to Figures 2A-2C below.

[0040] Alternatively, the image area can be just a portion of the image (e.g. an image block), and the image block can be divided into a plurality of separate area blocks. In this embodiment, the formulas provided herein can be used to accurately estimate the point spread function for the image block using a plurality of the area blocks. This embodiment is described in further detail with reference to Figures 3A-3F below.

[0041] Still alternatively, the new formulas provided herein allow for the simultaneous use of multiple photographs from a group of photographs to determine a group point spread function when each of the photographs in the group have approximately the same point spread function. For example, multiple photographs can have approximately the same point spread function when the blur is caused by lens aberrations in the image apparatus 10. In this embodiment, the multiple photographs can be considered the image area, and some or all of the photographs can be considered a separate area block. This embodiment is described in further detail with reference to Figures 8A and 8B below.

[0042] Figure 2A illustrates the captured image 14 of Figure 1 divided into a plurality of image blocks 222. The size, shape and number of the image blocks 222 can be varied. For example, in Figure 2A, each image block 222 is generally rectangular shaped and the image 14 is divided into a seven by seven grid that includes forty-nine, equally sized, image blocks 222. Alternatively, the shape and/or number image blocks 222 can be different than that illustrated in Figure 2A. For example, the image 14 can be divided into a five by five or a nine by nine grid. It should be noted that any of the image blocks 222 can be referred to as a first image block, a second image block, a third image block, etc.

[0043] Additionally, it should be noted that some of the image blocks 222 in Figure 2A have good texture, some of the image blocks 222 have very little texture, and some of the image blocks have many strong edges in certain directions.

[0044] As provided above, new formulas are provided herein that allow for the simultaneous use of multiple smaller area blocks 227 to determine an area point spread function for an image area 225 that includes the multiple area blocks 227. In one embodiment, the image area 225 is the entire photograph 14 or a portion of the entire photograph 14 (e.g. a ¼ of ½ of the photograph 14), and each area block 227 can be a separate one of the image blocks 222. In this embodiment, the new formulas provided can use some or all of the image blocks 222 to accurately estimate the point spread function of the image 14 or a portion of the image 14.

[0045] In one embodiment, the blurring model for an image area 225 that includes a plurality of area blocks 227 can be expressed as follows:

Bi = Li ^* K + Ni B₂ = L₂ ^* K + N₂

Equation (4)

[0046] In Equation (4), for a selected image area 225, (i) the subscript 1 represents one of the area blocks 227 (can also be referred to as the "first area block"); (ii) the subscript 2 represents another of the area blocks 227 (can also be referred to as the "second area block"); (iii) the subscript 3 represents yet another one of the area blocks 227 (can also be referred to as the "third area block"); and (iv) the subscript J represents the "J" area block 227, and J is equal to the number of the image blocks 222 that are being used to determine the area point spread function for the selected image area 225. In certain embodiments, each area block 227 in a given image area 225 is assumed to have been blurred with the same point spread function.

[0047] One non-exclusive example of a PSF cost function based on the blurring model of Equation (4) can be expressed as follows:

/

C

[0048] (K) = Y\\L *K - B. ^P + ^D_*K

7=1

Equation (5)

In this example, the PSF cost function is a generalized regularized least squares cost function for the general value of p. Moreover, in Equation 5 and elsewhere in this document, (i) c(K) \s the Point Spread Function estimation cost function;

(ii) Lj is the latent sharp image at area block j; (iii) Bj is the blurry image at area block j; (iii) "K" represents a PSF kernel; (iv) Ddenotes a regularization operator, (v) γ is a regularization weight that is a selected constant that helps to achieve the proper balance between the fidelity and the regularization terms; (vi) the subscript p denotes the norm or pseudonorm for the fidelity term(s), (vii) the superscript p denotes the power for the fidelity term(s), (viii) the subscript q denotes the norm or pseudonorm for the regularization term(s), (ix) the superscript q denotes the power for the regularization term(s), (x) j is a specific area block, and (xi) J is the total number of area blocks. As provided above, p and q can have a value of 2, 1 or less than 1 .

[0049] For the example, if a 2-norm is used in the fidelity and regularization terms of Equation (5), then a non-exclusive example of a closed form formula for the cost function of Equation 5) minimum is as follows:

F(K)

Equation (6)

where F is the Fourier transform operator.

[0050] In Equation 6, the power spectra of sharp area blocks 227 get added. Because the power spectra are non-negative, if any of them is non-zero, then the sum is non-zero. This makes the algorithm more stable.

[0051] In one embodiment, the present invention is directed to blind deconvolution method where the area point spread function is estimated from the same image area 225 that is going to be deblurred, and the best parts of the image area 225 are being selected to estimate the area point spread. Alternatively, the present invention can be used in situations where the latent sharp image is known, and the PSF is being estimated for test chart images. For this usage, the fact that if any of the power spectra are non-negative, the sum is non-zero can be used when designing the PSF estimation charts.

[0052] In certain embodiments, the test patterns Lj can be selected in such a way that their spectra have no common zeros, or at least as few of them as possible. In blind image deconvolution (both the latent sharp image and PSF are unknown), this can improve the chance that there will be fewer zero or near zero values in the latent sharp image power spectrum.

[0053] Essentially most of the problems with the accuracy of PSF estimation come from the first term in the denominator in the Equation (3) (block power spectrum) being zero or close to zero for some frequencies. The second term in the denominator (regularization) is added to prevent the denominator from being zero. When area blocks 227 are used, the PSF estimation can be more accurate because the first term in the denominator of Equation (6) tends to contain fewer zeros and values close to zero, and thus less regularization is needed.

[0054] In case of motion blur estimation (blind deconvolution), fewer zeros can be achieved by choosing the area blocks 227 in a smart way. In the case of optics blur estimation, when an artificially generated test chart is the ideal sharp image, the l_i to Lj that represent the area blocks of the test chart can be chosen in such a way to stay away from zero as much as possible. With the present algorithm, this term is the sum of power spectra of the area blocks 227, which are all non-negative. This means that if, for example, there are some zeros in the power spectrum of that they can be compensated for by choosing l_2 in such a way that it has zeroes in different locations. Thus, guaranteeing that adding them will not create some new zeroes that will still need to be taken care of.

[0055] It should be noted that the comments regarding zeroes in the power spectrum also apply to cost functions with a power p and/or a power q other than two (2). For these cost functions, a closed form formula equivalent to Equation 6 does not exist, however, when the cost function is minimized by some iterative method (like re-weighted regularized least squares, or alternating minimization derived from variable splitting method), formulas very similar to Equation 6 are used in each iteration.

[0056] The basic PSF formula of Equation (5) can be modified in various ways. For example, additional fidelity terms with derivatives (e.g. directional derivatives D_x and D_y) can be used in the cost function, individual area blocks 227 can get different weights ω , and/or the regularization terms can be different as provided below:

+ D *L * K-D *B

7=1

Equation (7)

[0057] In Equation (7), and elsewhere in this document, (i) C0i is a fidelity weight for the first fidelity term, (ii) (¾ is a fidelity weight for the second fidelity term, (iii) is a regularization weight for the first regularization term, (iv) J2 is a regularization weight for the second regularization term, and (v) q are individual block weights based on an estimated block "quality".

[0058] For the example, if a 2-norm is used in the fidelity term(s) (p=2) and regularization term(s) (q=2) of Equation (7), then a non-exclusive example of a closed form formula for the cost function of Equation (7) minimum is as follows:

Equation (8)

+ co₂G)∑Cj (L.) (L.) + + y₂G)

7=1

In Equation 8, G = F(D_x )F(D_x) + F(D_y )F(D_y).

[0059] It should be noted that other closed form formulas can be used for Equation (7). For example, when post-processed latent sharp image derivatives are used for estimating PSF, the derivative operator (factor ίύ-\ + )2 G) cannot be factored out and the formula is different. [0060] Additionally, it should be noted that other significant modifications are possible to PSF cost function of Equation (7), such as (i) assuming a Laplacian prior (this results in regularization terms with a power ((q=1 )); (ii) assuming a hyper-Laplacian prior (this results in regularization terms with a power of less than one ((q<1 )); and/or (iii) using a spatial mask in either the regularization or fidelity term(s) (or in both). Such cost functions do not have a closed form formula for their minimum, but can be minimized by various iterative numerical methods. One non-exclusive method uses a variable splitting algorithm. One non-exclusive example of a variable splitting algorithm is disclosed in PCT/US13/30227 and PCT/US13/30252. As far as permitted, the contents of PCT/US13/30227 and PCT/US13/30252 are incorporated herein by reference. Another non-exclusive example of a variable splitting technique is described in the paper by Sroubek and Milanfar (F. Sroubek, P. Milanfar: Robust Multichannel Blind Deconvolution via Fast Alternating Minimization, I EEE Trans Im. Proa, 2012), the contents of which are incorporated herein by reference).

[0061] As another example, a PSF cost function having features of the resent invention can be expressed in a more general form as provided below:

+

Equation (9) In Equation (9) and elsewhere (i) COrj , CO1 , CO2 , COn , CO12 , and CO22 are the fidelity weights for the respective fidelity terms; (ii) and are the first order partial derivative operators; (iii) ^xx , ^yy , and ^xy are the second order partial derivative operators; and (iv) the subscript p denotes the norm for the fidelity term(s), and (v) the superscript p denotes the power for the fidelity term(s).

[0062] The exact regularization term(s) utilized in Equation (9) can vary. In one non-exclusive embodiment, the regularization term(s) can be expressed as follows: γΛΚ + γ₂ ω_χ * κ + D_y * K Equation (10).

Thus,

Re gularizatbnTernis) = γ_ί K ^q + γ₂( Ό_χ * Κ + D_y * K ) In

Equation (10) and elsewhere, the gamma 1 and gamma 2

are regularization term weights. It should be noted that the fidelity term weights (CO- omega) and the regularization term weights (γ-gamma) can be chosen to be zero

("0") or another weight. As a result thereof, in Equation (9), this PSF cost function covers cases both with and without derivatives. In one embodiment, the PSF cost function can be assumed to have a Gaussian prior (q=2) to allow for a closed form formula for the minimum. Alternatively, the PSF cost function can (i) assume a Laplacian prior (this results in regularization terms with a norm of one ((q=1 )); or (ii) assuming a hyper-Laplacian prior (this results in regularization terms with a norm of less than one ((q<1 )). [0063] For the example, if a 2-norm is used in the fidelity term(s) (p=2) and regularization term(s) (q=2) of Equation (9), then a non-exclusive example of a closed form formula for the cost function of Equation (9) minimum is as follows:

1 /

Terml^ FiL FiB_j)

F(K) =

Term F(L_j)F(L_j ) + Term 2

Equation (1 1 )

In Equation (1 1 ), the value of Terml will depend on the fidelity term(s) used, and value of Term2 will depend upon the regularization term(s) used in Equation (9). It should be noted that the closed form formula of Equation (1 1 ) includes (i)

the denominator where F is a Fourier transform. The numerator term

can be described as (i) the sum over j of the products of the complex conjugate of the Fourier transform of the latent sharp image number j and the Fourier transform of the blurry image number j, or (ii) the sum over j of the Fourier transforms of the cross-correlations between the latent sharp image number j and blurry image number j. Further, the denominator term )

can be described as the sum of the power spectra of the latent sharp images (blocks).

[0064] Somewhat similarly, if other than a 2-norm is used in the fidelity term(s) (ρΦ2) and regularization term(s) (q+2) of Equation (9), then an iterative algorithm can be used for minimizing the cost function (like the one resulting from variable splitting). Equation (1 1 .1 ) can be used when the power p is not equal to two (2), and the cost function is minimized by variable splitting, as such formula needs to be evaluated in each iteration. The algorithm may include as one of the steps the evaluation of a closed form formula as provided in Equation (1 1 .1 ) as follows:

Equation (1 1 .1 )

In Equation (1 1 .1 ), the value of Terml will again depend on the fidelity term(s) used, and value of Term2 will again depend upon the regularization term(s) used in Equation (9). It should be noted that the closed form formula of Equation in the numerator, and

where F is a Fourier Transform.

Further, Equation (1 1 .1 ) is similar to Equation (1 1 ), except with Aj is used instead of Bj (the blurry image at area block j) one or more times. As used herein, Aj is a modified image at block j that can be computed using (i) the blurry image at block j (Bj), or (ii) a combination of the blurry image at block j (Bj) and the latent sharp image at area block j (Lj). As a non-exclusive example, if p=1 (Laplacian prior fidelity term), then Aj can be either the result of soft thresholding applied to the blurry image Bj (altered version of Bj), or it can be the result of soft thresholding applied to the residuum Bj - K^*l_j, depending on the choice of the auxiliary variable in the variable splitting algorithm.

[0065] For Equation (1 1 .1 ), the numerator term

can be described as (i) the sum over j of the products of the complex conjugate of the Fourier transform of the latent sharp image number j and the Fourier transform of the modified image number j, or (ii) the sum over j of the Fourier transforms of the cross-correlations between the latent sharp image number j and modified image number j. Further, the denominator term )

' can be described as the sum of the power spectra of the latent sharp images (blocks).

[0066] In one embodiment, if the regularization terms of Equation (10) are used in Equation (9), Terml and Term2 can be expressed as follows:

Terml = (O₀ + O_iF(D_x)F{D_x) + C0₂F{D )F(D ) +

m F(D F(DJ + m_l2F(D_v )F(D_xy ) + ω₂₂ F{D_yy )F(D_yy )

Equation (12)

Term! = γ_{ι +} γ₂ (F(D_x )F(D_x ) + F(D_y )F{D_y ))

Equation (13)

[0067] It should be noted that other equivalent equations can be used for Terml and/or Term2 depending upon the regularization term utilized. Further, the equation for Term2 can be arbitrary.

[0068] It also should be noted that other closed form formulas can be used for Equation (9). For example, when post-processed latent sharp image derivatives are used for estimating PSF, the derivative operator Terml cannot be factored out and the formula is different.

[0069] It should also be noted that other types of regularization term(s) can be utilized in the PSF cost functions provided herein. As non-exclusive examples, the regularization term(s) can include one or more of the terms (i) ||K| |q^q with any q (q=2, q=1 , q<1 ); (ii) ||D_a ^*K| |q^C| with any q (q=2, q=1 , q<1 ) and any direction of derivative a; (iii) ||D^*K| |q^q with any q (q=2, q=1 , q<1 ) and any regularization operator D (e.g. Laplacian, but it could represent also a derivative of order higher than 1 ) ; and/or (iv) || |grad(K) | | | _q with any q (q=2, q=1 , q<1 ) where |grad(K)| denotes the magnitude of the gradient of K.

[0070] It should be noted that Equation (9) can be modified to include individual wei hts for the area blocks as follows:

Re gulariz tbnTerm(s).

Equation (14)

[0071] As yet another non-exclusive example, a PSF cost function having features of the present invention can be expressed in a more general form as rovided below:

Equation (15)

[0072] In Equation (15), and elsewhere in this document, (i) K is the PSF kernel to be estimated, (ii) Κ_ρΓβν is the postprocessed PSF kernel from a previous iteration, (iii) Lj, Lj_iX , Lj _y are the jth latent sharp image and its derivatives with respect to x and y, (iv) Bj, Bj_iX , Bj _y are the jth blurry image and its derivatives with respect to x and y, (v) D_x , D_y are convolution kernels for derivatives with respect to x and y, (vi) CO is a general fidelity weight for the first fidelity term, (vii) Cj are individual block weights based on an estimated block

"quality", and (viii) J-_\ is a regularization weight for the first regularization term, (ix) γ₂ is a regularization weight for the second regularization term, (x) the subscript p denotes the norm for the fidelity term(s), (xi) the superscript p denotes the power for the fidelity term(s), (xii) the subscript q denotes the norm for the regularization term(s), and (xiii) the superscript q denotes the power for the regularization term(s). As provided above, p and q can have a value of 2, 1 or less than 1 . It should be noted that Lj, Lj_iX , can possibly be postprocessed.

[0073] For the example, if a 2-norm is used in the fidelity term(s) (p=2) and regularization term(s) (q=2) of Equation (15), then a non-exclusive example of a closed form formula for the cost function of Equation (15) minimum is as follows:

Equation (16)

\ 2

In Equation 16, Z = ( F { D _Z ) + F ( D _y )

[0074] It should be noted that other closed form formulas can be used for Equation (15).

[0075] As provided herein, the algorithms used for solving the blind deconvolution problem estimate both the point spread function kernel K and the latent sharp image L. In certain embodiments, these algorithms are iterative and usually they alternate, improving the current estimates of the point spread function and the latent sharp image. The steps of improving the current estimate of the point spread function kernel and the current estimate of the latent sharp image are typically in the form of solving an optimization problem. For blind deconvolution, (i) the point spread function estimation cost function is formulated and minimized to find a solution for the point spread function blur kernel estimate, and (ii) the latent sharp image estimation cost function is formulated and minimized to find a solution for the latent sharp image estimate. With this design, the point spread function estimation and the latent sharp image estimation can be solved in iterations by minimizing the point spread function estimation cost function, and by minimizing the latent sharp image estimation cost function in alternating fashion.

[0076] In addition to alternating minimization of the point spread function estimation cost function, and the latent sharp image estimation cost function, it is also common to use a multi-resolution approach, where first, the coarse versions of the latent sharp image L and the point spread function kernel K are estimated using a heavily downsampled blurry image, and then these estimates are refined at higher and higher resolution, until the full resolution of the given blurred image is reached. The number of iterations for each resolution can be varied to achieve desired level of deblurring.

[0077] The latent sharp image cost function used in the optimization problem can vary. One non-exclusive example of a suitable latent sharp image cost function c(L) can be expressed as follows: -^K*^LJ. + )

Equation (17)

In Equation (17), and elsewhere in this document, (i) K is the PSF kernel, (ii) L is the latent sharp image, (iii) Lj, Lj_iX , L_j>y are the jth latent sharp image and its derivatives with respect to x and y, (iv) B_j, B_jiX , B_{j y} are the jth blurry image and its derivatives with respect to x and y, (v) D_x , D_y are convolution kernels for derivatives with respect to x and y, (vi) COrj is the fidelity weight for the first fidelity term, (vii) COi is the fidelity weight for the second fidelity term, (viii) y is a regularization weight for the regularization term, (ix) the subscript p denotes the norm for the fidelity term(s), (x) the superscript p denotes the power for the fidelity term(s), (xi) the subscript q denotes the norm for the regularization term(s), and (xii) the superscript q denotes the power for the regularization term(s). As provided above, p and q can have a value of 2, 1 or less than 1 .

[0078] Importantly, other latent sharp image cost functions can be used than provided in Equation (17).

[0079] Referring back to Figure 2A, the present methods provided herein can be applied to an entire, several megapixel image 14 (seen as one huge image area 225), with the algorithm using only a few small area blocks 227 (i.e., even below 1 % of the total area of the image in certain embodiments) to calculate the area point spread function to improve speed. In this example, the best parts (e.g. the best area blocks 227) of the image area 225 can be selected to estimate the area point spread function with improved accuracy. On the other extreme, the algorithm could use many, relatively small area blocks 227, with close to one hundred percent (100%) of the image area 225 being used to calculate the area point spread function to improve accuracy and stability.

[0080] As provided herein, using multiple area blocks 227 to calculate the area point spread function for each image area 225 improves the accuracy of the area point spread function. One reason is the fewer zeros and near zero values in the combined area block power spectrum will allow for the use of less regularization. Another reason is that the PSF array size needs to be extended only to the size of the area block 227, not the entire image area 225 (which essentially means that additional information is being inserted about the size of the support of area PSF). [0081] Figure 2B illustrates the image area 225 of Figure 2A in a grid format with the forty-nine area blocks 227, and without the actual image. In Figure 2B, the area blocks 227 have been labeled A-WA for ease of discussion.

[0082] In one embodiment, all of the area blocks 227 of the image area 225 are used in the algorithm to estimate the area PSF cost function to improve accuracy and stability. Alternatively, to simplify calculations, in certain embodiments, for a selected image area 225, not all of the area blocks 227 (image blocks 222 in this example) are used to calculate the area PSF for the image area 225. Not using all area blocks 227 reduces the amount of data (size of the image area 225) that is processed and makes computations faster.

[0083] As alternative, non-exclusive examples, only one, two, three, four, five, six, seven, ten, fifteen, twenty, twenty-five, thirty, or thirty-five area blocks 227 can be used to estimate the area PSF for a selected image area 225. Stated in a different fashion, in alternative, non-exclusive examples, only approximately two, five, ten, fifteen, twenty, thirty, forty, fifty, sixty or seventy percent of the area blocks 227 in a selected image area 225 are used in the calculation of the area PSF for that image area 225. It should be noted that any percentage between one and one hundred (1 -100) of the area blocks 227 in the image area 225 can be used in the area PSF calculation.

[0084] In one embodiment, a default number and selection (fixed pattern) of area blocks 227 ("default area block method") are used to estimate the area PSF for the image area 225. As a non-exclusive example, for each image area 225, five of the area blocks 227 can be used to estimate the area PSF of that image area 225. As an example, the area blocks 227 labeled "Q", "S", "Y", ΈΑ", and "GA" can be the default area blocks 227 that are used to estimate the area PSF for that image area 225. Alternatively, a different number of area blocks 227 (less than five or greater than five) and/or other area blocks 227 (e.g. "A", "B", "C", etc) can be used as the default area blocks 227 for estimating the area PSF of that image area 225.

[0085] In another embodiment, the accuracy of the estimation of area point spread function can be improved with adaptive selection of the area blocks 227 ("adaptive area block method") in that image area 225. In this embodiment, the goal is to use only good area blocks 227 in the calculation of the area PSF kernel. With this design, the estimated area PSF will be more accurate, because bad area blocks 227 are not used in the calculation of the area PSF.

[0086] In one embodiment, one or more of the area blocks 227 (labeled A-WA in this example) can be evaluated and assigned a block quality based on a quality of that area block 227. For example, the control system 18, 20 (illustrated in Figure 1 ) can first process each of the area blocks 227 (A-WA in this example) to identify the area blocks 227 which would be best for determining the area point spread function for the image area 225.

[0087] More specifically, each of the area blocks 227 can be evaluated and assigned a feature value 229 based on the quality of features in that area block 227. In this example, the control system 18, 20 can assign the feature value 229 to each of the area blocks 227, with the feature value 229 being based on the quality of its power spectrum (e.g. strength, and isotropy). In this embodiment, (i) the feature value 229 of each area block 227 is larger (better) for area blocks 227 that have a larger power spectrum, and the feature value 229 of each area block 227 is smaller (worse) for area blocks 227 that have a smaller power spectrum. Thus, (i) the area blocks 227 with a large power spectrum have significant amount of texture (edges in multiple directions) and are assigned larger (better) feature values 229, and can be labeled as "good", and (ii) area blocks 227 with a low power spectrum contain little texture (mainly smooth) and are assigned smaller (worse) feature values 229, and can be labeled as "bad". Further, in certain embodiments, area blocks 227 with texture in only one direction (strong edges in only one direction) will be assigned smaller (worse) feature values 229, and can be labeled as "bad".

[0088] The scale used for the feature values 229 can vary. In one example, the scale is 0-200, with low feature values having a score that is closer to zero, and high features values having a score that is closer to two hundred. Alternatively, the values can be reversed with the low score being assigned to area blocks 227 with good features, and the high score being assigned to area blocks 227 with bad features. In these examples, the feature value 229 has a numerical value. Alternatively, the feature value 229 can have another type of value or grade (e.g. good or bad).

[0089] It should be noted that a separate feature value 229 is provided for each of the area blocks 227 in Figure 2B. In this example, the feature values 229 are made up to facilitate this discussion and are not actual feature values 229. In the example illustrated in Figure 2A, area blocks "LA", ΈΑ", "W", Έ", and "H" have the highest feature values, and area blocks "G", "N", "U", "F", and "JA" have the lowest feature values.

[0090] In this embodiment, if it is desired to use the five best area blocks 227, then area blocks 227 labeled "LA", "EA", "W", "E", and "H" are used in the PSF cost function (e.g. Equations 5, 7, or 9) to estimate the area PSF kernel.

[0091] In certain embodiments, the area blocks 227 with better feature values 229 can be given more weight than the area blocks 227 with worse feature values 229 in the PSF cost function. As provided above, each area block 227 can be assigned some numerical feature value 229 (or block quality measure) based on the quality of the area block 227. Subsequently, these feature values

229 can be transformed into block weights Cj based on the numerical feature value 229. In one embodiment, (i) the block weight Cj is bigger for area blocks 227 that have good feature values 229 (e.g., have stronger and more isotropic power spectrum) and, (ii) the block weight Cj is smaller for area blocks 227 that have bad feature values 229 (e.g., low texture). Additionally, the weights Cj can be normalized in some way so as not to mess up with the regularization weights. For example the block weights can be made to sum up to one (1 ) or to s (the number of area blocks 227 used).

[0092] After the block weights q are determined, they can be used in the PSF cost function (see e.g. Equations 7, 8, and 14-16) to give more weight to the area blocks 227 with better feature values 229 and less weight to the area blocks 227 with worse feature values 229. [0093] Figure 2C is a simplified illustration of the area point spread function 230 of the image area 225 (photograph 14) calculated with multiple of the best area blocks 227 (image blocks 222). In this embodiment, because the image area 225 is the entire image 14, the area point spread function 230 can also be referred to as the image point spread function.

[0094] Figure 3A illustrates the captured image 14 of Figure 1 divided into the plurality of image blocks 222 (separated with relatively thick lines) and each image block 222 is divided into a plurality of separate microblocks 324 (separated with thinner, dashed lines). The size, and number of the microblocks 324 in each image block 222 can be varied. For example, in Figure 3A, each microblock 324 is generally rectangular shaped and each image block 222 is divided into a five by five grid that includes twenty-five, equally sized, microblocks 324. Alternatively, the shape and/or number microblocks 324 can be different than that illustrated in Figure 3A. For example, one or more of the image blocks 222 can be divided into a three by three or a seven by seven grid.

[0095] For any given image block 222, any of the microblocks 324 can be referred to as a first microblock, a second microblock, a third microblock, etc. It should be noted that for each of the image blocks 222, some of the microblocks 324 can have good texture, some of the microblocks 324 can have very little texture, and some of the microblocks 324 can have many strong edges.

[0096] As provided above, new formulas are provided herein that allow for the simultaneous use of multiple smaller area blocks 327 to determine an area point spread function for an image area 325 that includes the multiple area blocks 327. In this embodiment, each image area 325 can be one of the image blocks 222, and each area block 327 can be a separate one of its microblocks 324. Thus, in the example illustrated in Figure 3A, the image 14 includes forty-nine separate image areas 325, and each image area 325 includes twenty-five separate area blocks 327.

[0097] Moreover, in this embodiment, for each image area 325 (image block 222 in this example), the new formulas (e.g. Equations 5, 7, and 9) can use some or all of its area blocks 327 (microblocks 324 in this example) to accurately estimate the area point spread function (also referred to as the "block PSF" in this example) of the respective image area 325. Subsequently, the area point spread function of two or more of the image areas 325 can be combined or blended to estimate a single image point spread function of the image 14, or estimate a spatially varying image point spread function for the blurry image 14.

[0098] Figure 3B illustrates one of the image areas 325 (image blocks 222 in this example), with the image areas 325 divided into twenty-five area blocks 327 (mircoblocks 324 in this example), and without the actual image. In Figure 3B, the area blocks 327 have been labeled A-Y for ease of discussion. It should be noted that a separate feature value 329 is also provided for each of the area blocks 327 in Figure 3B. In this example, the feature values 329 are made up to facilitate this discussion and are not actual feature values 329.

[0099] Similar to the example provided above, in one embodiment, for each image area 325, all of the area blocks 327 in that image area 325 can be used in the algorithm to estimate the area PSF to improve accuracy and stability. In this example, area blocks labeled A-Y can be used in the area PSF estimation.

[00100] Alternatively, to simplify calculations, for one or more of the image areas 325, a plurality but not all of its area blocks 327 are used to calculate the area PSF. Stated in another fashion, several, but not all of the area blocks 327 of a selected image area 325 are used for estimating the area PSF for that image area 325. Not using all area blocks 327 reduces the amount of data that is processed and makes computations faster. For each image area 325, the number of area blocks 327 used can vary similar as provided above in reference to the discussion of Figure 2B.

[00101] In one embodiment, a default number and selection (fixed pattern) of area blocks 327 ("default area block method") are used to estimate the area PSF for each of the selected image areas 325. As a non-exclusive example, for each image area 325, five of the area blocks 327 can be used to estimate the area PSF of that image area 325. In one example, the area blocks 327 labeled "G", "I", "M", "Q", and "S" can be the default area blocks 327 that are used to estimate the area PSF for that image area 325. Alternatively, a different number of area blocks 327 (less than five or greater than five) and/or other area blocks 327 (e.g. "A", "C", "D", etc) can be used as the default area blocks 327 for estimating the area PSF of that image area 325.

[00102] In another embodiment, the accuracy of the estimation of area PSF can be improved with adaptive selection of area blocks 327 ("adaptive area block method") in that image area 325. For example, for each image area 325, the control system 18, 20 (illustrated in Figure 1 ) can first process the area blocks 327 to identify the area blocks 327 which would be best for determining the area PSF for that respective image area 325.

[00103] Similar to the embodiment provided above, the control system 1 8, 20 can assign the feature value 329 to each of the area blocks 327, with the feature value 329 being based on the quality of its power spectrum (e.g. strength, and isotropy). In this embodiment, the feature value 329 of each area blocks 327 is larger (better) for area blocks 327 that have a stronger power spectrum with more high frequency content, and the feature value 329 of each area blocks 327 is smaller (worse) for area blocks 327 that have a weaker power spectrum.

[00104] In Figure 3B, the feature value 329 for each area block 327 is listed under the letter. For example, area block 327 "A" has a feature value of 3.2, while area block 327 "U" has a feature value of 122.3. In this example, area blocks 327 "U", "K", and "P" have the highest feature values 329 while area block 327 "A", "B", and "C" have the lowest feature values 329.

[00105] Figure 3C illustrates three area blocks 327, labeled "A", "B", and "C". In this example, the area blocks 327 are pretty smooth and do not have very many edges and texture. Thus, these area blocks 327 will have a low feature value 329 and can be classified as "bad area blocks".

[00106] Figure 3D illustrates three additional area blocks 327, labeled "U", "P", and "K". In this example, the area blocks 327 have good texture, with strong edges in multiple directions. Thus, these area blocks 327 will have a high feature value and can be classified as "good microblocks".

[00107] Referring back to Figure 3B, as provided herein, in certain embodiments, for each image area 325, only the area blocks 327 with the best feature values are used to calculate the area PSF. Selecting the area blocks 327 with the best feature values helps to improve the PSF estimation accuracy because these area blocks 327 have high texture in many directions. In one embodiment, for each image area 325, only a predetermined number of the area blocks 327 with the highest feature values are used to calculate the PSF. For example, if the predetermined number is three, then the three area blocks 327 with the highest feature values will be used. In Figure 3B, area blocks 327 "U", "K", and "P" have the three highest feature values 329 and these area blocks 327 will be used in the area PSF calculations. Alternatively, the predetermined number can be more than three or less than three.

[00108] In another embodiment, for each image area 325, only area blocks 327 having a feature value equal to or greater than a predetermined value threshold will be used to calculate the area PSF. For example, if the predetermined value threshold is forty, then the area blocks 327 having a feature value 329 equal to or greater than forty will be used. In Figure 3B, area blocks 327 "U", "K", "P", and "V" have a feature value 329 that is greater than forty, and these area blocks 327 will be used in the area PSF calculations. Alternatively, the predetermined threshold can be more than forty or less than forty.

[00109] In certain embodiments, the method provided herein utilizes the area PSF estimated for each of the individual image areas 325 (i) to estimate a single image point spread function for the entire blurry image 14, or (ii) to estimate a spatially varying image point spread function for the blurry image 14.

[00110] Alternatively, the method provided herein can estimate the area PSF for less than all of the image areas 325 illustrated in Figure 3A to reduce the number of computations. For example, in one embodiment, the method ("default block method") provided herein calculates the area PSF for only a default, predetermined number (less than the total number) and/or default pattern (fixed pattern) of image areas 325, and utilizes the area PSF estimated for the default, predetermined number and default pattern of individual image areas 325 (i) to estimate the single image point spread function for the entire blurry image 14, or (ii) to estimate the spatially varying image point spread function for the blurry image 14.

[00111] The default, predetermined number and default pattern of individual image areas 325 used to calculate the image PSF can be varied. As a nonexclusive example, Figure 4 is a simplified illustration of the forty-nine image blocks 222, with the individual area point spread functions 430 calculated for twenty-four image areas 325 (image blocks 222). In this example, the default pattern is every other image area 325 in every other row is calculated. Alternatively, the default, predetermined number can be greater or lesser than twenty-four, and the default pattern can be different than that illustrated in Figure 4.

[00112] In the example illustrated in Figure 4, the method provided herein calculates the area PSF 430 for only the twenty-four (of the forty-nine) image areas 325, and utilizes the area PSF 430 for only the twenty-four image blocks 222 to calculate the image point spread function. In this example, because the block PSF 430 is not calculated for all of the image areas 325, the computational requirements are less, and the speed of the calculations is increased.

[00113] It should be noted in the example illustrated in Figure 4, that the default area blocks 327 (microblocks illustrated in Figure 3A) are used to calculate the area PSF 430 for the default, predetermined number of individual image areas 325. This can result in inaccuracies due to bad area blocks 327 being used in the calculations.

[00114] In another embodiment, the adaptive area block method is used to calculate the area PSF for the default, predetermined number of individual image areas 325. This can result in better accuracies because the accuracy of the area PSFs is increased due to good area blocks 327 being used in the calculations. Figure 5 is a simplified illustration of the forty-nine image blocks 222, with the individual area point spread functions 530 calculated for twenty-four image areas 325 (image blocks 222) using the adaptive area block method described above. In this example, the default pattern is similar to that illustrated in Figure 4. It should be noted that comparing Figures 4 and 5, the area PSF's 530 illustrated in Figure 5 are probably more accurate because of the use of the adaptive area block method.

[00115] In yet another embodiment, an adaptive method ("adaptive image area selection method") is also used to select the image areas 325 in which the area PSFs will be calculated.

[00116] Figure 6 is a simplified illustration of the forty-nine image blocks 222, with the individual area point spread functions 630 calculated for twenty-four image areas 325 (image blocks 222) using the adaptive selection of the image areas 325, and the adaptive area block method described above to calculate the area PSFs 630 for these image areas 325. For example, the area PSFs 630 can be estimated for only the image areas 325 with the highest combined feature value (good blocks only). In this example, for each image area 325, the feature values 329 (illustrated in Figure 3B) for its area blocks 327 (illustrated in Figure 3B) are combined to calculate a combined feature value for each image area 325. For example, the image area 325 of Figure 3B has a combined feature value of 776.5.

[00117] Subsequently, the area PSF 630 can be calculated for only the image areas 325 with the highest combined feature values. Selecting the image areas 325 with the best combined feature values helps to improve the area PSF estimation accuracy because these image areas 325 have high texture in many directions. Further, the overall image PSF will be improved because of the improved accuracy of the area PSFs 630.

[00118] In one embodiment, the area PSF 630 for only a predetermined number of the image areas 325 with the highest combined feature values will be calculated, and only the area PSFs 630 from these image areas 325 are used to calculate the overall image PSF. For example, if the predetermined number is twenty-four, then the area PSF 630 for the twenty-four image areas 325 with the highest combined feature values will be calculated and used. Alternatively, the predetermined number can be more than three or less than twenty-four.

[00119] In another embodiment, the area PSF 630 will only be calculated for the image blocks 222 having a combined feature value equal to or greater than a predetermined combined threshold, and only the area PSFs 630 from these image areas 325 are used to calculate the overall image PSF.

[00120] The method used to calculate the overall image PSF from the multiple area PSFs 430, 530, 630 can vary. In certain embodiments, the overall image PSF is generated by synthesizing one or more (e.g. all) of the area PSFs 430, 530, 630. The overall image PSF can represent the blur for the entire image or the pattern of image PSFs 430, 530, 630 can be studied to determine a spatially varying overall image PSF. By dividing the image 14 into a plurality of image areas 325, the present method is better able to handle slowly spatially varying, image PSF of the image 14.

[00121] A method for synthesizing multiple PSFs is disclosed in PCT/US1 3/29655 filed on March 7, 2013 and entitled "Globally Dominant Point Spread Function Estimation". As far as permitted the content of PCT/US1 3/29655 are incorporated herein by reference.

[00122] Figure 7 is a simplified illustration of an estimated, overall image point spread function 732 for the captured image 14 of Figure 1 , generated by synthesizing area PSFs 630 of Figure 6.

[00123] Figure 8A is a simplified illustration of a group 850 of photographs 814A, 814B, 814C (each illustrated as a box), and Figure 8B is a simplified illustration of an estimated, area point spread function 830 for the group of photographs of Figure 8A. In this embodiment, each of the photographs 814A, 814B, 814C has been blurred with approximately the same PSF. This situation occurs when the blur is caused by lens aberrations in the image apparatus 1 0. The number of photographs in the group 850 can vary. For example, the number can be less than three or more than three.

[00124] As provided herein, the new formulas (Equations 5, 7, 9) provided above allow for the simultaneous use of multiple photographs 814A, 814B, 814C from the group 850 of photographs to determine the area point spread function 830 when each of the photographs 814A, 814B, 814C in the group 850 have approximately the same point spread function. [00125] In this embodiment, the group 850 of photographs can be considered the image area 825, and some or all of the photographs 814A, 814B, 814C can be considered a separate area block 827. Stated in another fashion, in this embodiment, the image area 825 (the group 850 of photographs) can divided into multiple area blocks 827 (individual photographs 814A, 814B, 814C or portions thereof). Subsequently, the area (group) PSF 830 of the image area 825 can be calculated using one of the PSF cost functions (e.g. Equations 5, 7, and 9) that sums the fidelity terms of at least two of the area blocks 825.

[00126] In one embodiment, all of the area blocks 827 of the image area 825 are used in the algorithm to estimate the area PSF cost function to improve accuracy and stability. Alternatively, to simplify calculations, not all of the area blocks 827 are used to calculate the area PSF for the image area 825. Similar to the embodiments provided above, the selection of the area blocks 827 can be adaptive.

[00127] While the current invention is disclosed in detail herein, it is to be understood that it is merely illustrative of the presently preferred embodiments of the invention and that no limitations are intended to the details of construction or design herein shown other than as described in the appended claims.

Claims

What is claimed is:

1 . A method for estimating an area point spread function of an image area, the method comprising the steps of:

dividing the image area into a plurality of area blocks; and estimating the area point spread function of the image area using a

PSF cost function that sums the fidelity terms of at least two of the area blocks.

2. The method of claim 1 wherein the image area is an entire photograph, and wherein the step of dividing includes dividing the photograph into the plurality of area blocks.

3. The method of claim 1 wherein the image area is a portion of a photograph, and wherein the step of dividing includes dividing the portion of the photograph into the plurality of area blocks.

4. The method of claim 1 wherein the image area is a group of photographs that are each blurred with approximately the same area point spread function, and wherein the step of dividing includes dividing the group of photographs into the plurality of area blocks.

5. The method of claim 1 wherein the step of estimating an area point spread function includes minimizing the PSF cost function having the form as follows:

<¾2∑||# L, *K -D„* B, + Re gulariz tbnTern s). where (i) c(K) is the PSF cost function, (ii) L^. is the latent sharp image at area block j, (iii) K is the PSF kernel, (iv) B- j_s the blurry image at area block j, (v) D_x and D_y are first order partial derivative operators, (vi) D_{xx j} D_{yy ;} and D_xy are second order partial derivative operators, (vii) <% , , co2 , ro , ω₁₂ , and (B22 are fidelity weight parameters for the fidelity terms, (viii) s is a total number of area blocks used estimate the area point spread function of the image area, and (ix) p is a power and a norm for the fidelity term(s).

6. The method of claim 5 further comprising the step of minimizing the PSF cost function using a closed form formula that includes (i) (Z,) (Z,)

in the denominator; where F is a Fourier Transform.

7. The method of claim 5 wherein at least one of the fidelity weight parameters has a value of zero.

8. The method of claim 5 wherein the step of estimating an area point spread function includes at least one of the RegularizationTerm(s) being selected from the group that includes (i) ||K||_q ^q, (ii) ||D_a ^*K||_q ^q, (iii) ||D^*K||, and/or (iv)

II |grad(K)| | | _q ^q; where (i) q is a power and a norm for the RegularizationsTerm(s) and has a value of two, one, or less than one, (ii) D is any regularization operator, (iii) a is any direction of derivative, and (iv) |grad(K)| denotes a magnitude of a gradient of K.

9. The method of claim 5 further comprising the steps of estimating a feature value for at least some of the area blocks, and determining which area blocks have the best feature values, and wherein the area blocks with the best feature values are used in the step of estimating the area point spread function.

10. The method of claim 5 further comprising the step of minimizing the PSF cost function using a closed form formula that includes (i) in

the denominator; where F is a Fourier Transform, and is a modified image at area block j.

1 1 . The method of claim 1 further comprising the steps of estimating a feature value for at least some of the area blocks, and wherein the step of estimating the area point spread function includes giving more weight in the PSF cost function to area blocks with better feature values and giving less weight to the area blocks with worse feature values.

12. The method of claim 1 further comprising the steps of estimating a feature value for at least some of the area blocks, and determining which area blocks have the best feature values, and wherein the area blocks with the best feature values are used in the step of estimating the area point spread function.

13. The method of claim 12 wherein the step of estimating a feature value includes the step of evaluating a power spectrum for at least some of the area blocks, and wherein the step of determining which area blocks have the best feature values includes estimating which of the area blocks have the largest power spectrum.

14. The method of claim 1 further comprising the steps of estimating a feature value for at least some of the area blocks, and determining which area blocks have the best feature values, and wherein a predetermined number of the area blocks with the best feature values are used in the step of estimating the area point spread function.

15. A system for estimating an area point spread function of an image area, the system comprising:

a control system that includes a processor that (i) divides the image area into a plurality of area blocks; and (ii) estimates the area point spread function of the image area using a PSF cost function that sums the fidelity terms of at least two of the area blocks.

16. The system of claim 15 wherein the image area is an entire photograph, and wherein the control system divides the photograph into the plurality of area blocks.

17. The system of claim 15 wherein the image area is a portion of a photograph, and wherein the control system divides the portion of the photograph into the plurality of area blocks.

18. The system of claim 15 wherein the image area is a group of photographs that are each blurred with approximately the same area point spread function, and wherein the control system divides the group of photographs into the plurality of area blocks.

19. The system of claim 15 wherein the PSF cost function is expressed as follows:

where (i) c(K) is the PSF cost function, (ii) L_j is the latent sharp image at area block j, (iii) K is the PSF kernel, (iv) B^. j_s the blurry image at area block j, (v) D_x and D_y are first order partial derivative operators, (vi) D_xx , D_{yy j} and are second order partial derivative operators, (vii) <% , coi , <¾ , ro , coi ₂ , and ω₂2 are fidelity weight parameters for the fidelity terms, (viii) s is a total number of area blocks used estimate the area point spread function of the image area, and (ix) p is the prior.

20. The system of claim 19 wherein the control system minimizes the PSF cost function using a closed form formula that includes (i)

in the numerator, and (ii)

in the denominator; where F is a Fourier

Transform.

21 . The system of claim 19 wherein at least one of the fidelity weight parameters has a value of zero.

22. The system of claim 19 wherein at least one of the RegularizationTerm(s) is selected from the group that includes (i) ||K| |q^q, (ii)

| | Da*K| |q^q, (iii) ||D^*K||, and/or (iv) | | |grad(K)| | | _q ^q; where (i) q is a power and a norm for the RegularizationsTerm(s) and has a value of two, one, or less than one, (ii) D is any regularization operator, (iii) a is any direction of derivative, and (iv) |grad(K)| denotes a magnitude of a gradient of K.

23. The system of claim 19 wherein the control system estimates a feature value for at least some of the area blocks, and determines which area blocks have the best feature values, and wherein the area blocks with the best feature values are used by the control system for estimating the area point spread function.

24. The system of claim 15 wherein the control system minimizes the PSF cost function using a closed form formula that includes

here F is a Fourier

Transform, and is a modified image at area block j.

25. The system of claim 15 wherein the control system estimates a feature value for at least some of the area blocks, and determines which area blocks have the best feature values, and wherein the area blocks with the best feature values are used by the control system for estimating the area point spread function .

26. The system of claim 25 wherein the control system estimates a feature value by evaluating a power spectrum of the area blocks, and determines which area blocks have the largest power spectrum.

27. The system of claim 15 wherein the control system estimates a feature value for at least some of the area blocks, and determines which area blocks have the best feature values, and wherein a predetermined number of the area blocks with the best feature values are used by the control system for estimating the area point spread function.

28. The system of claim 15 wherein the control system estimates a feature value for at least some of the area blocks, and estimates the area point spread function by giving more weight in the PSF cost function to area blocks with better feature values and giving less weight to the area blocks with worse feature values.

29. The system of claim 15 further comprising a capturing system that captures an image of a scene.