A Tutorial on Image Restoration

Peyman Milanfar* Electrical Engineering Department and CfAO UCSC August, 2003

* Based on slides from Michael Elad

1

Outline 1. The image restoration problem 2. Elementary solutions 3. Statistical estimators 4. Space Varying Restoration 5. Robust Estimation 6. Non-Linear Restoration 7. What Next 2

The Forward Model Original image - X Blur H

0.0144 0.0281 0.0351 0.0281 0.0144

0.0281 0.0547 0.0683 0.0547 0.0281

0.0351 0.0683 0.0853 0.0683 0.0351

0.0281 0.0547 0.0683 0.0547 0.0281

Measured image - Y

0.0144 0.0281 0.0351 0.0281 0.0144

Noise - N

X,Y,N of size M × M 3

1.2 A Section 300

X

250

200

150

Y

100

50

0

0

20

40

60

80

100

120

140

160

180

200

4

1.3 Frequency Domain Low-pass filter

F

=

=Hf(u,v)

F

=

=Xf (u,v)

F

=

= Nf (u,v)

F

=

= Yf (u,v)

Yf=HfXf+Nf Point-by-point multiply 5

1.8 Importance • The image restoration problem is common in imaging systems. • Other important problems can be solved using the same tools developed here. • The problem and its solutions form a very interesting mathematical field, called “Inverse Problems”. 6

Chapter 2 Elementary Solution 7

2.1 Basic Idea Yf=HfXf+Nf Neglect the noise

H f−1Yf = X f + H f−1 N f Element-by-Element inverse

Restored

Hf

Inverse filter 8

2.2 Avoiding Singularities Note: If Hf has zero elements, then at those frequencies, the inverse filter does not exist. Thus use the following inverting equation instead:

H  ≈ 1 * * Hf Hf +C  Hf C H

* f

−1 f

for for

2

H f >> C 2

H f >1. This is silly, since we amplify mostly noise! 13

2.7 Change of Strategy Previous Strategy - Based on Hf: Where Hf is high, apply exact inverse, and where low, apply Hf/C New Startegy

- Based on [Signal/Noise]*Hf : |Hf |SNR>>1 - inverse |Hf |SNR f f  =  SNR2 2 * 2  C Hf Hf ⋅ SNR T

-T

T

• For measurements falling in [-T,T], the behavior is as in the classic ρ(x)=x2 case • Outliers (|x|>T), do not impact the estimation at all! 49

5.5 Summary • In the robust estimation paradigm, we are not restricting ourselves to ρ(x)=x2, but rather using other ρ(x) functions • By using less-convex functions, we get better treatment for outliers • The problem with the resulting approach is that it might yield non-convex optimization problems 50

Chapter 6 Robust Smoothness 51

6.1 The Basic Idea Instead of using

ε ε

2 MAP

2 MAP

(X ) =

(X ) =

2

Y − H X + λ DX

2

Let us use

Y − H X + λρ{DX} 2

where ρ(DX) stands for applying ρ(x) per each entry and summing all the results 52

6.2 A simple 1D Example Assume that for a signal x we know that

x 0 = 0, x10 = 1 and we would like to find x2, x3, … , x9, such that the following penalty is minimized: 10

ε 2 (X ) = ∑ ρ{x k − x k −1} k =1

53

6.3 The Problem 1.0

?

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

0

1

2

3

4

5

6

7

8

9

10

54

6.4 The Example’s Results

10 jumps of 0.1

4 jumps of 0.25

1 jump of 1

0.10

1.00

3.16

0.25

1.00

2.00

1.00

1.00

1.00 55

6.5 Robust Smoothness ε

2 MAP

(X ) =

Y − H X + λρ{DX} 2

Looking at the image DX, small values correspond to smooth regions, and thus should be treated as before. High values correspond to edges. Instead of suppressing edges, we reduce the cost for these points by using an appropriate ρ(x) function 56

6.6 Iterative Solution ε ∂ε

2 MAP

2 MAP

(X ) =

Y − H X + λρ{DX} 2

(X ) = 2H T [H X − Y ] + 2λDTρ' {DX}

∂X

{

}

X k +1 = X k − µ H [H X k − Y ] − λD ρ' {DX k } T

T

X0 can be chosen to be Y 57

6.7 Iterative System Xk

Y HT

H D

ρ’

DT

λ

µ Xk+1 58

6.8 Effective Weight Weight approach X k + 1 = ..... − λ D T WD X k T Robust approach X k +1 = ..... − λD ρ' {D X k }

In order to get that they are equivalent we should require ρ' {D X k } = WD X k

ρ' (x ) w (x ) = x

which means that we have to re-compute the weights at every iteration 59

6.9 Robust Functions ρ( x )

ρ' ( x )

0.5x 2 x

 0.5x 2  0.5T 2

ρ' ( x ) / x

x sign( x )

x ≤T

x x ≤T  otherwise  0 otherwise x

sign( x ) x 60

6.11 Example Using ρ(x)=|x| and X0=Y, Error=27.6 700 600 500 400

ε (X )

300

2

200 100

ˆ MAP X−X

2

60

80

0 0

20

40

100

61

6.10 Local Minimum For NON-CONVEX ρ(x)

ε

2 MAP

(X ) =

Y − H X + λρ{DX} 2

is not convex as well! The solution may be different for different initializations 62

Chapter 7 What Next ? 63

7.1 The Chosen Prior Definition: A ‘prior’ is a function f(X) that gives very low values for ‘good’ images and high values for ‘bad’ images

We have seen so far:

1. f (X) = DX

2

2. f (X) = [DX] W[DX] 3. f (X) = ρ{DX} T

64

7.2 Better Priors Priors should exploit properties of images in each application: e.g. 1. Presence of point sources, 2. Multi-resolution behavior, 3. Edges directionality, 4. Relation between colors 5. If possible, higher level description 65

7.3 Learning Priors Instead of GUESSING how images should look, why not let them tell their story. Take a large set of images and learn from them what is the best prior (or at least, what are the parameters for a pre-specified parametric prior) Very well suited to astronomical applications.

66

7.4 More General Problems 1. Image Scaling problem 2. Super-Resolution problem 3. Image sequence restoration problem 4. Other imaging degradation models for blur and noise 5. Other inverse problems having the same structure 67

Topic : Multiframe Resolution Enhancement

68

Resolution Enhancement Idea • Given multiple low-resolution moving images of a scene (a video), generate a high resolution image (or video). frame1 , frame2 , L, frameN-1 , frame N 144444424444443 High Resolution Frame

frame1 , frame2 , L, L , frame N -1 , frame N 144 42444 3 144424443 High Resolution Frame1

High Resolution Frame2

“Trading off time resolution or view diversity to gain spatial resolution” 70

Resolution Enhancement Model • A simple model relating the low-resolution blurry image to the high resolution crisper image. “PSF”

x1 x3

y1

x2

y1 = a1 x1 + a 2 x 2 + a3 x3 + a 4 x 4 + e1

x4

y2 = 0 ⋅ x1 + a1 x2 + 0 ⋅ x3 + a3 x4 + e2 y3 = 0 ⋅ x1 + 0 ⋅ x2 + a1 x3 + a2 x4 + e3 y4 = 0 ⋅ x1 + 0 ⋅ x2 + 0 ⋅ x3 + a1 x4 + e4 71

Low vs High Res Pixels

x2 enhancement

72

The Mathematical Model y k = Ak X + e k

k-th frame

for

1≤ k ≤ N

Ak = Dk H k Fk Downsampling

Warping

Blurring

Α k = [Τ k,1

Τ k, 2

L

Τ k ,N 2 ] Upper-banded, “nearly” Toeplitz

BTTB system*

Y = AX + e

*PCG-based methods developed by N. Nguyen, 2000

73

The Warping Operation as a Matrix (F) Z

X Per every point in X find a matching point in Z

 x1    M      xj  =      M    x N  

1

0

1

F[j,i]=1

0

  z1   M     z j     M  1  z N  74

Additional Intricacies • • • • •

Sampling is not a point operation – there is a blur Motion (warping) must be estimated! Motion may include perspective warp, local motion, etc. The system is typically underdetermined and ill-conditioned. 10’s or 100’s of thousands of unknown variables and data. 75

Some Examples: License Plate Reading

Digital Video Camera from 2nd story window 81

More Examples

Infrared Camera (Night Vision) Data Courtesy Wright Labs USAF

82

83 Data Courtesy Vigilant Technology

Before

84

After

85

Detail Before

86 Data Courtesy Vigilant Technology

Detail After

87 Data Courtesy Vigilant Technology

Conclusions • The field of image restoration has been around for more than 40 years, with more than 10,000 papers, and numerous active researchers in engineering and elsewhere. • There is much more to say because • New applications continue to pop up, • The results are still not satisfactory, • As computers improve, complicated solutions become practical

88