# Patent

## Title
Methods and systems for image processing

## Abstract

Aspects of the present invention are related to systems and methods for image processing, in particular, de-noising and de-blurring.

## Background
FIELD OF THE INVENTION
Embodiments of the present invention comprise methods and systems for image processing, in particular, noise removal and blur removal.
SUMMARY
Some embodiments of the present invention comprise methods and systems for iteratively estimating a jointly de-noised and de-blurred image from an input image. Exemplary embodiments comprise computation of non-local weights using a current image estimate and generating a new image estimate using the non-local weights, the current estimate and the input image. In some embodiments, a high-resolution image may be estimated from an input image comprising a combination of low-resolution images.
Some embodiments of the present invention comprise methods and systems for image de-noising using non-local weights. In some embodiments, a high-resolution image may be estimated from an input image comprising a combination of low-resolution images.
Some embodiments of the present invention comprise methods and systems for iterative image de-noising using non-local weights. In some embodiments, a high-resolution image may be estimated from an input image comprising a combination of low-resolution images.
Some embodiments of the present invention comprise methods and systems for iterative image de-blurring using local constraints.
The foregoing and other objectives, features, and advantages of the invention will be more readily understood upon consideration of the following detailed description of the invention taken in conjunction with the accompanying drawings.



BRIEF DESCRIPTION OF THE SEVERAL DRAWINGS
FIG. 1 is a picture depicting an exemplary image-capture model;
FIG. 2 is a chart showing exemplary embodiments of the present invention comprising forming a re-sampled image from one or more low-resolution images, generating a weight matrix based on the re-sampled image pixel reliability, and iteratively estimating a de-blurred image;
FIG. 3 is a chart showing exemplary embodiments of the present invention comprising de-blurring using local variations and regularization;
FIG. 4 is a chart showing exemplary embodiments of the present invention comprising de-noising prior to de-blurring;
FIG. 5 is a picture depicting an exemplary image exhibiting spatially repetitive patterns;
FIG. 6 is a picture illustrating a repetitive pattern and the horizontal and vertical period of the repetition;
FIG. 7 is a picture illustrating neighboring pixels associated with a center pixel;
FIG. 8 is a picture illustrating neighboring pixels and associated pixel patches associated with a center pixel;
FIG. 9 is a chart showing exemplary embodiments of the present invention comprising determination of non-local weights and a non-local weighted averaging of an input image;
FIG. 10 is a chart showing exemplary embodiments of the present invention comprising determination of non-local weights which may progressively improve with iterative estimation of a de-noised image;
FIG. 11 is a chart showing exemplary embodiments of the present invention comprising non-local regularization constraints; and
FIG. 12 is a chart showing exemplary embodiments of the present invention comprising non-local regularization constraints and dynamic range expansion.



DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS
Embodiments of the present invention will be best understood by reference to the drawings, wherein like parts are designated by like numerals throughout. The figures listed above are expressly incorporated as part of this detailed description.
It will be readily understood that the components of the present invention, as generally described and illustrated in the figures herein, could be arranged and designed in a wide variety of different configurations. Thus, the following more detailed description of the embodiments of the methods and systems of the present invention is not intended to limit the scope of the invention but it is merely representative of the presently preferred embodiments of the invention.
Elements of embodiments of the present invention may be embodied in hardware, firmware and/or software. While exemplary embodiments revealed herein may only describe one of these forms, it is to be understood that one skilled in the art would be able to effectuate these elements in any of these forms while resting within the scope of the present invention.
Throughout the following description, an image may be denoted using matrix/vector notion, for example, X, or equivalently as a 2-dimensional (2-D) discrete array of pixel values, for example, X (i, j) or X. An operator may be denoted using matrix/vector notation, for example, H, or equivalently as a 2-D operator, for example, a filter, denoted H (i, j). Convolution may be denoted *.
Some embodiments of the present invention comprise methods and systems for constructing a high-resolution image from one or more low-resolution images. A low-resolution image may be captured by a camera or sensor system with a limited, or lower than desired, spatial resolution. A low-resolution image may be subject to one or more types of degradation inherent in the image-capture process by which the low-resolution image may be obtained. Exemplary types of image degradation may comprise blurring, noise and nonlinearities. Blurring may involve suppression of high-frequency detail in the low-resolution image, and blurring may be due to camera optics, the imaging sensor and other factors. Noise may be a form of random signal variation that may be caused by the imaging sensor, by the system electronics, by analog-to-digital conversion and other factors. Exemplary nonlinearities may include saturation and other factors.
FIG. 1 depicts an exemplary capture-process model 10 for an ith low-resolution image capture. An ideal, high-resolution image 12, which may be denoted X, associated with a scene may be degraded 14 by lens and sensor blur, which may be represented by a combined point-spread-function, denoted HPSF. The degraded image 15 may be subject to sensor motion degradation 16, which may be denoted Mi. The resulting degraded image 17 may be down sampled 18, and the down-sampling operator may be denoted Di. Additive noise, which may be denoted Vi, present in the ith captured low-resolution image 22 may be added 20 to the down-sampled image 19. The model illustrated in FIG. 1 may be described by:

Yi=DiMiHPSFX+Vi,

where the ith captured low-resolution image 22 may be denoted Yi.

In some embodiments of the present invention, the down-sampling operation and sensor motion may be known. In these embodiments, the effect of the known down-sampling operation and sensor motion may be undone by combining the observed, low-resolution images into a single, re-sampled image, which may be denoted Yresampled, according to:




  
    
      
        
          
            Y
            resampled
          
          =
          
          ⁢
          
            
              ∑
              
                
                  i
                  =
                  1
                
                ,
                ...
                ⁢
                
                  
                
                ,
                I
              
            
            ⁢
            
              
                M
                i
                
                  -
                  1
                
              
              ⁢
              
                D
                i
                
                  -
                  1
                
              
              ⁢
              
                Y
                i
              
            
          
        
        ,
      
    
  
  
    
      
        
          =
          
          ⁢
          
            
              
                H
                PSF
              
              ⁢
              X
            
            +
            V
          
        
        ,
      
    
  




where V may be the up-sampled, motion-compensated noise field, and Mi−1 and Di−1 may denote the inverse sensor motion and down-sampling operations, respectively, associated with the ith capture. The re-sampled image, Yresampled, may still be degraded by the camera blur, HPSF, and the noise, V.

Some embodiments of the present invention comprise methods and systems for removing the effect of the blur described by HPSP from the re-sampled image, Yresampled. An estimate, which may be denoted {circumflex over (X)}, of the ideal, high-resolution image may be computed by minimization of a cost function formulated as follows:

{circumflex over (X)}=arg min Φ(X),

with

Φ(X)=μYresampled−HPSFX∥2+λC(X),

where C (X) may be a regularizing, also considered stabilizing, constraint.

In some embodiments of the present invention, iterative minimization methods may be used to iteratively compute estimates of X. In some of these embodiments, the steepest descent approach may be used to generate, at each iteration, a next estimate, {circumflex over (X)}k+1, from a current estimate, {circumflex over (X)}k, by updating the current estimate in the opposite direction of the gradient of the cost function. In some embodiments, the initial estimate, {circumflex over (X)}0, may be based directly on the re-sampled image. In some of these embodiments, the initial estimate may be determined according to:

{circumflex over (X)}0=Yresampled.

In some embodiments, when λ=0, the estimate of X may be the unconstrained least-squares (LS) solution. In some of these embodiments, an iterative minimization scheme may be given according to:

{circumflex over (X)}k+1={circumflex over (X)}k−β[HPSFT(HPSF{circumflex over (X)}k−Yresampled)].

where β may be the step-size, or update weight. In some embodiments, the update weight, β, may be held constant. In alternative embodiments, an optimal value of the update weight, β, may be computed at each iteration. In some of these embodiments, the optimal value may be determined to minimize the cost function given the current steepest descent direction. Some embodiments of the present invention may comprise computation of improved descent directions according to a conjugate gradient minimization method. Alternative embodiments may use other minimization methods known in the art.

In some embodiments, when λ>0, the regularization constraint, C ({circumflex over (X)}), may be based on a linear high-pass filter combined with the L2 norm according to the Tikhonov-Miller approach. In some embodiments, the linear high-pass filter may comprise a 2-dimensional (2-D) Laplacian filter, which may be denoted L. In these embodiments, the regularization constraint may be determined according to:

C({circumflex over (X)})=∥L{circumflex over (X)}∥2.

Some embodiments may comprise an iterative minimization according to:

{circumflex over (X)}k+1={circumflex over (X)}k−β[HPSFT(HPSF{circumflex over (X)}k−Yresampled)]−βλLTL{circumflex over (X)}k 

using the steepest descent approach, where β may be the step-size, or update weight. Some embodiments may comprise a discrete Laplacian filter which may comprise weighting of neighboring pixels in a small, local, spatial window.

The regularization constraint may be based on a Total Variation (TV) norm in alternative embodiments of the present invention. In some embodiments, the TV norm may be computed as the L1 norm of the gradient magnitude. In these embodiments, the regularization constraint may be determined according to:

C({circumflex over (X)})=∥Δ{circumflex over (X)}∥1.

where an exemplary discrete version may be determined as follows:





  
    C
    ⁡
    
      (
      X
      )
    
  
  =
  
    
      ∑
      
        i
        ,
        j
      
    
    ⁢
    
      
        D
        
          i
          ,
          j
        
      
      ⁡
      
        (
        X
        )
      
    
  




where (i, j) is a pixel location and





  
    
      
        D
        
          i
          ,
          j
        
      
      ⁡
      
        (
        X
        )
      
    
    ≡
    
      
        
          ∑
          
            l
            ,
            
              m
              ∈
              S
            
          
        
        ⁢
        
          
            [
            
              
                X
                ⁡
                
                  (
                  
                    i
                    ,
                    j
                  
                  )
                
              
              -
              
                X
                ⁡
                
                  (
                  
                    
                      i
                      +
                      l
                    
                    ,
                    
                      j
                      +
                      m
                    
                  
                  )
                
              
            
            ]
          
          2
        
      
    
  
  ,




where S denotes a local neighborhood proximate to the pixel location.

In some embodiments, the local variation, Di,j (X), at a pixel location, (i, j), may be defined using a local neighborhood of the four connected neighboring pixels. In these embodiments, S={(1,0), (0,1), (−1,0), (0,−1)}. In alternative embodiments, the local variation may be defined using a local neighborhood of the eight connected neighboring pixels. In yet alternative embodiments, other local neighborhoods may be used. Using steepest-descent-based minimization, some embodiments may comprise an iterative update according to:




  
    
      
        X
        ^
      
      
        k
        +
        1
      
    
    ⁡
    
      (
      
        i
        ,
        j
      
      )
    
  
  =
  
    
      
        
          X
          ^
        
        k
      
      ⁡
      
        (
        
          i
          ,
          j
        
        )
      
    
    -
    
      β
      ⁡
      
        [
        
          
            
              H
              PSF
              T
            
            ⁡
            
              (
              
                
                  -
                  i
                
                ,
                
                  -
                  j
                
              
              )
            
          
          *
          
            (
            
              
                
                  
                    H
                    PSF
                  
                  ⁡
                  
                    (
                    
                      i
                      ,
                      j
                    
                    )
                  
                
                *
                
                  
                    
                      X
                      ^
                    
                    k
                  
                  ⁡
                  
                    (
                    
                      i
                      ,
                      j
                    
                    )
                  
                
              
              -
              
                
                  Y
                  resampled
                
                ⁡
                
                  (
                  
                    i
                    ,
                    j
                  
                  )
                
              
            
            )
          
        
        ]
      
    
    -
    
      0.5
      ⁢
      β
      ⁢
      
        
      
      ⁢
      
        λ
        ⁡
        
          [
          
            
              ∑
              
                l
                ,
                
                  m
                  ∈
                  S
                
              
            
            ⁢
            
              
                
                  α
                  
                    i
                    ,
                    j
                    ,
                    l
                    ,
                    m
                  
                
                ⁡
                
                  (
                  
                    
                      X
                      ^
                    
                    k
                  
                  )
                
              
              ⁢
              
                (
                
                  
                    
                      
                        X
                        ^
                      
                      k
                    
                    ⁡
                    
                      (
                      
                        i
                        ,
                        j
                      
                      )
                    
                  
                  -
                  
                    
                      
                        X
                        ^
                      
                      k
                    
                    ⁡
                    
                      (
                      
                        
                          i
                          +
                          l
                        
                        ,
                        
                          j
                          +
                          m
                        
                      
                      )
                    
                  
                
                )
              
            
          
          ]
        
      
    
  




where * denotes convolution and





  
    
      α
      
        i
        ,
        j
        ,
        l
        ,
        m
      
    
    ⁡
    
      (
      
        X
        ^
      
      )
    
  
  =
  
    
      1
      
        
          D
          
            i
            ,
            j
          
        
        ⁡
        
          (
          
            X
            ^
          
          )
        
      
    
    +
    
      
        1
        
          
            D
            
              l
              ,
              m
            
          
          ⁡
          
            (
            
              X
              ^
            
            )
          
        
      
      .
    
  




The αi,j,l,m may comprise data-adaptive weights for the local neighbors in S. Alternative embodiments may comprise iterative update according to minimization using optimized steepest descent, conjugate gradient and other minimization methods.

Alternative embodiments of the present invention may comprise methods and systems for removing the effect of the blur described by HPSF from a re-sampled image, Yresampled, formed from multiple low-resolution images captured with multiple different exposure times. In these embodiments, a re-sampled image may be formed according to:




  
    
      Y
      resampled
    
    =
    
      
        ∑
        
          
            i
            =
            1
          
          ,
          ...
          ⁢
          
            
          
          ,
          I
        
      
      ⁢
      
        
          τ
          i
        
        ⁢
        
          M
          i
          
            -
            1
          
        
        ⁢
        
          D
          i
          
            -
            1
          
        
        ⁢
        
          Y
          i
        
      
    
  
  ,




where τi may be based on the exposure time of the ith captured low-resolution image and the maximum exposure time among all images, i=1, . . . , I.

In some embodiments, associated with a re-sampled image, Yresampled, a diagonal weight matrix, which may be denoted W, may be defined, wherein entries W (i, j) on the diagonal may be equal to the reliability of pixel (i, j) in the re-sampled image, Yresampled. In some embodiments, reliability values may be between 0.0 and 1.0, wherein values closer to 0.0 may indicate data samples with lower reliability than those with values closer to 1.0. In alternative embodiments, other values may be used. In some embodiments, pixels with saturated value, or value close to saturation, may be considered unreliable and may be assigned a reliability value close to 0.0.
A cost function comprising dynamic-range expansion may be given as follows for some embodiments:

Φ(X)=(Yresampled−HPSFX)TW(Yresampled−HPSFX)+λC(X).

where W may be the weight matrix comprising the reliability values. Using steepest-descent-based minimization and the TV constraint, some embodiments may comprise an iterative update according to:





  
    
      
        
          X
          ^
        
        
          k
          +
          1
        
      
      ⁡
      
        (
        
          i
          ,
          j
        
        )
      
    
    =
    
      
        
          
            X
            ^
          
          k
        
        ⁡
        
          (
          
            i
            ,
            j
          
          )
        
      
      -
      
        β
        ⁡
        
          [
          
            
              
                H
                PSF
                T
              
              ⁡
              
                (
                
                  
                    -
                    i
                  
                  ,
                  
                    -
                    j
                  
                
                )
              
            
            *
            
              (
              
                
                  W
                  ⁡
                  
                    (
                    
                      i
                      ,
                      j
                    
                    )
                  
                
                ⁢
                
                  (
                  
                    
                      
                        
                          H
                          PSF
                        
                        ⁡
                        
                          (
                          
                            i
                            ,
                            j
                          
                          )
                        
                      
                      *
                      
                        
                          
                            X
                            ^
                          
                          k
                        
                        ⁡
                        
                          (
                          
                            i
                            ,
                            j
                          
                          )
                        
                      
                    
                    -
                    
                      
                        Y
                        resampled
                      
                      ⁡
                      
                        (
                        
                          i
                          ,
                          j
                        
                        )
                      
                    
                  
                  )
                
              
              )
            
          
          ]
        
      
      -
      
        0.5
        ⁢
        β
        ⁢
        
          
        
        ⁢
        
          λ
          ⁡
          
            [
            
              
                ∑
                
                  l
                  ,
                  
                    m
                    ∈
                    S
                  
                
              
              ⁢
              
                
                  
                    α
                    
                      i
                      ,
                      j
                      ,
                      l
                      ,
                      m
                    
                  
                  ⁡
                  
                    (
                    
                      
                        X
                        ^
                      
                      k
                    
                    )
                  
                
                ⁢
                
                  (
                  
                    
                      
                        
                          X
                          ^
                        
                        k
                      
                      ⁡
                      
                        (
                        
                          i
                          ,
                          j
                        
                        )
                      
                    
                    -
                    
                      
                        
                          X
                          ^
                        
                        k
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            l
                          
                          ,
                          
                            j
                            +
                            m
                          
                        
                        )
                      
                    
                  
                  )
                
              
            
            ]
          
        
      
    
  
  ,




where * denotes convolution, and the re-sampled image may be based on the exposure times of the low-resolution images and the cost function may comprise the weight matrix. These embodiments may be described in relation to FIG. 2 and FIG. 3. Multiple low-resolution images may be received 30 and a re-sampled image may be formed 32 from the low-resolution images. The re-sampled image may be formed 32 according to:





  
    Y
    resampled
  
  =
  
    
      ∑
      
        
          i
          =
          1
        
        ,
        ...
        ⁢
        
          
        
        ,
        I
      
    
    ⁢
    
      
    
    ⁢
    
      
        τ
        i
      
      ⁢
      
        M
        i
        
          -
          1
        
      
      ⁢
      
        D
        i
        
          -
          1
        
      
      ⁢
      
        Y
        i
      
    
  




as described above. A weight matrix may be generated 34 based on the re-sampled image, and an initial estimate of the de-blurred image may be formed 36. In some embodiments, the initial estimate may be based on the re-sampled image. An updated estimate may be generated 38 using a TV constraint and minimization procedure. An iteration-termination condition may be examined 40, and if the condition is met 41, the iteration process may terminate 42. If the condition is not met 43, then a new estimate maybe generated 38.

FIG. 3 describes one iteration according to exemplary embodiments of the present invention comprising a TV constraint and reliability weight matrix. The blur operator may be applied 58 using the known point-spread-function 54, HPSF, the re-sampled image 52, Yresampled, the weight matrix 50, W, and the current estimate 56, {circumflex over (X)}k. The local variations may be computed 60 using the current estimate 56, and the regularization operator may be applied 62 using the computed 60 local variations and the current estimate 56. An update direction may be computed 64 using the results of the blur operation and the regularization operation, and the current estimate 56 may be updated generating a new, current estimate 66.
One of the main concerns during de-blurring may be the potential for noise amplification. Noise may be inherently present in captured image data and ignoring noise during de-blurring may render a resulting high-resolution image less useful. The use of regularizing, or stabilizing, constraints during de-blurring may be one way to control or suppress the noise in the data and prevent its amplification. The Total Variation constraint may be effective in controlling the noise. However, this constraint, similar to most regularizing constraints, acts locally and has the potential to re-introduce a small amount of blur. The strength of the constraint may be selected appropriately to balance the dangers of noise amplification and blurring.
Some embodiments of the present invention described in relation to FIG. 4 may comprise, prior to de-blurring 76, de-noising 74 of a re-sampled image 73 formed by combining 72 multiple low-resolution images 70. The de-noised image 75 may be de-blurred 76 to form the final estimate 78 of the high-resolution image. Some embodiments of the present invention may comprise de-noising methods known in the art applied prior to the de-blurring methods and systems above-described herein. Alternative embodiments of the present invention may comprise de-noising according to systems and methods of the present invention described below.
Some embodiments of the present invention comprise methods and systems for de-noising and noise suppression which may allow improved de-noising capability with a decreased potential for re-introducing blur. Some of these embodiments may comprise non-local means.
Some embodiments of the present invention may make use of prior knowledge about repetitive patterns in image data to restrict a non-local means determination. One exemplary scenario wherein the image data exhibits repetitive patterns may be captured images of an LCD (Liquid Crystal Display) panel, where the LCD panel itself may be controlled to display a particular test pattern. The resulting images may contain an array of LCD pixels and sub-pixels. The pattern of LCD pixels visible in an image may be very regular due to tight control over the setup of the camera and LCD panel. An example of such an image 80 is shown in FIG. 5. The image may show a pattern of LCD pixels that may be repeated horizontally and vertically, with a constant period between one pattern and the next, as illustrated in FIG. 6 for the exemplary image 80 shown in FIG. 5. The image 80 comprises a pattern 82 that may repeat with a horizontal period 84, which may be denoted lx, and a vertical period 86, which may be denoted ly. In some embodiments of the present invention, at each pixel in the image, pixels in different areas of the image that are highly similar, and therefore represent excellent candidates for de-noising purposes, may be located.
Some images may comprise small bright or dark irregularities, which may correspond to pixel defects in an LCD panel. These (sub)pixel defects are exceptions in the regular pattern and must be preserved during de-noising. Regular de-noising techniques have a potential to weaken the visibility of these LCD panel defects. Some embodiments of the present invention comprise data-adaptive weights that may avoid this problem.
Embodiments of the present invention may not introduce additional blur, and hence may have an advantage over many existing de-noising techniques. Performing de-noising according to embodiments of the present invention before de-blurring may allow the reduction of the strength of the regularization constraint during de-blurring, since the power of the noise will have been reduced. This may reduce the risk of blur due to regularization, and also may reduce the risk of other artifacts that may occur in some cases when regularization may be applied.
Some embodiments of the present invention comprise non-iterative methods and systems for de-noising an image. Some of these embodiments may comprise weighted averaging using non-local neighbor pixels and data adaptive weights. In these embodiments, an image may be processed pixel-by-pixel, and the pixel considered for processing, for example, the pixel at location (i, j), may be considered the center pixel. A neighbor index set N{ . . . , −3, −2, −1, 0, 1, 2, 3, . . . }, where 0 may refer to the center pixel being processed and other index values refer to neighboring pixels that may be used in the de-noising process, may be defined. A patch of pixels around the center pixel, and likewise around the neighboring pixels, may be defined. The patch size may be defined to be (2A+1)×(2B+1) pixels. Alternative patch sizes may be used. The neighbor index set N may also refer to neighboring patches in addition to neighboring pixels.
The nth neighboring pixel may be offset vertically by Ln pixels and horizontally by Mn pixels in the image. In some embodiments of the present invention, these offsets may be defined using one direction vector. In alternative embodiments of the present invention, these offsets may be defined using two direction vectors. In some embodiments, using one or two direction vectors may depend on whether the pattern repeats along only one direction or along two directions in the image plane. A first direction vector may be defined as r=[rx,ry]T, and a second direction vector may be defined as s=[sx,sy]T. For the exemplary image 80 shown in FIG. 5 and FIG. 6, these direction vectors may be defined as r=[lx,0]T 91 and s=[0,ly]T 93 as illustrated in FIG. 7 relative to a center pixel 90. The offsets for neighboring pixels may be defined as linear integer combinations of the direction vectors according to:

[Mn,Ln]T=pns+qnr 

where pn and qn are integers. Thus, the horizontal and vertical offsets may be determined according to Mn=pnsx+qnrx and Ln=pnsy+qnry, respectively. For example, the eight-nearest neighboring offset pixels may be determined using:

(pn,qn)ε{(−1,−1),(0,−1),(1,−1),(−1,0),(0,0),(1,0),(−1,1),(0,1),(1,1)}.

FIG. 7 illustrates a center pixel 90 and the eight nearest neighboring offset pixels 92, 94, 96, 98, 100, 102, 104, 106 corresponding to (pn,qn)=(0,1), (1,0), (1,1), (−1,0), (−1,1), (1,−1), (0,−1) and (−1,−1), respectively.
In some embodiment of the present invention, neighbors may be restricted to pixels that are located on a regular lattice defined by one or two direction vectors. Only one direction vector may be needed if the image pattern repeats only in one direction.
A de-noising weight for the nth neighboring pixel of the pixel located at (i,j), also considered the center pixel, may be determined, in some embodiments of the present invention, according to:




  
    
      
        w
        
          i
          ,
          j
          ,
          n
        
      
      ⁡
      
        (
        
          Y
          resampled
        
        )
      
    
    =
    
      exp
      ⁢
      
        {
        
          
            -
            
              1
              
                h
                2
              
            
          
          ⁢
          
            
              ∑
              
                u
                =
                
                  -
                  A
                
              
              A
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  
                    -
                    B
                  
                
                B
              
              ⁢
              
                
              
              ⁢
              
                
                  [
                  
                    
                      
                        Y
                        resampled
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                          
                          ,
                          
                            j
                            +
                            v
                          
                        
                        )
                      
                    
                    -
                    
                      
                        Y
                        resampled
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                            +
                            
                              L
                              n
                            
                          
                          ,
                          
                            j
                            +
                            v
                            +
                            
                              M
                              n
                            
                          
                        
                        )
                      
                    
                  
                  ]
                
                2
              
            
          
        
        }
      
    
  
  ,




where h may be a parameter associated with the de-noising strength, wherein a larger value of h may correspond to a higher de-noising strength. In these embodiments, the de-noising weights may be based on the pixel differences between corresponding pixels in the center patch and a neighboring patch. Alternative embodiments of the present invention may use alternative function forms to the exponential function, exp. In some embodiments of the present invention, the weights may be determined using a mathematical formulation. In alternative embodiments of the present invention, a pre-computed mapping, for example, a lookup-table, may be used to determine the weights.

FIG. 8 illustrates an exemplary image 120. An exemplary center pixel 131 and a patch 121 of pixels associated with the center pixel 131 are shown in FIG. 8. Eight neighboring pixels 132-139 and associated patches 122-129 are also shown.
Some embodiments of the present invention may be described in relation to FIG. 9. Non-local weights 154 may be computed 152 from an input image 150. The non-local weights 154 may be used in a non-local weighted averaging 156 of the input image 150 to form a de-noised version 158 of the input image. In some embodiments of the present invention, the input image 150 may comprise a re-sampled image formed from multiple, low-resolution images.
In some embodiments of the present invention, weights may be determined for the neighboring pixels associated with a center pixel, and a de-noised value associated with the center pixel may be determined according to:




  
    
      
        X
        ^
      
      denoised
    
    ⁡
    
      (
      
        i
        ,
        j
      
      )
    
  
  =
  
    
      1
      
        
          ∑
          
            n
            ∈
            N
          
        
        ⁢
        
          
        
        ⁢
        
          
            w
            
              i
              ,
              j
              ,
              n
            
          
          ⁡
          
            (
            
              Y
              resampled
            
            )
          
        
      
    
    ⁢
    
      
        ∑
        
          n
          ∈
          N
        
      
      ⁢
      
        
      
      ⁢
      
        
          
            w
            
              i
              ,
              j
              ,
              n
            
          
          ⁡
          
            (
            
              Y
              resampled
            
            )
          
        
        ⁢
        
          
            
              Y
              resampled
            
            ⁡
            
              (
              
                
                  i
                  +
                  
                    L
                    n
                  
                
                ,
                
                  j
                  +
                  
                    M
                    n
                  
                
              
              )
            
          
          .
        
      
    
  




Some embodiments of the present invention comprise iterative de-noising using non-local neighbor pixels and data-adaptive weights. In these embodiments, the weights may be determined at each iteration using the current estimate of the de-noised image according to:




  
    
      
        w
        
          i
          ,
          j
          ,
          n
        
      
      ⁡
      
        (
        
          
            X
            ^
          
          k
        
        )
      
    
    =
    
      exp
      ⁢
      
        {
        
          
            -
            
              1
              
                h
                2
              
            
          
          ⁢
          
            
              ∑
              
                u
                =
                
                  -
                  A
                
              
              A
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  
                    -
                    B
                  
                
                B
              
              ⁢
              
                
                  [
                  
                    
                      
                        
                          X
                          ^
                        
                        k
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                          
                          ,
                          
                            j
                            +
                            v
                          
                        
                        )
                      
                    
                    -
                    
                      
                        
                          X
                          ^
                        
                        k
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                            +
                            
                              L
                              n
                            
                          
                          ,
                          
                            j
                            +
                            v
                            +
                            
                              M
                              n
                            
                          
                        
                        )
                      
                    
                  
                  ]
                
                2
              
            
          
        
        }
      
    
  
  ,




and an updated de-noised estimate may be generated from the current weight estimates and the re-sampled image according to:





  
    
      
        X
        ^
      
      
        k
        +
        1
      
    
    ⁡
    
      (
      
        i
        ,
        j
      
      )
    
  
  =
  
    
      
        
          X
          ^
        
        k
      
      ⁡
      
        (
        
          i
          ,
          j
        
        )
      
    
    -
    
      
        β
        [
        
          
            1
            
              
                ∑
                
                  n
                  ∈
                  N
                
              
              ⁢
              
                
              
              ⁢
              
                
                  w
                  
                    i
                    ,
                    j
                    ,
                    n
                  
                
                ⁡
                
                  (
                  
                    
                      X
                      ^
                    
                    k
                  
                  )
                
              
            
          
          ⁢
          
            
              ∑
              
                n
                ∈
                N
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  w
                  
                    i
                    ,
                    j
                    ,
                    n
                  
                
                ⁡
                
                  (
                  
                    
                      X
                      ^
                    
                    k
                  
                  )
                
              
              ⁢
              
                
                  Y
                  resampled
                
                ⁡
                
                  (
                  
                    
                      i
                      +
                      
                        L
                        n
                      
                    
                    ,
                    
                      j
                      +
                      
                        M
                        n
                      
                    
                  
                  )
                
              
            
          
        
        ]
      
      .
    
  




In these embodiments, the de-nosing weights may be continuously improved as the iterations progress.
In some embodiments of the present invention, the initial estimate, {circumflex over (X)}0, may be based directly on the re-sampled image. In some of these embodiments, the initial estimate may be determined according to:

{circumflex over (X)}0=Yresampled.

Some embodiments of the present invention may be described in relation to FIG. 10. In these embodiments, non-local weights 164 may be calculated 162 from a re-sampled image 160 to generate a first de-noised image estimate 168 using non-local weighted averaging 166 of the re-sampled image 160. In subsequent iterations, non-local weights 172 may be computed 170 from the current de-noised image estimate 168. The non-local weights 172 based on the current de-noised image estimate and the re-sampled image 160 may be used to form a new, de-noised image estimate 168 using non-local weighted averaging 166 of the re-sampled image 160.
During the first iteration, the de-noising weights may be computed from the input image 160; however, during subsequent iterations, the de-noising weights may be computed on the output of the previous iteration 168. The de-noising weights may be, thus, iteratively improved.
The above-described iterative de-noising scheme may be derived indirectly by introducing a cost function that may be subsequently minimized with a gradient descent method. The cost function may be defined:




  
    
      
        
          
            Γ
            ⁡
            
              (
              
                X
                ^
              
              )
            
          
          =
          
            
              ∑
              
                i
                ,
                j
              
            
            ⁢
            
              
            
            ⁢
            
              
                [
                
                  
                    
                      X
                      ^
                    
                    ⁡
                    
                      (
                      
                        i
                        ,
                        j
                      
                      )
                    
                  
                  -
                  
                    
                      1
                      
                        
                          ∑
                          
                            n
                            ∈
                            N
                          
                        
                        ⁢
                        
                          
                        
                        ⁢
                        
                          
                            w
                            
                              i
                              ,
                              j
                              ,
                              n
                            
                          
                          ⁡
                          
                            (
                            
                              X
                              ^
                            
                            )
                          
                        
                      
                    
                    ⁢
                    
                      
                        ∑
                        
                          n
                          ∈
                          N
                        
                      
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          
                            w
                            
                              i
                              ,
                              j
                              ,
                              n
                            
                          
                          ⁡
                          
                            (
                            
                              X
                              ^
                            
                            )
                          
                        
                        ⁢
                        
                          
                            Y
                            resampled
                          
                          ⁡
                          
                            (
                            
                              
                                i
                                +
                                
                                  L
                                  n
                                
                              
                              ,
                              
                                j
                                +
                                
                                  M
                                  n
                                
                              
                            
                            )
                          
                        
                      
                    
                  
                
                ]
              
              2
            
          
        
        ,
      
    
    
      
        (
        19
        )
      
    
  




with de-noising weights wi,j,n({circumflex over (X)}) defined:





  
    
      w
      
        i
        ,
        j
        ,
        n
      
    
    ⁡
    
      (
      
        X
        ^
      
      )
    
  
  =
  
    exp
    ⁢
    
      
        {
        
          
            -
            
              1
              
                h
                2
              
            
          
          ⁢
          
            
              ∑
              
                u
                =
                
                  -
                  A
                
              
              A
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  
                    -
                    B
                  
                
                B
              
              ⁢
              
                
              
              ⁢
              
                
                  [
                  
                    
                      
                        X
                        ^
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                          
                          ,
                          
                            j
                            +
                            v
                          
                        
                        )
                      
                    
                    -
                    
                      
                        X
                        ^
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                            +
                            
                              L
                              n
                            
                          
                          ,
                          
                            j
                            +
                            v
                            +
                            
                              M
                              n
                            
                          
                        
                        )
                      
                    
                  
                  ]
                
                2
              
            
          
        
        }
      
      .
    
  




An iterative gradient descent algorithm for de-noising may be constructed by determining the gradient, also considered the derivative, of Γ({circumflex over (X)}) with respect to {circumflex over (X)}.
An alternative cost function may be defined:




  
    
      Γ
      ⁡
      
        (
        
          X
          ^
        
        )
      
    
    =
    
      
        ∑
        
          i
          ,
          j
        
      
      ⁢
      
        
      
      ⁢
      
        
          1
          
            
              ∑
              
                n
                ∈
                N
              
            
            ⁢
            
              
            
            ⁢
            
              
                w
                
                  i
                  ,
                  j
                  ,
                  n
                
              
              ⁡
              
                (
                
                  X
                  ^
                
                )
              
            
          
        
        ⁢
        
          
            ∑
            
              n
              ∈
              N
            
          
          ⁢
          
            
          
          ⁢
          
            
              
                
                  w
                  
                    i
                    ,
                    j
                    ,
                    n
                  
                
                ⁡
                
                  (
                  
                    X
                    ^
                  
                  )
                
              
              ⁡
              
                [
                
                  
                    
                      X
                      ^
                    
                    ⁡
                    
                      (
                      
                        i
                        ,
                        j
                      
                      )
                    
                  
                  -
                  
                    
                      Y
                      resampled
                    
                    ⁡
                    
                      (
                      
                        
                          i
                          +
                          
                            L
                            n
                          
                        
                        ,
                        
                          j
                          +
                          
                            M
                            n
                          
                        
                      
                      )
                    
                  
                
                ]
              
            
            2
          
        
      
    
  
  ,




yielding, through gradient determination, an iterative gradient descent algorithm for de-noising.

A de-noised estimate may still contain blur, and may be subsequently de-blurred. In some embodiments, the de-blurring may comprise methods known in the art. In alternative embodiments, the de-blurring may comprise de-blurring according to embodiments of the present invention described herein using cost-function minimization and regularization.
Some embodiments of the present invention comprise methods and systems for removing blur using a non-local regularization constraint based on repetitive patterns in a high-resolution image associated with multiple, low-resolution observations. Some embodiments of the present invention comprising a non-local regularization constraint may provide improved noise suppression. Compared to systems using separate de-noising prior to de-blurring, some embodiments of the present invention comprising a non-local regularization constraint may provide lower computation cost.
An estimate, which may be denoted {circumflex over (X)}, of an ideal, high-resolution image may be computed by minimization of a cost function formulated as follows to achieve simultaneous de-blurring and de-noising:



  
    
      
    
    ⁢
    
      
        
          X
          ^
        
        =
        
          arg
          ⁢
          
            
          
          ⁢
          
            min
            ⁡
            
              [
              
                
                  
                    
                    
                      
                        Y
                        resampled
                      
                      -
                      
                        
                          H
                          PSF
                        
                        ⁢
                        
                          X
                          ^
                        
                      
                    
                    
                  
                  2
                
                +
                
                  λ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    Γ
                    ⁡
                    
                      (
                      
                        X
                        ^
                      
                      )
                    
                  
                
              
              ]
            
          
        
      
      ,
      
        

      
      ⁢
      
        
      
      ⁢
      where
    
  




  
    
      Γ
      ⁡
      
        (
        
          X
          ^
        
        )
      
    
    =
    
      
        ∑
        
          i
          ,
          j
        
      
      ⁢
      
        
          [
          
            
              
                X
                ^
              
              ⁡
              
                (
                
                  i
                  ,
                  j
                
                )
              
            
            -
            
              
                1
                
                  
                    ∑
                    
                      n
                      ∈
                      N
                    
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      w
                      
                        i
                        ,
                        j
                        ,
                        n
                      
                    
                    ⁡
                    
                      (
                      
                        X
                        ^
                      
                      )
                    
                  
                
              
              ⁢
              
                
                  ∑
                  
                    n
                    ∈
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    
                      w
                      
                        i
                        ,
                        j
                        ,
                        n
                      
                    
                    ⁡
                    
                      (
                      
                        X
                        ^
                      
                      )
                    
                  
                  ⁢
                  
                    
                      X
                      ^
                    
                    ⁡
                    
                      (
                      
                        
                          i
                          +
                          
                            L
                            n
                          
                        
                        ,
                        
                          j
                          +
                          
                            M
                            n
                          
                        
                      
                      )
                    
                  
                
              
            
          
          ]
        
        2
      
    
  



and wi,j,n({circumflex over (X)}) may be data-adaptive de-noising weights, based on the similarity between a center-pixel patch and non-local neighboring patches, given by:





  
    
      w
      
        i
        ,
        j
        ,
        n
      
    
    ⁡
    
      (
      
        X
        ^
      
      )
    
  
  =
  
    exp
    ⁢
    
      
        {
        
          
            -
            
              1
              
                h
                2
              
            
          
          ⁢
          
            
              ∑
              
                u
                =
                
                  -
                  A
                
              
              A
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  
                    -
                    B
                  
                
                B
              
              ⁢
              
                
              
              ⁢
              
                
                  [
                  
                    
                      
                        X
                        ^
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                          
                          ,
                          
                            j
                            +
                            v
                          
                        
                        )
                      
                    
                    -
                    
                      
                        X
                        ^
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                            +
                            
                              L
                              n
                            
                          
                          ,
                          
                            j
                            +
                            v
                            +
                            
                              M
                              n
                            
                          
                        
                        )
                      
                    
                  
                  ]
                
                2
              
            
          
        
        }
      
      .
    
  




An iterative update equation for minimizing the cost function may be determined using gradient descent. In some embodiments of the present invention, the iterative update may be made according to:




  
    
      
        X
        ^
      
      
        k
        +
        1
      
    
    ⁡
    
      (
      
        i
        ,
        j
      
      )
    
  
  =
  
    
      
        
          X
          ^
        
        k
      
      ⁡
      
        (
        
          i
          ,
          j
        
        )
      
    
    -
    
      β
      ⁡
      
        [
        
          
            
              H
              PSF
              T
            
            ⁡
            
              (
              
                
                  -
                  i
                
                ,
                
                  -
                  j
                
              
              )
            
          
          *
          
            (
            
              
                
                  
                    H
                    PSF
                  
                  ⁡
                  
                    (
                    
                      i
                      ,
                      j
                    
                    )
                  
                
                *
                
                  
                    
                      X
                      ^
                    
                    k
                  
                  ⁡
                  
                    (
                    
                      i
                      ,
                      j
                    
                    )
                  
                
              
              -
              
                
                  Y
                  resampled
                
                ⁡
                
                  (
                  
                    i
                    ,
                    j
                  
                  )
                
              
            
            )
          
        
        ]
      
    
    -
    
      β
      ⁢
      
        
      
      ⁢
      
        λ
        [
        
          
            
              
                X
                ^
              
              k
            
            ⁡
            
              (
              
                i
                ,
                j
              
              )
            
          
          -
          
            
              1
              
                
                  ∑
                  
                    n
                    ∈
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    w
                    
                      i
                      ,
                      j
                      ,
                      n
                    
                  
                  ⁡
                  
                    (
                    
                      
                        X
                        ^
                      
                      k
                    
                    )
                  
                
              
            
            ⁢
            
              
                ∑
                
                  n
                  ∈
                  N
                
              
              ⁢
              
                
              
              ⁢
              
                
                  
                    w
                    
                      i
                      ,
                      j
                      ,
                      n
                    
                  
                  ⁡
                  
                    (
                    
                      
                        X
                        ^
                      
                      k
                    
                    )
                  
                
                ⁢
                
                  
                    
                      X
                      ^
                    
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        i
                        +
                        
                          L
                          n
                        
                      
                      ,
                      
                        j
                        +
                        
                          M
                          n
                        
                      
                    
                    )
                  
                
              
            
          
        
        ]
      
    
  




where * denotes convolution.

FIG. 11 describes one iteration according to exemplary embodiments of the present invention comprising a non-local regularization constraint based on repetitive patterns in a high-resolution image associated with multiple, low-resolution observations. A blur operator may be applied 188 using a known point-spread-function 184, HPSF, the re-sampled image 182, Yresampled, and the current estimate 186, {circumflex over (X)}k. Non-local weights may be computed 190 using the current estimate 186, and a non-local regularization operator may be applied 192 using the computed 190 non-local weights and the current estimate 186. An update direction may be computed 194 using the result of the blur operation and the regularization operation, and the current estimate 186 may be updated generating a new, current estimate 196.
In some embodiments of the present invention, a re-sampled image, Yresampled, may be formed from multiple low-resolution images captured with multiple different exposure times. In these embodiments, the re-sampled image may be determined according to:




  
    
      Y
      resampled
    
    =
    
      
        ∑
        
          
            i
            =
            1
          
          ,
          ...
          ⁢
          
            
          
          ,
          I
        
      
      ⁢
      
        
      
      ⁢
      
        
          τ
          i
        
        ⁢
        
          M
          i
          
            -
            1
          
        
        ⁢
        
          D
          i
          
            -
            1
          
        
        ⁢
        
          Y
          i
        
      
    
  
  ,




where τi may be based on the exposure time of the ith captured low-resolution image and the maximum exposure time among all images, i=1, . . . , I.

In some embodiments, associated with a re-sampled image, Yresampled, a diagonal weight matrix, which may be denoted W, may be defined, wherein entries W (i, j) on the diagonal may be equal to the reliability of pixel (i, j) in the re-sampled image, Yresampled. In some embodiments, reliability values may be between 0.0 and 1.0, wherein values closer to 0.0 may indicate data samples with lower reliability than those with values closer to 1.0. In alternative embodiments, other values may be used. In some embodiments, pixels with saturated value, or value close to saturation, may be considered unreliable and be assigned a reliability value close to 0.0.
An estimate, which may be denoted {circumflex over (X)}, of an ideal, high-resolution image may be computed by minimization of a cost function formulated to achieve simultaneous de-blurring and de-noising and accounting for dynamic-range expansion may be given as follows for some embodiments:



  
    
      
        X
        ^
      
      =
      
        arg
        ⁢
        
          
        
        ⁢
        
          min
          ⁡
          
            [
            
              
                
                  
                    (
                    
                      
                        Y
                        resampled
                      
                      -
                      
                        
                          H
                          PSF
                        
                        ⁢
                        X
                      
                    
                    )
                  
                  T
                
                ⁢
                
                  W
                  ⁡
                  
                    (
                    
                      
                        Y
                        resampled
                      
                      -
                      
                        
                          H
                          PSF
                        
                        ⁢
                        X
                      
                    
                    )
                  
                
              
              +
              
                λ
                ⁢
                
                  
                
                ⁢
                
                  Γ
                  ⁡
                  
                    (
                    
                      X
                      ^
                    
                    )
                  
                
              
            
            ]
          
        
      
    
    ,
    
      

    
    ⁢
    
      
    
    ⁢
    where
  




  
    
      Γ
      ⁡
      
        (
        
          X
          ^
        
        )
      
    
    =
    
      
        ∑
        
          i
          ,
          j
        
      
      ⁢
      
        
          [
          
            
              
                X
                ^
              
              ⁡
              
                (
                
                  i
                  ,
                  j
                
                )
              
            
            -
            
              
                1
                
                  
                    ∑
                    
                      n
                      ∈
                      N
                    
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      w
                      
                        i
                        ,
                        j
                        ,
                        n
                      
                    
                    ⁡
                    
                      (
                      
                        X
                        ^
                      
                      )
                    
                  
                
              
              ⁢
              
                
                  ∑
                  
                    n
                    ∈
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    
                      w
                      
                        i
                        ,
                        j
                        ,
                        n
                      
                    
                    ⁡
                    
                      (
                      
                        X
                        ^
                      
                      )
                    
                  
                  ⁢
                  
                    
                      X
                      ^
                    
                    ⁡
                    
                      (
                      
                        
                          i
                          +
                          
                            L
                            n
                          
                        
                        ,
                        
                          j
                          +
                          
                            M
                            n
                          
                        
                      
                      )
                    
                  
                
              
            
          
          ]
        
        2
      
    
  



and wi,j,n ({circumflex over (X)}) may be data-adaptive de-noising weights, based on the similarity between a center-pixel patch and non-local neighboring patches, given by:





  
    
      w
      
        i
        ,
        j
        ,
        n
      
    
    ⁡
    
      (
      
        X
        ^
      
      )
    
  
  =
  
    exp
    ⁢
    
      
        {
        
          
            -
            
              1
              
                h
                2
              
            
          
          ⁢
          
            
              ∑
              
                u
                =
                
                  -
                  A
                
              
              A
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  
                    -
                    B
                  
                
                B
              
              ⁢
              
                
              
              ⁢
              
                
                  [
                  
                    
                      
                        X
                        ^
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                          
                          ,
                          
                            j
                            +
                            v
                          
                        
                        )
                      
                    
                    -
                    
                      
                        X
                        ^
                      
                      ⁡
                      
                        (
                        
                          
                            i
                            +
                            u
                            +
                            
                              L
                              n
                            
                          
                          ,
                          
                            j
                            +
                            v
                            +
                            
                              M
                              n
                            
                          
                        
                        )
                      
                    
                  
                  ]
                
                2
              
            
          
        
        }
      
      .
    
  




An iterative update equation for minimizing the cost function may be determined using gradient descent. The iterative update in some embodiments of the present invention may be made according to:




  
    
      
        X
        ^
      
      
        k
        +
        1
      
    
    ⁡
    
      (
      
        i
        ,
        j
      
      )
    
  
  =
  
    
      
        
          X
          ^
        
        k
      
      ⁡
      
        (
        
          i
          ,
          j
        
        )
      
    
    -
    
      β
      ⁡
      
        [
        
          
            
              H
              PSF
              T
            
            ⁡
            
              (
              
                
                  -
                  i
                
                ,
                
                  -
                  j
                
              
              )
            
          
          *
          
            (
            
              
                W
                ⁡
                
                  (
                  
                    i
                    ,
                    j
                  
                  )
                
              
              ⁢
              
                (
                
                  
                    
                      
                        H
                        PSF
                      
                      ⁡
                      
                        (
                        
                          i
                          ,
                          j
                        
                        )
                      
                    
                    *
                    
                      
                        
                          X
                          ^
                        
                        k
                      
                      ⁡
                      
                        (
                        
                          i
                          ,
                          j
                        
                        )
                      
                    
                  
                  -
                  
                    
                      Y
                      resampled
                    
                    ⁡
                    
                      (
                      
                        i
                        ,
                        j
                      
                      )
                    
                  
                
                )
              
            
            )
          
        
        ]
      
    
    -
    
      β
      ⁢
      
        
      
      ⁢
      
        λ
        [
        
          
            
              
                X
                ^
              
              k
            
            ⁡
            
              (
              
                i
                ,
                j
              
              )
            
          
          -
          
            
              1
              
                
                  ∑
                  
                    n
                    ∈
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    w
                    
                      i
                      ,
                      j
                      ,
                      n
                    
                  
                  ⁡
                  
                    (
                    
                      
                        X
                        ^
                      
                      k
                    
                    )
                  
                
              
            
            ⁢
            
              
                ∑
                
                  n
                  ∈
                  N
                
              
              ⁢
              
                
              
              ⁢
              
                
                  
                    w
                    
                      i
                      ,
                      j
                      ,
                      n
                    
                  
                  ⁡
                  
                    (
                    
                      
                        X
                        ^
                      
                      k
                    
                    )
                  
                
                ⁢
                
                  
                    
                      X
                      ^
                    
                    k
                  
                  ⁡
                  
                    (
                    
                      
                        i
                        +
                        
                          L
                          n
                        
                      
                      ,
                      
                        j
                        +
                        
                          M
                          n
                        
                      
                    
                    )
                  
                
              
            
          
        
        ]
      
    
  




where * denotes convolution and W (i, j) denotes the reliability weight.

FIG. 12 describes one iteration according to some embodiments of the present invention comprising a non-local regularization constraint based on repetitive patterns in a high-resolution image associated with multiple, low-resolution observations comprising differing exposure times. A blur operator may be applied 208 using a known point-spread-function 204, HPSF, the re-sampled image 202, Yresampled, the reliability weights 200, W, and the current estimate 206, {circumflex over (X)}k. Non-local weights may be computed 210 using the current estimate 206, and a non-local regularization operator may be applied 212 using the computed 210 non-local weights and the current estimate 206. An update direction may be computed 214 using the result of the blur operation and the regularization operation, and the current estimate 206 may be updated generating a new, current estimate 216.
The terms and expressions which have been employed in the foregoing specification are used therein as terms of description and not of limitation, and there is no intention in the use of such terms and expressions of excluding equivalence of the features shown and described or portions thereof, it being recognized that the scope of the invention is defined and limited only by the claims which follow.

## Claims

1. A method for forming an image, said method comprising:
associated with a first image, determining a pattern-repetition direction vector comprising a horizontal component and a vertical component;
for a first pixel in a reconstructed image, determining a first plurality of weights associated with a first plurality of pixels in said reconstructed image, wherein
each pixel in said first plurality of pixels is located an associated integral number of horizontal-component pixel values from said first pixel in the horizontal direction and an associated number of vertical-component pixel values from said first pixel in the vertical direction;
each pixel in said first plurality of pixels comprises an associated pixel value; and
said weights are determined using said reconstructed image; and

updating a first pixel value associated with said first pixel in said reconstructed image based on the value of a pixel in said first image at a location associated with said first pixel in said reconstructed image and a combination of said associated pixel values of said first plurality of pixels using said first plurality of weights, wherein said determining a pattern-repetition direction vector, said determining a first plurality of weights and said updating are performed by a computing device.

2. A method as described in claim 1, wherein:
said first image is produced by combining a second image and a third image; and
said first image has a first-image spatial resolution greater than a second-image spatial resolution associated with said second image and a third-image spatial resolution associated with said third image.

3. A method as described in claim 2, wherein said second-image spatial resolution and said third-image spatial resolution are approximately equal.

4. A method as described in claim 2, wherein said combining a second image and a third image comprises:
a first inverse motion operator associated with said second image;
a second inverse motion operator associated with said third image;
a first inverse down-sampling operator associated with said second image; and
a second inverse down-sampling operator associated with said third image.

5. A method as described in claim 2 further comprising:
associating a reliability value with the value of said pixel in said first image at said location associated with said first pixel location in said reconstructed image; and
wherein:
said combining a second image and a third image is based on a first exposure time associated with said second image and a second exposure time associated with said third image; and
said updating is further based on said reliability value.

6. A method as described in claim 1 further comprising initializing said reconstructed image to said first image.

7. A method as described in claim 1, wherein said determining a first plurality of weights comprises, for a second pixel in said first plurality of pixels and said first pixel, calculating the pixel-value differences between corresponding pixels in a first neighborhood proximate to said first pixel and a corresponding second neighborhood proximate to said second pixel.

8. A method as described in claim 7 further comprising:
squaring said differences;
summing said squared differences;
adjusting said summed, squared difference according to a first parameter; and
exponentiating the adjusted, summed, squared difference.

9. A method for forming an image, said method comprising:
producing a combined image by combining a first image and a second image using a first exposure time associated with said first image and a second exposure time associated with said second image, wherein said combined image has a combined-image spatial resolution greater than a first spatial resolution associated with said first image and a second spatial resolution associated with said second image;
associating a reliability value with a first pixel in said combined image; and
updating a first pixel value in a reconstructed image at a location corresponding to said first pixel in said combined image based on said reliability value, a value associated with said first pixel in said combined image and a local variation measure of said reconstructed image in a region proximate to said location, wherein said combining, said associating and said updating are performed by a computing device.

10. A method as described in claim 9, wherein said first spatial resolution and said second spatial resolution are approximately equal.

11. A method as described in claim 9, wherein said region comprises a four-connected neighborhood.

12. A method as described in claim 9, wherein said combining comprises:
a first inverse motion operator associated with said first image;
a second inverse motion operator associated with said second image;
a first inverse down-sampling operator associated with said first image; and
a second inverse down-sampling operator associated with said second image.

13. A method as described in claim 9, wherein said updating comprises minimization of a first cost function comprising a first constraint based on a total variation norm.

14. A method for forming an image, said method comprising:
associated with a first image, determining a pattern-repetition direction vector comprising a horizontal component and a vertical component;
for a first pixel in said first image, determining a first plurality of weights associated with a first plurality of pixels in said first image, wherein each pixel in said first plurality of pixels is located an associated integral number of horizontal-component pixel values from said first pixel in the horizontal direction and an associated number of vertical-component pixel values from said first pixel in the vertical direction and comprises an associated pixel value; and
determining a first de-noised pixel value in a de-noised image at a location associated with said first pixel in said first image by combining said associated pixel values of said first plurality of pixels using said first plurality of weights, wherein said determining a pattern-repetition direction vector, said determining a first plurality of weights and said determining a first de-noised pixel value are performed by a computing device.

15. A method as described in claim 14, wherein:
said first image is produced by combining a second image and a third image; and
said first image has a first-image spatial resolution greater than a second-image spatial resolution associated with said second image and a third-image spatial resolution associated with said third image.

16. A method as described in claim 15, wherein said combining a second image and a third image comprises using a first exposure time associated with said second image and a second exposure time associated with said third image.

17. A method as described in claim 15, wherein said second-image spatial resolution and said third-image spatial resolution are approximately equal.

18. A method as described in claim 15, wherein said combining comprises:
a first inverse motion operator associated with said second image;
a second inverse motion operator associated with said third image;
a first inverse down-sampling operator associated with said second image; and
a second inverse down-sampling operator associated with said third image.

19. A method as described in claim 14, wherein said determining a first plurality of weights comprises, for a second pixel in said first plurality of pixels and said first pixel, calculating the pixel-value differences between corresponding pixels in a first neighborhood proximate to said first pixel and a corresponding second neighborhood proximate to said second pixel.

20. A method as described in claim 19 further comprising:
squaring said differences;
summing said squared differences;
adjusting said summed, squared difference according to a first parameter; and
exponentiating the adjusted, summed, squared difference.

21. A method as described in claim 14, further comprising:
associating a reliability value with said first de-noised pixel; and
updating a first pixel value in a reconstructed image at a location corresponding to said first de-noised pixel in said de-noised image based on said reliability value, a value associated with said first de-noised pixel in said de-noised image and a local variation measure of said reconstructed image in a region proximate to said location corresponding to said first de-noised pixel.

22. A method for forming an image, said method comprising:
associated with a first image, determining a pattern-repetition direction vector comprising a horizontal component and a vertical component;
for a first pixel in said first image, determining a first plurality of weights associated with a first plurality of pixels in said first image, wherein
each pixel in said first plurality of pixels is located an associated integral number of horizontal-component pixel values from said first pixel in the horizontal direction and an associated number of vertical-component pixel values from said first pixel in the vertical direction;
each pixel in said first plurality of pixels comprises an associated pixel value; and
said weights are determined using a de-noised image associated with said first image; and

updating a first de-noised pixel value at a location in said de-noised image corresponding to said first pixel in said first image by combining said associated pixel values of said first plurality of pixels using said first plurality of weights, wherein said determining a pattern-repetition direction vector, said determining a first plurality of weights and said updating are performed by a computing device.

23. A method as described in claim 22 further comprising initializing said de-noised image to said first image.

24. A method as described in claim 22, wherein said determining a first plurality of weights comprises, for second pixel in said de-noised image corresponding to a second pixel location in said first image from said first plurality of pixels and for a first pixel location in said de-noised image corresponding to said first pixel in said first image, calculating the pixel-value differences between corresponding pixels in a first neighborhood proximate to said first pixel location in said de-noised image and a corresponding second neighborhood proximate to said second pixel location in said de-noised image.

25. A method as described in claim 24 further comprising:
squaring said differences;
summing said squared differences;
adjusting said summed, squared difference according to a first parameter; and
exponentiating the adjusted, summed, squared difference.

26. A method as described in claim 22, wherein:
said first image is produced by combining a second image and a third image; and
said first image has a first-image spatial resolution greater than a second-image spatial resolution associated with said second image and a third-image spatial resolution associated with said third image.

27. A method as described in claim 26, wherein said combining a second image and a third image comprises using a first exposure time associated with said second image and a second exposure time associated with said third image.

28. A method as described in claim 26, wherein said second-image spatial resolution and said third-image spatial resolution are approximately equal.

29. A method as described in claim 26, wherein said combining comprises:
a first inverse motion operator associated with said second image;
a second inverse motion operator associated with said third image;
a first inverse down-sampling operator associated with said second image; and
a second inverse down-sampling operator associated with said third image.

30. A method as described in claim 22, further comprising:
associating a reliability value with said first de-noised pixel; and
updating a first pixel value in a reconstructed image at a location corresponding to said first de-noised pixel in said de-noised image based on said reliability value, a value associated with said first de-noised pixel in said de-noised image and a local variation measure of said reconstructed image in a region proximate to said location.